What makes biochemical networks tick?
A graphical tool for the identification of oscillophores
Boris N. Goldstein
1
, Gennady Ermakov
1
, Josep J. Centelles
3
, Hans V. Westerhoff
2
and Marta Cascante
3
1
Institute of Theoretical and Experimental Biophysics, Russian Academy of Sciences, Pushchino, Moscow Region, Russia;
2
BioCentrum Amsterdam, Departments of Molecular Cell Physiology (IMC, VUA) and Mathematical Biochemistry (SILS, UvA),
Amsterdam, the Netherlands;
3
Department of Biochemistry and Molecular Biology, Faculty of Chemistry and CeRQT at Barcelona
Scientific Parc, University of Barcelona, Spain
In view of the increasing number of reported concentration
oscillations in living cells, methods are needed that can
identify the causes of these oscillations. These causes always
derive from the influences that concentrations have on
reaction rates. The influences reach over many molecular
reaction steps and are d efined by the d etailed molecular
topology of the network. So-called Ôautoinfluence pathsÕ,
which quantify the influence of one molecular species upon
itself through a particular path through the network, can
have positive or negative values. The former bring a ten-
dency towards instability. In this molecular context a new

mathematical models [6,7]. The need for such mathematical
models is appreciated even more w hen studying biochemical
oscillations and their synchronization [7–13].
The behavior of potential biochemical oscillators may
depend on the kinetic properties of t heir surroundings,
interacting with t he oscillator through c ommon metabolites
(e.g [8,14]). Other systems, such as the cell cycle of tumor
cells may be more autonomous [9]. Most intracellular
oscillations involve more than five components that interact
in a nonlinear man ner [8]. This makes them unsuitable for
intuitive analysis, a phenomenon encountered more fre-
quently in Systems Biology [8]. New theoretical approaches
are needed that streamline the study of such cases of
Systems Biology, dissecting the system into various inter-
acting kinetic regimes, whilst relating to molecular mecha-
nisms.
Various types of approach can be helpful here. Graph-
theoretic approaches can help dissect the dynamics of
enzyme reactions [15,16] and this is what made others and
ourselves [20,25,26] examine whether these approaches can
also do this for networks. Earlier w e have applied graph
theory in order t o simplify the King–Altman–Hill [15,16]
analysis of steady-state enzyme reactions [17,18]. This
approach was later extended to presteady-sta te enzyme
kinetics [19], to s tability a nalysis o f enzyme systems [20],
and to the analysis of concentration oscillations in enzyme
cycles [21].
In this paper, the graph-theoretical stability analysis
developed b y C larke [ 22] as modified b y I vanova [21,23,24]
is the starting point for a more comprehensive approach

which oscillations occur. Presence or absence of steady
states on the border of the phase space [21,23] then suffices
to predict the occurrence of limit-cycle oscillations. We
illustrate our method by applying it to two biochemical
systems, which include oscillophores of two different classes.
Results
Paths: graphical representation of kinetic influences
We represent kinetic schemes b y dual g raphs, combining
reaction-centered and substance-centered graphs [21].
Accordingly, our kinetic schemes for biochemical networ ks
have two kinds of vertices, i.e. one kind for species (here
shown by open circles) a nd one kind for reactions (shown by
closed circles). The circles are connected by arrows. For
example, the reaction x
i
+ x
j
fi x
m
is represented by the
following reaction-centered graph:
where x
i
, x
j
and x
m
are c hemical s pecies (substances ) and v
r
is the rate of the r

i
Á x
a
jr
j
ð1Þ
where k
r
is the k inetic constant. T his implies t hat we do not
dissect biochemical networks i nto t he net enzyme-catalyzed
reactions, but into the unidirectional elementary r eaction
steps underlying the enzyme kinetics. The terms x
i
, x
j
and x
m
include the concentrations of both m etabolites and enzyme-
forms. Rates v
r
are always positive. As a consequence of the
dissection down to t he molecular p rocesses, the stoichio-
metric coefficients equal one or zero, i.e.
a
ir
; b
ir
¼ 1or0 ð2Þ
with 1 for participating and 0 for nonparticipating species.
Reactions involving more than one molecule of a single

À!
a
ir
ðrÞ
Graph 2:
Similarly to the procedure developed b y Clarke [22], we
linearize the system of Eqn (3) in the vicinity of the steady
state. We do this to investigate the stability of this state. In
this way we obtain the influence a small change in the
concentration of substance j, i.e. Dx
j
,hasonthetime
displacement of t he conc entration of species i from its
steady-state value:
dDx
i
dt
¼
X
r;j
ðÀa
ir
þ b
ir
Þ
@v
r
@x
j
Dx

v
r
x
j
Dx
j
¼
X
j
b
ij
Dx
j
ð6Þ
Coefficients b
ij
are the elements of the Jacobian (matrix) B
representing the direct influences of x
j
on x
i
:
b
ij
¼
X
r
ðÀa
ir
þ b

j
, t he sign o f its contribution depending on the direction
of the reaction, as specified by Eqn (7). The terms that
multiply an a and a b therewith represent the positive
influence that a substrate of a reaction has on the product o f
the r eaction. The terms of b
ij
that multiply two a’s, represent
the negative influences of two substances on each other
when both are consumed in that reaction. Indeed, each
element of the Jacobian corr esponds to one or a number of
such direct influences of one metabolite on another, direct in
the sense that the influence is through s ingle reaction steps.
A number of s uch reaction steps may operate in parallel
(but not in series) for each element of the Jacobian. In
addition, one reaction step may convey more than one
influence.
Graph 1.
3878 B. N. Goldstein et al.(Eur. J. Biochem. 271) Ó FEBS 2004
To predict t he dynamics of the system it is i ndeed helpful
to classify these influences into positive and negative ones.
Again, depending on whether the reaction stoichiometry
is an a or a b, two types of influence are seen in Eqn (7).
They are shown graphically in Eqns (8) and (9) together
with their corresponding contributions to b
ij
:
x
j
À!

a
jr

v
r
À
a
ir
x
i
; b
ij
¼ a
jr
ÁÀa
ir
Á
v
r
x
j
¼À
v
r
x
j
ð9aÞ
Although they are similar to Graph (1), Eqns (8) and (9a)
have meanings that differ from t he meaning of Graph (1):
They do not represent chemical conversions but rather the

ij
.
They correspond to the influence of a substrate on a nother
substrate o f the same reaction r. If i ¼ j, Eqn (9a) defines a
so-called negative half-step instead of a negative step (we
omit ÔinfluenceÕ for brevity):
x
i
À!
a
ir

v
r
À
a
ir
x
i
; b
ii
¼ a
ir
ÁÀa
ir
Á
v
r
x
i

r
x
i
ð9cÞ
hence its name Ôhalf-stepÕ. This is the (negative) influence a
substrate has on its own removal. It is obtained for all
substrates of any elementary reaction.
The main point of the present section is that any Jacobian
matrix element equals the sum of a number of direct parallel
influence steps (one-step influence paths) in the kinetic
scheme, i.e. the sum of paths through reaction-centered
graphs of the type o f Graph 1 (these paths may contain
parallel and antiparallel arrows). The sign of that element
therefore depends on the both the sign and the magnitudes
of these influen ce paths (see below). If all its i nfluence paths
are positive, the Jacobian matrix element will be positive
and for the J acobian matrix eleme nt to b e negative at l east
one influence path must be negative. These are properties
that we shall use below.
How graphical structures relate to instability
For the linear system given in Eqn (6) the so-called
characteristic polynomial p(k)is:
pðkÞ¼detðB À kIÞ¼0 ð10Þ
Here B is again the Jacobian with elements b
ij
and I is the
unit matrix. The polynomial Eqn (10) can be expanded as
follows
pðkÞ¼k
m

2
Á

X
i;j
b
ii
b
jj
À
X
i;j
b
ij
b
ji

;
a
3
¼ðÀ1Þ
3
Á

X
i;j;k
b
ii
b
jj

p
equals the sum of all
possible Ôp
th
order simplest combinations of minus auto-
influence pathsÕ. An autoinfluence path (or, shorter, a cycle)
is defined as a cyclic path of any length through the
diagram, such that any reactio n and any species occurs only
once on that path. Autoinfluence paths of lengths 1, 2 , 3, etc.
correspond to the terms b
ii
, b
ij
b
ji
, b
ij
b
jk
b
ki
, etc., respectively,
inEqn(12).Theycontain1,2,3,etc.speciesand1,2,3,etc.
influence steps, respectively. Autoinflu ence paths of length 1
are h alf-steps, graphically represented as i n Eqn (9c). An
autoinfluence path runs from some species k back to species
k an d travels through positive i nfluence steps [reaction
nodes with equally directed arrows, as in E qn (8)] or
negative influnce s teps [reaction nodes with oppositely
directed arrows, as in E qn (9a)]. Consequently, a Ôminus

network, is the basis of the graphical analyses of the
characteristic equation and of t he method we develop here.
Inspection of Eqn 12 shows t hat in a ll coefficients a
p
the
term consisting of negative half-steps (Eqn 9c) only, which
corresponds to products of Jacobian elements b
ii
only, is
always positive: of the term of o rder p the sign is ( )1)
p
multiplied by ()1)
p
. I ndeed, a ll these terms always constitute
negative combined autoinfluence paths.
This graphical procedure allows us to determine all
coefficients of the characteristic polynomial for systems o f
simple reactions. The graphical determination of character-
istic polynomial co efficients for complex stoichiometries has
been elaborated by Ivanova [23].
The concentrations are restricted by balance constraints
(conserved sum concentrations, such as NADH + NAD)
and by the requirement that they be positive. These
restrictions define upper and lower limits for the values
the concentrations can assume (i.e. borders of the so-called
phase space). Any negative a
i
coefficient implies that the
system is unstable [22]. Such instability could lead t o infinite
growth (explosion) of some concentrations, unless the

structures, i.e. ÔsubgraphsÕ (see below), in t he kinetic s cheme
that contribute t erms to the coefficients of the characteristic
polynomial of a predictable sign and that hence help
determine th e stability properties of t he system. The aim of
this paper is to identify ÔnegativeÕ subgraphs, because they
can induce instability; their positive combined autoinfluence
bestows them with oscillophoretic potential.
The instability condition that a
p
be negative for s ome
p < m translates to the c ondition that the positive
combined autoinfluence paths of order p should outweigh
the negative combined autoinfluence paths of that same
order. From this, an Ôinstability ruleÕ follows. This is stated
as ÔInstability is promoted (counteracted) by positive
autoinfluence paths.Õ This c onnotes with instability be ing
generated by positive feedback loops.
Subgraphs favoring instability
As mentioned above we d eal here with the formulation t hat
decomposes biochemical networks into truly elementary
reactions. A network then c onsists of a great many such
reactions (each represented as a black node in our reaction
equations), each of which connects a number of species
(represented as white nodes). The entire network may
become unstable when part of it w ould by itself be unstable.
Consequently it can be useful to identify parts of the larger
network that are unstable.
The g raphical representation o f a subnetwork with e qual
numbers of species a nd reactions is here called a subgraph.
It is useful to consider subgraphs because all combined

), there is both a negative half step
(because, as usual, x
i
stimulates its own removal; Eqn 9c)
and a positive l oop [becau se x
i
now also stimulates its own
production; compare Eqn (8) with i ¼ j ]. Adding these two,
Eqn (7) shows t hat they cancel each other:
Àa
1
3 b
ii;r;1
¼ðÀa
ir
a
ir
þ b
ir
a
ir
ÞÁ
v
r
x
i
¼ 0; ð13Þ
where the symbol ’ means ÔcontainsÕ. The value of zero is
obtained because a ll stoichiometric coefficients equal one (i.e.
the reaction x

the substrate and not as the product. Therefore, a
1
is
constructed only from the corresponding negative half-steps
with the values [compare Eqns (7) and (9c)]:
Àa
1
3 b
ii;2
¼Àa
ir
a
ir
Á
v
r
x
i
¼À
v
r
x
i
ð14Þ
This corresponds to a single cycle (from x
i
back onto itself)
with a single n egative i nfluence step, i.e. it is negative in
terms of autoinfluence and promotes stability.
The sum total f or one-reaction subgraphs is thereby

ii
b
jj
¼ÀðÀa
ir
a
ir
ÞÁ ðÀa
js
a
js
Þþa
js
b
js
ÀÁ
Á
v
r
v
s
x
i
x
j
¼Àða
ir
a
ir
a

ir
b
jr
ÞÁ
v
r
v
s
x
i
x
j
¼
v
r
v
s
x
i
x
j
ð17Þ
In the first facto r of Eqn (16), which corresponds to the
direct self-influence t erm b
ii
, one recognizes the influence
that x
i
has on itself through its own d egradation v
r

back onto itself which is
negative. O ne of th em multiplies w ith the negative half-step
of x
j
back on to itself and constitutes n egative a utoinfluence
(even number of cycles and even number o f negative
influences). The other multiplies w ith the positive loop of x
j
back onto itself through v
s
: two cycles with one negative
influence constituting a positive (destabilizing) autoinflu-
ence. These two autoinfluence paths of order two cancel
each other.
In Eqn (17) one re cognizes the influence that x
j
has o n x
i
(+a
js
b
is
) b ecause t he former is the substrate of the reaction
v
s
that produces the latter, as well as the a nalogous influence
x
i
has on x
j

, cannot be negative, and hence
cannot cause instability. Then only the negative autoinflu-
ence path of Eqn (16) remains, which then cancels the
positive autoinfluence of Eqn (17). Revolving around the
cycle then does not lead to an increase in the number of
molecules. Elementary reactions producing more types of
product than types of substrate are essential for the
occurrence of instability, due to the restriction s on
stoichiometries that d erive f rom ou r descent t o t he
molecular level.
Another branched cycle, Graph (4), having two bran-
ches, contributes to –a
2
the same positive, destabilizing term:
Graph 3.
Graph 4.
Ó FEBS 2004 Graphical tool to identify oscillophores (Eur. J. Biochem. 271) 3881
Negative Graphs (3) and (4) involve positive paths with
branching that can be interpreted as positive feedback
interactions (autocatalysis). They can also be interpreted as
product activation in some enzyme reactions, because a
reaction product here stimulates the same re action. Another
example of such positive paths of influence occurs in the
case of the antiport of two ligands by a protein molecule
through the membrane [24].
Three-reaction subgraphs. Inthesamewayweidentify
the negative ( instability generating) subgraphs with three
species and three reactions, by pointing out that they
have positive p aths of influence. We divide these g raphs
into two classes, i.e. those with positive influence steps

b
23
À b
11
b
23
b
32
À b
22
b
13
b
31
À b
33
b
12
b
21
Þ
The first term here corresponds to all h alf-steps multiplied,
the second and the third terms correspond to the circular
paths running through all three reactions and all three
species. The other three terms correspond each to a single
half-step, multiplied by a circular p ath running through two
reactions and two species.
Using the graphical rules mentioned above, we now
consider the left hand subgraph in Graph (6) as an example.
This subgraph contains the following three simplest com-

12
b
12
b
33
) would run through the two
reactions on the right hand side (negative), a nd one
negative half-step on the left (positive). However, the
latter would touch the cycle, hence this one does not
count. The third term, i.e. b
11
b
23
b
32
is also ÔemptyÕ for this
diagram.
3. The third type of combined autoinfluence is negative: a
combined influence of three separate anti half-steps: )1.
The sum of these t wo positive and one negative combined
autoinfluences contributes a negative term into the coeffi-
cient a
3
, and therewith promotes instability. The subgraph
on the left of Graph (6) is therewith negative.
In fact all of t he subgraphs in G raphs (5) a nd (6)
contribute the positive term
v
1
v

obtain sustained oscillations, the efflux reaction should be
reversible, leading, for example, to an inhibitory enzyme
complex (see below for an example). In the latter case the
reversible steady-state efflux equals zero and the negative
subgraph is upheld.
n-Reaction subgraphs. We can now formulate the proper-
ties of any negative subgraph t hat contains an arbitrary
equal number of species and reactions. Such a negative
graph should be constructed at least of two even cycles,
formed by a branched reaction. Moreover, species of the
negative subgraphs should not be connected with other
parts of the full scheme by outgoing irreversible reactions.
The outgoing i rreversible reactions cause t he correspond-
ing opposite stationary fluxes to be equal, canceling the
negative s ubgraph with a positive graph of th e same
absolute value (see b elow for a n example). Therefore,
only damped oscillations can be obtained in such a case.
Additional reversible reactions, leading through the
species of the n egative graph to dead-end species , do
not eliminate the negative graph.
Graph 5.
Graph 6.
3882 B. N. Goldstein et al.(Eur. J. Biochem. 271) Ó FEBS 2004
Examples of biochemical oscillations
Two interacting enzymes. Here we discuss one of the
simplest biochemical osc illators. Its kinetic scheme (Fig. 1)
contains a negative graph [Graph (4)] of second order (two
species and two reactions). This sug gests that the system of
Fig. 1 may oscillate. A more detailed analysis should then
be undertaken to determine whether it will actu ally oscillate.

is conserved:
½E
2
-Pþ½E
2
¼E
2
or x
4
þ x
2
¼ 1 ð18Þ
The reaction participants and their normalized concentra-
tions x
i
(i.e. their concentrations divided by the total
concentrations of E
2
) are shown i n Fig. 1. The following
kinetic equations correspond to Fig. 1:
dx
1
dt
¼ k
2
x
2
x
3
À k

1
, x
2
, x
3
.
We shall now analyze the characteristic polynomial o f
the t hird order f or Eqns (19). We know that the
coefficient a
1
is always positive, and it is readily seen that
the coefficient a
3
is also positive h ere. The coefficient a
2
contains a negative t erm, which corresponds to the
negative graph highlighted in Fig. 1 by the heavy lines.
This negative term equals –v
1
v
2
/x
1
x
2
. At steady state v
2
¼
v
4

Fig. 1, one recognizes that irreversible effluxes of species
of the negative graph must be present in order for the
steady state condition to be satisfied. For the same
reason, all the biochemical schemes involving a negative
graph of second order can only induce damped oscilla-
tions. Damped oscillations calculated for F ig. 1 are
shown i n F ig. 2.
Substrate inhibited bifunctional enzyme. Many kinetic
graphs that generate oscillations with only positive auto-
influence p aths are known from the literature. Some of them
have been classified [25]. A lthough the graphs with positive
paths implemented here e.g. Gra phs (3–5) are simpler t han
the kinetic graphs in [25], t hey can represent biochemical
reality. For example, t he second of the subgraphs in
Graph (5) has been used to analyze the network topological
basis for oscillatory antiport of t wo different ions across the
cell membrane [24].
Less studied are the graphs that include negative paths,
such as those in G raph (6). Two negative p aths in the c ycle
of subgraphs here reflect the competition of two r eactions
for a single species. For example, the competition of protein
X and the enzyme E
2
for the acetyl group in pyruvate
dehydrogenase c omplex has been shown to be important for
the prediction of o scillatory behavior [28].
The phenomenon of substrate inhibition is often associ-
ated with the potential for oscillations [33]. An earlier graph-
Fig. 2. Calculated time dependence of the normalized E
1

Experimental Bi ophysics, RAS, M oscow Region, Russia).
Fig. 1. Reaction scheme of two enzymes d emodify ing e ach other. Filled
circles represent reactions, open circles represent substances. T he rate
of the reaction t hat combines E
2
with P to yield E
2
-P is given as v
3
.IfP
represents a p hosphate group, reaction nu mbe r 3 could be a p rotein
kinase, and v
1
should t hen represent the dephosphorylation of E
2
-P, as
catalyzed by E
1
. This reaction re leases P a nd E
2
. In this reaction E
1
is
used but im mediately r eleased a s it i s a catalyst. The rate of reaction 2 i s
v
2
, which is catalyzed by E
2
(which then functions as a protein p hos-
phatase) and dephosphorylates E

fi S
2
fi S
1
,as
catalyzed by the bifunctional enzyme (E
1
/E
2
;E
1
and E
2
are
two states of a single protein). H ere t he arrows between
symbols correspond to the preferential reaction orientation.
In reaction 1, E
1
catalyzes the f orward reaction S
1
fi S
2
and E
2
catalyzes reaction 2 , which runs in the opposite
direction, S
2
fi S
1
. Reaction 2 may be coupled to a s ource

3
. G raph (7) does not contain irreversible effluxes from
the s pecies (E
2
,E
3
and S
2
) of the negative graph that i s
shown b y h eavy ar rows, a nd contai ns only the influx to S
2
,
i.e. 1 fi S
2
(from reaction 1 to S
2
). Therefore it retains its
oscillophoretic potential. Inhibitory reversible reactions,
added to the negative graph, do not interfere with that
potential.
The subgraph highlighted by the heavy arrows in the
full Graph (7), is one of the negative g raphs identified
in this paper [i.e. the second left of the subgraphs in
Graph (6)]. This negative graph is the b ranched cycle
with one positive loop, represen ting the E
3
catalyzed
reaction tha t ma kes E
2
out of E

not shown.
Oscillations can be expected if we add reversible inhibi-
tion of E
3
by substrate S
1
to Graph ( 7). The reversible
inhibitions of E
3
by both S
2
and S
1
do not eliminate the
negativity of the negative subgraph, because these reversible
steps d o not contribute additional terms to the terms of the
negative graph. Their contributed effluxes are eq ual to
influxes. However, t he number of s pecies of reactions
becomes larger with this new inhibition. Accordingly, a
positive graph with four species and four reactions, as well
as a negative graph with three species and three reactions,
are obtained in Graph (7). This is a sufficient condition for
oscillations to arise.
We shall now show how a necessary condition for
oscillations to occur follows from the absence of steady
states on the border of the phase space. The full s ystem
contains seven species variables:
x
1
¼½E


These species are interdependent through the following
three balance constraints:
x
2
þ x
4
þ x
6
þ x
7
¼ S ¼ constant
x
1
þ x
3
¼ E ¼ constant ð20Þ
x
5
þ x
6
þ x
7
¼ E
0
¼ constant
These constraints reflect the conserved total concentra-
tions of the substrates (we shall use S ¼ 3.3 relative units) of
the b ifunctional e nzyme (E ¼ 0.2 relative units), and of the
modifying enzyme ( E¢ ¼ 0.31 relative units). The following

6
and v
7
relate to the reversible reaction E
3
+S
1
b « E
3
S
1
.The
equalities [Eqn (21)] together with the constraints [Eqn (20)]
allow u s t o obtain all seven concentration values for one of
the steady states in the phase space [E
1
] ¼ 0, [S
2
] ¼ 0,
[E
3
S
2
] ¼ 0, [E
2
] ¼ E, and [S
1
], [E
3
], [E

¼Àk
1
x
1
x
2
þ k
2
x
3
x
4
À k
6
x
2
x
5
þ k
7
x
7
dx
3
dt
¼ k
3
x
1
x

x
6
ð22Þ
In addition to referring to the a bsence of steady states on the
border of t he phase space, the procedure b y Clarke [22]
enables us to identify qualitatively phase trajectories that
lead to a stable limit c ycle. The c haracteristic polynomial o f
the system in Eqn (22) reads:
k
4
þ a
1
k
3
þ a
2
k
2
þ a
3
k þ a
4
¼ 0 ð23Þ
If in this polynomial a
4
> 0 for all concentration values
and a
3
< 0 in the unstable steady state, oscillations can be
obtained. Analysis of negative and positive subgraphs and

is unraveling more and more of the molecular specifics
that underlie cell fun ction.
Our method to classify potential biochemical oscillators is
based on t he graphical analysis of t he kinetic schemes. Our
approach is similar in s ome aspects to the procedure
described previously [25,26]. However, the representation o f
the k inetic schemes i n terms of dual grap hs [ 21] is different,
and has enabled us t o simplify the identification a nd the
classification of oscillophoretic networks.
Because above we were most concerned with demon-
strating the basis of our method, we here summarize how
the approach may be implemented in the context of a
known reaction network. First the network kinetics
should b e drawn o ut i n a detailed molecular s cheme
making all molecular interactions, such as the binding of
a ligand to a n enzyme, exp licit. Then one should try to
recognize subgraphs of known s ign in that s cheme. Here
one may resort to the subgr aphs identified i n this p aper,
or to subgraphs that may appear i n future w ork
analyzing networks more extensively. Alternatively, one
may u se the method of making an inventory o f t he
autoinfluences within each subgraph and determine
whether there are more positive ones than negative ones,
in which case the subgraph is negative (unstable). Having
identified the (negative) subgraphs with oscillop horetic
potential, one may then analyze their effect quantitatively
and compare the results to those obtained for through
analysis of all other subgraphs of the same order in the
same network, as was illustrated for the two examples in
this paper. The network outside the former s ubgraph

hierarchy, different methods for the analysis of th e d ynam-
ics are needed.
We demonstrated how reaction networks that are
formulated down to the detail of simple unimolecular and
bimolecular reactions can be organized into topologies. The
latter can then be examined for their potential to induce
oscillations. Oscillophoretic topologies involve branched
directed cycles, c onstructed of a n even number of negative
paths a nd any number of positive paths. Our approach has
the advantage th at it considers positive and negative
interactions in a unified manner.
The implication o f the identification o f an oscillophoretic
subgraph is that if such a subgraph is found in a large
network, then that network may be unstable and give rise to
oscillations; the presenc e of an oscillophoretic subgraph is a
necessary condition for the network t o engage in t he
oscillations. However, i t is not a sufficient condition.
Ó FEBS 2004 Graphical tool to identify oscillophores (Eur. J. Biochem. 271) 3885
Whether the overall network actually engages in an
oscillation when an oscillophoretic subgraph is present
depends on the precise parameter values. To estimate the
parameter domain where oscillatory phenomena can be
observed, the numerical value of the negative graph should
be compared with the values of other graphs of the same
order in the system. In practice this means that to produce
oscillations, reactions involved in the negative graph should
be rapid enough a s compared with their surrou nding
reactions. We here p erformed such an analysis for two
examples, one with positive and one with negative inter-
actions.

component properties that contribute to the tendency of a
system to engage in more complex behavior such as limit-
cycle oscillations. Actual and subtle interactions of the
components then determine whe ther or not the oscillations
actually occur.
Acknowledgements
This work was supported by a grant from Ministry of Science and
Technology of the Spanish Government (SAF 2002–02785), INTAS
grant (97–1504), and the Netherlands’ Organization for Scientific
Research. We t hank T. Sukhomlin for disc ussions.
References
1. Betz, A. & Sel’kov, E. (1969) Control of phosphofructokinase
(PFK) activity in conditions simulating those of glycolysing yeast
extract. FEBS Lett. 3, 5–9.
2.Richard,P.,Bakker,B.M.,Teusink,B.,VanDam,K.&
Westerhoff, H.V. (1996) Acetaldehyde mediates the synchroni-
zation of sustained glycolytic oscillations i n p opulations of yeast
cells. Eur. J. Biochem. 235, 238–241.
3. Berridge, M.J., Bootman, M.D. & Lipp, P. (1998) Calcium – a life
and death signal. Nature 395, 645–648.
4. Schuster, S., Marhl, M. & Ho
¨
fer, T. (2002) Modeling of simple
and complex calcium oscillations. Eur. J. Biochem. 269, 1333–
1355.
5. Nurse, P. (2001) The cell cycle and development. Introduction.
Novartis Fou nd Symp 237, 1–2.
6. Heinrich, R., Rapoport, S.M. & Rapoport, T.A. (1977) Metabolic
regulation and mathematical models. Progr. Biophys. Mol. Biol.
32, 1–82.

17. Volkenstein, M.V. & Goldstein, B.N. (1966) A new method for
solving the problems of stationary kinetics of enzymological
reactions. Biochim. Biophys. Acta 115, 471–477.
18. Volkenstein, M.V. & Goldstein, B.N. (1966) Allosteric enzyme
models and their an alysis by the theory of graphs. Bioc him. Bio-
phys. A cta. 115, 478 –485.
19. Goldstein, B.N. (1983) Analysis of cyclic enzyme reaction
schemes by graph – the oretic method. J. Theor. Biol. 103 , 247–
264.
20. Goldstein, B.N. & Shevelev, E. (1985) Stability of multienzyme
systems with feedback regulation: a graph theoretical approach.
J.Theor. Biol. 112, 493–503.
21. Goldstein, B.N. & Ivanova, A. (1987) Hormonal regulation
of 6-phosphofructo-2-kinase/fructose-2,6-biphosphatase: kinetic
model. FEBS Lett. 217, 212– 215.
22. Clarke, B.L. (1980) Stability of complex reaction networks. Adv.
Chem. Phys. 43 , 1–115.
23. Ivanova, A.N. (1979) Conditions for the unique steady state of
kinetic systems as connected with struc ture of the reaction
schemes. Kinetika i K ataliz (in R ussian) 20, 1 019–1028.
24. Goldstein, B.N., Holmuhamedov, E.L., Ivanova, A.N. & Fur-
man, G.A. (1987) Ionophore-induced osc illations in erythrocytes.
Mol. Biol. 21, 1 10–117.
25. Eiswirth, M., Freund, A. & Ross, J. (1991) Operational procedure
toward the classification of chemical oscillators. J. Phys. Chem. 95,
1294–1299.
3886 B. N. Goldstein et al.(Eur. J. Biochem. 271) Ó FEBS 2004
26. Eiswirth, M., Freund, A. & Ross, J. ( 1991) Operational procedure
toward the classification of chemical oscillators. Adv. Chem. Phys.
190, 127–150.

mitochondrial cycles of energization. FEBS Le tt. 519, 41–44.
38.Wolf,J.,Sohn,H.Y.,Heinrich,R.&Kuriyama,H.(2001)
Mathematical analysis of a mechanism for autonomous metabolic
oscillations in continuous of Saccharomyces cerevisiae. FEBS Lett.
499, 230–234.
39. Westermark, P.O. & Lansner, A. (2003) A model of phospho-
fructokinase and glycolytic oscillation s in the pancreatic cell.
Biophys. J . 85, 126–139.
40. Lahav, G., Rosenfeld, N., Sigal, A., Geva-Zatorsky, N., Levine,
A.J., Elowitz, M.B. & Alon, U. (2004) Dynamics of the p53-
Mdm2 feedback loop in individual cells. Nat. Genet. 36, 113–114.
41. Brady, N.R., Elmore, S.O., Van Beek, J.H.G.M., Krab, K.,
Courtoy, P.J., Hue, L. & Westerhoff, H.V. (2004) Coordinated
behavior of mitochondria in both space and time: a reactive
oxygen species-activated wave of mitochondrial depolarization.
Biophys. J . in press.
42. Westerhoff, H.V. & Van Dam, K. (1987) Thermodynamics
and Control of Biological Free Energy Transduction. Elsevier,
Amsterdam.
Ó FEBS 2004 Graphical tool to identify oscillophores (Eur. J. Biochem. 271) 3887