Chromokinetics of metabolic pathways
Jo¨ rg W. Stucki
Department of Pharmacology, University of Bern, Bern, Switzerland
Some methods to study and intuitively understand steady-
state flows in complicated metabolic pathways are dis-
cussed. For this purpose, a suitable decomposition of
complex metabolic schemes into smaller subsystems is used.
These independent subsystems are then interpreted as basic
colors of a chromatic coloring scheme. The mixture of these
basic colors allows an intuitive picture of how a steady state
in a metabolic pathway can be understood. Furthermore,
actions of drugs can be more easily investigated on this
basis. An anaerobic variant of pyruvate metabolism in rat
liver mitochondria is presented as a simple example. This
experiment allows measurement of the percentage that each
basic color contributes to the steady states resulting from
different experimental conditions. Possible implementations
of existing algorithms and rational design of new drugs
are discussed. A
MATHEMATICA
program, based on a
new algorithm for finding all basic colors of stoichiometric
networks, is included.
Keywords: stoichiometric network analysis; extremal cur-
rents; elementary flux modes; steady-state flux analysis;
rational drug design.
When confronted with a realistic metabolic pathway, e.g.
glycolysis and fatty acid breakdown supplemented with the
citric acid cycle plus oxidative phosphorylation, most
investigators rapidly lose track because of the sheer com-
plexity of such schemes. Therefore, a way is sought to split

quantitative analysis from a computer simulation of the
differential equations, and owing to the close analogy with
coloring methods, I call this analysis ÔchromokineticsÕ.Note,
however, that this analysis applies to steady states only and
cannot cope with transients. I would like to carry over the
intuitive picture of coloring to the reader without insisting
too much on the powerful mathematical machinery stand-
ing behind it. I suspect that with the help of this
chromokinetic idea, metabolic pathways can be understood
more intuitively and predictions can also be made. The aim
of this article is to illustrate this interpretation, step-by-step,
with a simple metabolic example.
The concept behind these basic colors is not new. It came
already to kinetic schemes in different guises: Clarke called
them Ôextremal currentsÕ [1,2]; later, Schuster called some-
thing similar Ôelementary modesÕ [3,4]. In the present report,
I stick to the original terminology introduced by Clarke [1].
The results and definitions published by Clarke are precise,
exhaustive and nonintuitive. Here I just add some color to
these abstract pictures.
A simple example
Anaerobic mitochondrial pyruvate metabolism
First, a simple example is needed. Isolating small autonomic
metabolic units from a whole organism can bring about
simplification. One example of this is mitochondria iso-
lated from rat liver. However, mitochondrial metabolism
still contains far too many reactions to allow coherent
illustration of the procedure. Hence, further simplification is
needed, and this can be achieved by external constraints
Correspondence to J. W. Stucki, Department of Pharmacology,

freedom. Each of these can be conceived of as running
separately, like a little independent clockwork toy. This
requires that we have to take care of conserved moieties. For
example, NADH and NAD
+
cannot pass the inner
mitochondrial membrane and are present only in minute
amounts in the mitochondria. Therefore, we must make
sure that the NADH produced can readily be consumed by
an appropriate oxidizing reaction, thus reproducing
NAD
+
. If not, the clockwork toy would come to a rapid
stop and no measurable net flow would result. The same
applies for the conserved pool of metabolite moieties linked
with CoA plus free CoASH. Figure 2 shows the unique four
partial diagrams, which the scheme admits as decomposi-
tions of the grand view. The calculation of these partial
diagrams will be shown below.
These diagrams are the basic colors used for the mixture
of any steady state. By defining the percentage that each of
these colors contributes, or weights, we can unambiguously
characterize the resulting steady state. To fix the chromatic
idea, we assign, to each partial diagram, a color, viz.: j
1
,
blue; j
2
,red;j
3

chromatic scheme by studying the action of a drug. In the
present scheme, we used fluorocitrate for inhibiting aconi-
tase. This substance is not a particularly useful drug, but is,
in fact, a deadly poison. Recalling that our ischemic
Fig. 1. Simplified model of mitochondrial pyruvate metabolism. This reaction scheme is a simplified version of a previously published model of
mitochondrial pyruvate metabolism, where computer-simulated and measured fluxes were compared [14]. The transport of CO
2
,ATP,ADPandP
i
through the inner mitochondrial membrane was omitted in this scheme. In other words, these molecules were treated as (constant) external
reactants (shown in italics). Furthermore, only the entry of citrate into the Krebs cycle is included and further metabolic changes were neglected, as
supported by the experimental results in Table 1. The r
i
at the arrows is used for the numbering scheme in Scheme 1. The elementary reactions are
listedinAppendix1;seealsoÔConcluding remarksÕ.
2746 J. W. Stucki (Eur. J. Biochem. 271) Ó FEBS 2004
mitochondria are already in a state which can be considered
dead for all practical purposes, because oxidative phos-
phorylation is blocked, still further inhibition of vital
reactions will cause no significant harm, but will be very
useful for illustrating the analysis. From Fig. 2 we can
readily deduce the balance equations for the contribution
(or fixed weight owing to stoichiometry) of each basic color,
j
i
, to the overall consumption and production of the
metabolites in the incubation medium:
Pyruvate ¼ 4j
1
þ 3j

square solution of an overdetermined system. The results,
j
i
, calculated from our experiment, are also included in
Table 1.
The result represented in this form is stunning indeed. To
appreciate this, suppose that one had to describe the
outcome of the incubation with standard practice, without
having access to chromatic decomposition methods. One
would then write, for example, that pyruvate consumption
declines monotonically with increasing fluorocitrate con-
centration, whereas production of 3-hydroxy-butyrate first
increases, then decreases, etc. The final statement would
then somehow say that further experiments, probably based
on radioactive tracer methods, are needed to better under-
stand what is going on.
By contrast, we claim that we can already explain what is
going on, on the basis of the meagre data record in Table 1,
Fig. 2. Independent currents in anaerobic pyruvate metabolism. These diagrams are graphical representations of the independent currents into which
thereactionschemeinFig.1canbedecomposed.Theycorrespondtothebasiccolorsj
1
–j
4
and were calculated as described in the text (Scheme 2).
The numbers on the arrows are the weights of the j
i
on the different metabolites involved. Thus, for example, the pyruvate metabolized in (A) is
4 · j
1
; malate produced 2 · j

Pyruvate (used) 12.56 7.92 2.00
Acetoacetate
0.34 0.23 0.00
3-Hydroxybutyrate
1.31 1.50 0.73
Malate
5.57 2.45 0.66
Citrate
0.67 0.69 0.01
Color (l moles/20 min)
j
1
(blue)
0.34 (9.41%) 0.23 (9.27%) 0.00 (0%)
j
2
(red)
1.23 (34.07%) 1.53 (61.69%) 0.72 (100%)
j
3
(green)
0.59 (16.34%) 0.72 (29.03%) 0.00 (0%)
j
4
(white)
1.45 (40.16%) )0.0 (0.0%) )0.0 (0%)
Ó FEBS 2004 Chromokinetics of metabolic pathways (Eur. J. Biochem. 271) 2747
even without sophisticated recourse to radioactive tracer
techniques. The coloring scheme of steady-state kinetics
allows us to clearly see what has happened. At the low

of these experimental results. Most fortunately our scheme
has only four basic colors, and thus 3D space just allows
translation of the chromatic scheme into a simple platonic
solid: a tetrahedron. In Clarke’s terminology this would be
called the convex current polytope resulting from cutting
the current cone with the hyper plane.
Figure 3A shows this current polytope in 3D space. If we
happened to have more than four colors, say five, we would
need four dimensions for a drawing, which could then no
longer be visualized. Of course, our example has been
tailored such that it can be represented in 3D. But it is
important to stress that we deal here with a constructed
special case and that, in general, the current polytopes
cannot be visualized (see below in ÔConcluding remarksÕ).
In our experiments, we have three different steady states:
that of the control and two with different concentrations of
fluorocitrate. Each one of these steady states must be
located at a defined point within the strict interior of the
current polytope. As three points define a plane, we can thus
construct such an Ôexperimental planeÕ and cut it with the
tetrahedron, and also calculate its corresponding color by
mixing the basic colors involved. This is illustrated in
Fig. 3B. Of course, this plane is an artificial construct, as by
selecting a more detailed range of fluorocitrate concentra-
tions, the steady states would probably no longer all lie
within that plane but rather follow a curve off the
Ôexperimental planeÕ. However, this simplification is useful
for illustrating the basic idea. In Fig. 3C, the Ôexperimental
planeÕ in 3D is then finally rotated into the 2D plane. It is
apparent that each of the three steady states has a color,

makes this decomposition a complex mathematical problem.
Fig. 3. Current polytope and experimental steady states are represented as colored geometrical objects. The current polytope of anaerobic mito-
chondrial pyruvate metabolism is represented as a tetrahedron. For the coloring of the faces in (A) a corresponding RGB value was associated with
the spatial coordinates. In (B) the Ôexperimental planeÕ was calculated using the percentage contributions of the three experimental conditions in
Table 1. This plane was then colored in place within the tetrahedron again by assigning a corresponding RGB value furnished by the spatial
coordinates. In (C), the same triangular plane was projected and colored into the 2D plane. The experimental points were added as yellow circles of
increasing size corresponding to the control, 0.33 m
M
and 0.66 m
M
fluorocitrate successively. This geometrical representation clearly demonstrates
how all of the three steady states are located within the accessible region of the current polytope (see Concluding remarks).
2748 J. W. Stucki (Eur. J. Biochem. 271) Ó FEBS 2004
Combinatorial method
Confronted with the chromatic decomposition problem
for the first time, the immediate reaction is to attempt
an exhaustive trial-and-error approach. Via a rigorous
enumeration, one might then try out all possible combi-
nations of arrows in the reaction scheme and sort out the
ones allowing a steady state. Clarke has published such a
brute force method in a key publication [1]. The algorithm
described therein only needs to be supplemented by a
small routine which generates all the combinations using a
backtrack algorithm [5]. I will not reproduce the pro-
cedures here in detail, as they will be of little use for
realistic metabolic pathways. The major problem with this
brute force approach is that it will rapidly explode. This
has to do with the fact that the number of combinations
grows faster than exponentially with growing network
size.

the following statistics: four basic colors (extremal currents
with positive components only), which were repeatedly
found 15 times; 27 singular submatrices; and 174 sub-
matrices that contained positive as well as negative com-
ponents. It is important to remain firm and to insist on
non-negative components only. Trying to interpret negative
components as backward reactions will rapidly lead to a
quagmire of sheer confusion. Thus, we strictly followed the
practice, decreed by Clarke, of formulating backward
reactions as separate reactions (which in our scheme we
actually did not need to do).
Summarizing the results of the calculations, we obtain
a matrix containing the j
i
values of the basic colors in
terms of the reactions, r
i
, as columns (Scheme 2). The
rows j
1
to j
4
of this color matrix correspond to the partial
diagrams in Fig. 2. The numbers above the arrows in the
partial diagrams in Fig. 2 are the weights given in the
colors matrix, whereas the experimental data in Table 1
yield the percentage contributions or barycentric coordi-
nates of each basic color, j
i
. Clarke calls this matrix ÔEÕ

The observation, that out of the many combinations only
an exceedingly small minority of solutions are useful basic
colors, inspires a bottom-up instead of a top-down
approach, as used above. The advantage of this has to do
with the basic organization of metabolic pathways, which
seem to resemble a collection of hubs rather than a
completely connected system.
I have decided to be brief on using mathematical
notations in this report. However, sometimes the use of
mathematical notations is unavoidable, lest the impression
of a lukewarm presentation is given. Thus, we express, for
the time being, changes of the concentrations as:
dc/dt ¼ f(v), (Eqn 2)
The reaction velocities, v
i
, may be nonlinear functions of the
concentrations (or of the chemical activities) of the species
involved. Delegating nonlinearity into the v
i
values, in this
manner, we note that the concentrations, c
i
, are linear
functions of these velocities and we put:
Ó FEBS 2004 Chromokinetics of metabolic pathways (Eur. J. Biochem. 271) 2749
dc/dt ¼ Nv, (Eqn 3)
where N is the stoichiometry matrix, as explained above.
The very definition of the steady state declares that the
concentration changes, dc/dt, must vanish as time proceeds
to infinity. This allows:

the Euclidean space spanned by the reaction velocities;
and, second, only one cone vector can be located within a
plane. All other vectors are not basic colors, but linear
combinations thereof. On the basis of these geometrical
properties, Wagner found a new algorithm for the
calculation of the current polytope from the kernel of
the stoichiometry matrix [6].
The
MATHEMATICA
program given in Appendix 2 imple-
ments this algorithm. Steps (1) and (2) have been already
mentioned above. Step (3) copies the kernel ÔkÕ into an initial
tableau and makes a bitmap with entries True for zeroes in
the tableau and False elsewhere. Step (4) sets up the filter
functions containing the necessary tests for rejection of
unfeasible solutions. In (5) the tableau is processed column-
wise, such that zeroes are produced by all possible linear
combinations of suitable rows of the tableau. If the resulting
combinations are neither identical nor themselves a linear
combination of already existing rows of the tableau, they
are attached to the tableau sequentially. The necessary test
is based on the geometrical properties mentioned above.
Finally, step (6) displays all basic colors of the network.
They correspond to the ones already given in Eqn (3).
Note that the above example is trivial insofar as only
pairwise combinations were needed. Other models may
require much more rich orchestration, such as trios,
quartets, etc., of combinations. An instructive example is
the model of the Belusov–Zhabotinsky reaction, discussed
by Clarke [2], which needs up to quintets. In addition, still-

systems as complicated as the metabolism of Escherichia coli
bacteria. This algorithm is based on a general solution of a
linear programming problem. A simplification thereof was
worked out in detail by Schuster and collaborators [3,4]. By
contrast to the solution based on the nullspace of the
stoichiometry matrix, this algorithm constructs an exhaust-
ive series of tableaux via a Gaussian elimination procedure
along with selection rules akin to those described above.
Detailed numerical examples, illustrating the application of
this algorithm, are presented in various publications by
Schuster [3,4] and notably in the valuable book of Heinrich
& Schuster [7]. A computer implementation of this algo-
rithm is available in the C language [4]. This algorithm
has further been implemented with the
MATLAB
language, in
a package called
FLUXANALYZER
, by S. Klamt and which is
available on the Internet. Moreover,
METATOOLS
by S. Sch-
uster and collaborators is another program in C, which can
also be downloaded from the Internet.
In his publications, Schuster emphasizes the so-called
Ôelementary flux modesÕ. These should not be confused with
the basic colors or vertices of the current polytope. The
elementary flux modes are identical to the basic colors only
2750 J. W. Stucki (Eur. J. Biochem. 271) Ó FEBS 2004
if there are solely irreversible reactions occurring in the

vation conditions. By using the Clarke algorithm, the
stoichiometry matrix needs the deletion of all dependent
species, otherwise the algorithm would attempt to find
solutions of a singular matrix. In simple cases, such as in the
conservation condition, NAD
+
+NADH¼ constant,
these are quite simple to detect. But already in more
complicated models, which at first sight may look trivial,
conservation conditions become a real problem. One
innocently looking example is the glycosome of Trypano-
somes [8]. There, the conservation of the organic phosphates
bound to different carbohydrates already requires some
analysis and could easily escape detection. One of the
biggest advantages of the calculation of the kernel of the
stoichiometry matrix is that all of these complications are
circumnavigated because conservation conditions can be
completely ignored. The basis vectors of the kernel have
automatically eliminated these redundancies and have
reduced the problem to full rank, i.e. containing only
independent species. The explanation for this is that the
NullSpace algorithm uses the row echelon form of the
matrix. Applied to the above-mentioned reaction scheme in
glycosomes, including reversible reactions, exactly three
basic colors are found [9]. When one is interested only in the
proper decomposition of the system into independent and
dependent species, one may use the program
GEPASI
by
P. Mendes that can be downloaded from the Internet.

calculation shows an 81%, 92% and 107% carbon recovery
in the control, with 0.33 m
M
and 0.66 m
M
fluorocitrate,
respectively. It can be estimated that % 50% of the missing
pyruvate units are caused by the fumarase reaction and the
remaining 50% have disappeared in the further metabolism
of the Krebs cycle beyond citrate. In the presence of
fluorocitrate, only fumarate production could explain the
difference, because aconitase is inhibited. Therefore, in
order to arrive at an exact result, the fumarase, as well as the
oxoglutarate dehydrogenase, reactions should be included
in the scheme, and the production of fumarate, ketoglut-
arate and succinate should be measured in the incubation
medium. However, by doing this we would immediately
leave 3D space and the current polytope could no longer be
visualized. This can be exemplified as follows: considering
the full scheme, shown previously [14], under anaerobic
conditions reveals that the Krebs cycle is interrupted at the
succinate dehydrogenase step as no oxidation of FADH
2
is
then possible. In other words, the Krebs cycle is now cut
into two trees: one hanging down from oxaloacetate to
fumarate; and the other from citrate to succinate. A
calculation of the current polytope for this case results
in a 7D geometrical object with 13 vertices. A similar
complication arises when including the exchange carriers

matrix. It seems, however, possible to make a qualitative
analysis based on an adjacency matrix, which contains the
known protein interactions. Using a variant of the methods
described above, useful information about the regulatory
functions of protein interactions can then be obtained. Such
investigations are currently in progress in our laboratory.
Acknowledgements
This work has been supported by grants from the Swiss National
Science Foundation. I am indebted to Dr Clemens Wagner and Dr
Robert Urbanczik for inspiring discussions. The technical expertise of
Mrs Lilly Lehmann, in performing the experiments, is gratefully
acknowledged.
References
1. Clarke, B.L. (1980) Stability of complex reaction networks.
In Advances in Chemical Physics XLIII (Prigogine,I.&Rice,S.A.,
eds), pp. 1–215. John Wiley & Sons, New York.
2. Clarke, B.L. (1981) Complete set of steady states for the general
stoichiometric dynamical system. J. Chem. Phys. 75, 4970–4979.
3. Schuster, S. & Ho
¨
fer, T. (1991) Determining all extreme semi-
positive conservation relations in chemical reaction systems: a test
criterion for conservativity. J. Chem. Soc. Faraday Trans. 87,
2561–2566.
4. Schuster, R. & Schuster, S. (1993) Refined algorithm and com-
puter program for calculating all non-negative fluxes admissible in
steady states of biochemical reaction systems with or without some
flux rates fixed. CABIOS 9, 79–85.
5. Stucki, J.W. (1978) Stability analysis of biochemical systems.
Progr. Biophys. Mol. Biol. 33, 99–187.

dehydrogenase (see Concluding remarks). The transport reactions across the inner mitochondrial membrane were simplified to influx and efflux
without detailed formulations as symports or antiports. The efflux of oxaloacetate was omitted altogether because only neglible concentrations
appear in the incubation medium. In order to allow the electrogenic passage of ATP into the mitochondria through the adeninenucleotide
translocase, the membrane potential was collapsed with valinomycin. The same applies for other nonelectroneutral transports.
No. Reaction Enzyme(s) EC No.
R
1
{}fi Pyruvate Influx
R
2
Pyruvate + CoASH + NAD
+
fi AcetylCoA + NADH + H
+
+CO
2
Pyruvate dehydrogenase 1.2.4.1.
R
3
2 AcetylCoA fi Acetoacetate + 2 CoASH Acetyltransferase + 2.3.1.9.
acetoacetylCoA hydrolase 3.1.2.11.
R
4
Acetoacetate + NADH + H
+
fi 3-OH-Butyrate + NAD
+
3-Hydroxybutyrate dehydrogenase 1.1.1.30.
R
5

11
Citrate fi { } Efflux
R
12
Malate fi { } Efflux
2752 J. W. Stucki (Eur. J. Biochem. 271) Ó FEBS 2004
Appendix 2.
Ó FEBS 2004 Chromokinetics of metabolic pathways (Eur. J. Biochem. 271) 2753
A color finder written in
MATHEMATICA
. This short program was written in
MATHEMATICA
4.
MATHEMATICA
code entries into
the computer are displayed as boldface Courier font. On a Macintosh Cube computer (450 Mhz) the calculation of the basic
colors for anaerobic pyruvate metabolism in mitochondria took 0.02 s. A more involved example, the reaction scheme for the
Belusov–Zhabotinsky reaction [2], took 0.8 s and produced all 34 basic colors. No attempt was made for further speed
optimization, such as matrix preprocessing or jumping out of prematurely finished loops. The major motivation to publish
this program here for the first time is to give the reader an interactive tool with which he or she can explore their own
problems. To do this, the user first has to enter his own stoichiometry matrix into
MATHEMATICA
’s front end, as in step (1). It
is mandatory to include reverse (backward) reactions as separate additional columns to the stoichiometry matrix. Then, the
MATHEMATICA
code from steps (2) to (6) has to be copied exactly. Step (5) contains timing and printing facilities, which may
be omitted. For large systems, containing many hundreds of basic colors, only their number may be of interest. This is
accessible with the last statement of the program. It is good practice to start a new problem with a new
MATHEMATICA
kernel.