Reduction of a biochemical model with preservation
of its basic dynamic properties
Sune Danø
1
, Mads F. Madsen
1
, Henning Schmidt
2
and Gunnar Cedersund
2
1 Department of Medical Biochemistry and Genetics, University of Copenhagen, Denmark
2 Fraunhofer Chalmers Research Centre for Industrial Mathematics, Gothenburg, Sweden
Systems biology aims to understand the behaviour of
biological systems, and in particular how the system’s
behaviour emerges from the interactions among its
components. Consequently, there is an increased focus
on biochemical dynamics and its relation to the under-
lying biochemical reaction network. This connection
is most often described by means of mathematical
models with varying degrees of detail. A primary
advantage of a detailed, biochemically formulated
model is that a one-to-one comparison can be made
between model and biochemistry. A major disadvan-
tage of such a full-scale model stems from its large
number of parameters. A large number of parameters,
compared with the information available from experi-
ments, makes the model unidentifiable. This means that
there are an infinite number of parameter combinations
Keywords
core model; glycolysis; Hopf bifurcation;
model optimization; model reduction

applications of the two methods result in models of eight (lumping), six
(elimination) and three (lumping followed by elimination) dimensions. All
models have similar dynamic properties and pin-point the same interactions
as being crucial for generation of the oscillations. The advantage of the
novel method is that it is algorithmic, and does not require input in the
form of biochemical knowledge. The lumping approach, however, is better
at preserving biochemical properties, as we show through extensive analy-
ses of the models.
Abbreviations
ACA, acetaldehyde; ADH, alcohol dehydrogenase; BPG, 1,3-bisphosphoglycerate; DHAP, dihydroxyacetone phosphate; ENVA, elimination of
nonessential variables; F6P, fructose 6-phosphate; FBP, fructose 1,6-bisphosphate; G6P, glucose 6-phosphate; GAP, glyceraldehyde-3-
phosphate; GAPDH, glyceraldehyde-3-phosphate dehydrogenase; Glc, glucose; HK, hexokinase; LASCO, lumping and subsequent
constrained optimization; ODE, ordinary differential equation; PFK, phosphofructokinase; PGK, phosphoglycerate kinase; PK, pyruvate kinase;
Pyr, pyruvate; trioseP, triosephosphates.
4862 FEBS Journal 273 (2006) 4862–4877 ª 2006 The Authors Journal compilation ª 2006 FEBS
which gives rise to virtually identical agreements with
the data [1]. Therefore, it is impossible to decide which
parts of the model’s predictions are well-supported,
and which are more or less arbitrary. In this way,
model complexity renders many of the advantages of
the one-to-one correspondence useless [2]. Large num-
bers of variables and reactions are also associated with
problems regarding, for example, numerics and model
analysis [3].
A number of model-reduction techniques for com-
plex chemical kinetics have been developed in order to
deal with such problems. As reviewed by Okino &
Mavrovouniotis [3], most model-reduction techniques
fall into three classes: lumping methods, techniques
based on sensitivity analysis, and timescale-based tech-

methods is explained by the fact that no single method
is superior in all cases. The applicability of a method
depends on both the objective with the reduction, and
on the type of complexity in the original model. There-
fore, test case studies comparing the consequences
of different model-reduction methods are of interest.
We chose the cyanide-induced glycolytic oscillations
observed in starved yeast cells, because this is a partic-
ularly well-studied biochemical model system. The
experimental system has been thoroughly characterized
in terms of both biochemistry and dynamics [16–23],
and this has led to a number of mathematical models
of this system [4,24–30].
Our 20D model [30] is a full-scale model that des-
cribes the system in detailed biochemical terms. It is in
quantitative agreement with almost all experimental
observations, but it suffers from the above-mentioned
general problems of detailed models. In contrast, we
have shown that the persistent oscillations can be des-
cribed as a 2D phenomenon [20,31]. Even though the
structure of the biochemical reaction network is not at
all present in this 2D model (Eqn 1), it is possible to
obtain biochemical interpretations of the two modes
involved in the oscillatory dynamics [31]. This raises
the general question to what extent can a full-scale
model be reduced to a smaller biochemically meaning-
ful model, with preserved basic dynamic properties?
This study addresses this question for the specific test
case of the 20D model developed in Hynne et al. [30].
The most basic dynamic property of the 20D model

being the value of p at the
bifurcation point.) The real parameter x
0
is the fre-
quency of oscillations at the bifurcation point, and the
imaginary part of the complex parameter r determines
the l-dependency of the frequency at the stationary
state. Re(r) determines the l-dependency of the linear
stability, and hence the direction of the bifurcation.
The complex nonlinearity parameter g determines the
properties of the limit cycle, which is born in the Hopf
bifurcation. In particular, the Hopf bifurcation is
supercritical if Re(g) < 0 and subcritical if Re(g)>0.
S. Danø et al. A case study in model reduction
FEBS Journal 273 (2006) 4862–4877 ª 2006 The Authors Journal compilation ª 2006 FEBS 4863
The Stuart–Landau equation is the 2D description dis-
cussed above, and the parameters x
0
, r and g can be
calculated from a full-scale model at a Hopf bifurca-
tion [34]. As such, it provides a firm connection
between the full-scale model and its basic dynamic
properties.
In this study we evaluate two model-reduction meth-
ods. The first is based on lumping and subsequent con-
strained optimization (LASCO); it is the optimization
step which involves new stages.
The second is the elimination of nonessential varia-
bles (ENVA). This is a new model-reduction method
with a philosophy similar to that of sensitivity analy-

bifurcation, when the mixed flow glucose concentration
[Glc
x
]
0
is used as bifurcation parameter [20] (Glc, glu-
cose). If these two properties are preserved, the model
is said to be in qualitative agreement with the 20D
model (as well as the experimental observations).
Good quantitative agreement is also said to be found
when the Stuart–Landau parameters, i.e. parameters
x
0
, r and g of Eqn (1), are in reasonable agreement
with those of the 20D model.
Figure 1 provides an overview of the models devel-
oped here. Two model-reduction strategies are applied.
The first, LASCO, is the elimination of nonessential
reactions, commonly known as lumping. The other,
ENVA, is the elimination of nonessential variables.
Starting with the 20D model [30] (Fig. 2), we use
LASCO to arrive at the 20L8D model, and ENVA to
arrive at the 20E6D model. The 20L8D model was fur-
ther reduced by ENVA, resulting in the 20LE3D
model.
We now describe the three model reductions in
detail, and present the resulting models.
Construction of the 20L8D model by LASCO
A traditional approach to model reduction is lumping.
In essence, a simpler model structure is obtained by

so that the dynamic properties of the 20L8D model are as similar
as possible to those of the original 20D model. Details of this pro-
cedure are given in the text. Application of ENVA reduced the
20L8D model to the 20LE3D model.
A case study in model reduction S. Danø et al.
4864 FEBS Journal 273 (2006) 4862–4877 ª 2006 The Authors Journal compilation ª 2006 FEBS
7D model [4]. We have adopted this model structure
here, as an example of a lumped model structure.
(Brusch et al. [37] have previously developed a modi-
fied version of the model by Wolf and Heinrich with
the same purpose as we have here. We, however, find
this model unsuitable due to problems concerning the
formulation of the model.) In order to be able to map
the 20D model onto the reduced model in a straight-
forward manner, we made the following modifications
to the model structure: extracellular glucose was intro-
duced as an additional species, glucose transporter kin-
etics (r ¼ GlcTrans) and glucose flows in and out of
the reactor (r ¼ inGlc) were added, a glycogen-produ-
cing side branch was added (r ¼ storage), and the
removal of extracellular acetaldehyde (ACA) (r ¼ out-
ACA) was changed so that it is now formally com-
posed of the ACA leaving the reactor with the
outflow, and the ACA being removed by reactions in
the extracellular medium. This resulted in the 20L8D
model structure shown in Table 1 and Fig. 3. The cor-
responding ordinary differential equations (ODEs) are
constructed from Table 1 according to
y
s

and v
sr
is the stoichiometric coefficient of species s in
reaction r.
Neither the original 7D nor the new 8D model
structure have dynamic properties similar to the 20D
model (or the experiments) when the parameter values
of Wolf and Heinrich are inserted. We therefore per-
formed a parameter optimization in order to achieve
this. The 8D model structure has a limited number of
intrinsic parameters: K
trans
, q, K
i
, y
vol
and [Glc
x
]
0
.
(Intrinsic parameters are those that are not scalar mul-
tipliers of the rate expressions [30,35,36]). This allows
us to perform parameter optimization in an efficient
and unique manner, which we now explain. We first
parameterize the velocity parameters (i.e. the non-
intrinsic parameters) k
0
, V
1

0
) [Glc
x
])
GlcTrans: Glc
x
fi Glc V
1
½Glc
x

K
trans
þ½Glc
x

HK–PFK: Glc + 2 ATP fi 2 trioseP
+ 2 ADP
V
2
½Glc½ATP
1þ
½ATP
K
i

q
GAPDH: trioseP + NAD
+
fi BPG

9
([ACA] ) [ACA
x
])
outACA: ACA
x
fi (k
0
+ k
10
) [ACA
x
]
S. Danø et al. A case study in model reduction
FEBS Journal 273 (2006) 4862–4877 ª 2006 The Authors Journal compilation ª 2006 FEBS 4865
operating point and the intrinsic parameters [30,37].
For example, V
1
¼
m
GlcTrans
ðK
trans
þ½Glc
x
]Þ
½Glc
x

. We then fix the

as the bifurcation
parameter as in Hynne et al. [30]). Because we are left
with only three free parameters, and are constrained
by a 2D Hopf manifold, it is possible to obtain a com-
plete overview of the parameter space. This allows us
to choose a unique parameter set which results in opti-
mal agreement with the dynamic properties of the 20D
model. In this sense, the resulting 20L8D model is
unique. Details of the optimization are given in the
supplementary material (Doc. S1). The final set of
parameters (11 velocity parameters and seven intrinsic
parameters) is given in Table S1 and Table 6.
In the cause of the optimization we noticed a
remarkable property of the reduced model. When con-
strained by the operating point and the Hopf manifold,
the frequency of oscillation and the right critical eigen-
vector (i.e. the complex, right eigenvector associated
with the complex conjugate eigenvalues which have
zero real part) are both constant under the variation of
the remaining free parameters (K
trans
, q and K
i
). This
implies that the constraints and the model structure in
combination dictate the frequency of the oscillation as
well as the relative amplitudes and phases of the species
(properties of the right critical eigenvector).
The frequency of oscillation does, however, change
if the operating point is changed. The 20L8D model is

‘biochemical intuition’ to choose the relevant reduced-
model structure, although this need not be the case
[5,39,40]. We present ENVA as an alternative
approach to model reduction, where the system’s basic
dynamic properties are used as a guide for the elimin-
ation of variables that are not essential for the dynam-
ics. We eliminate a variable by fixing the metabolite
concentration at its steady-state value at a particular
operating point of the original model. A systematic
search is performed among the possible models in
order to identify the model of lowest dimension, which
retains the basic dynamic properties of the full system.
The basic dynamic properties of each of the reduced
Fig. 3. Reaction network of the 20L8D model. Extracellular species
and reactions are shown in red.
A case study in model reduction S. Danø et al.
4866 FEBS Journal 273 (2006) 4862–4877 ª 2006 The Authors Journal compilation ª 2006 FEBS
models are evaluated by calculating the eigenvalues of
the Jacobian matrix at a particular stationary state,
common to all models. In this case, the basic dynamic
property is oscillation. Models with an oscillatory
mode are identified as those with one or more sets of
complex conjugate eigenvalues, and oscillatory models
are identified as those with complex conjugate eigen-
values with positive real parts. This is standard nonlin-
ear dynamics theory [41].
Because we seek model(s) of the lowest dimension,
which retains the basic dynamic properties of the full
system, it is not necessary to search all 2
N

model. Its structure consists of the reactions involving
one or more of these species (Table 3, Fig. 4).
The eliminated variables must be represented in the
reduced model. This can be done in several ways. For
example, the quasi steady-state approximation can be
applied for each of the eliminated variables, or they
can simply be fixed at their steady-state values. In this
study we want the models to become as simple as poss-
ible, and we therefore take the last approach, which
results in simpler rate equations. When N ) n ¼ 14
species are fixed at their steady-state values, we intro-
duce 14 conservation-of-mass relations.
In order to ensure self-consistency, we must make
sure that all these conservation-of-mass relations are
fulfilled within the framework of the model. Some of
these relations are external to the reduced model
(Table 3) and do not call for any action. Others have
both internal and external parts. For example,
d½NADH
dt
¼ 0 ¼ v
glycerol
þ v
ADH
À v
GAPDH
has the external part v
ADH
and the internal parts
v

Complex
eigenvalues
11110000100000000000 )6.28 ± 4.98 i 11110000100000000001 )6.28 ± 4.98 i
11101100000000000000 )3.25 ± 9.86 i 11101110000000000000 )0.42 ± 7.53 i
11111100000000000000 )3.25 ± 9.86 i 11101101000000000000 )3.25 ± 9.86 i
11111000100000000000 )6.08 ± 4.42 i 11101100100000000000 )2.73 ± 9.57 i
11110100100000000000 )7.12 ± 5.32 i 11101100010000000000 )5.08 ± 10.6 i
11110001100000000000 0.95 ± 7.91 i 11101100001000000000 )3.25 ± 9.86 i
11110000110000000000 )7.23 ± 12.9 i 11101100000100000000 )3.27 ± 9.84 i
11110000101000000000 )6.28 ± 4.98 i 11101100000010000000 )3.25 ± 9.86 i
11110000100100000000 )6.28 ± 4.97 i 11101100000001000000 )3.25 ± 9.86 i
11110000100010000000 )6.28 ± 4.98 i 11101100000000100000 )3.25 ± 9.86 i
11110000100001000000 )6.28 ± 4.98 i 11101100000000010000 )3.25 ± 9.86 i
11110000100000100000 )6.28 ± 4.98 i 11101100000000001000 )3.25 ± 9.86 i
11110000100000010000 )6.28 ± 4.98 i 11101100000000000100 )3.25 ± 9.86 i
11110000100000001000 )6.28 ± 4.98 i 11101100000000000010 )3.25 ± 9.86 i
11110000100000000100 )6.28 ± 4.98 i 11101100000000000001 )3.25 ± 9.86 i
11110000100000000010 )6.28 ± 4.98 i 01110011010000000000 )168 ± 0.197 i
S. Danø et al. A case study in model reduction
FEBS Journal 273 (2006) 4862–4877 ª 2006 The Authors Journal compilation ª 2006 FEBS 4867
d½PEP
dt
¼ 0 ¼ v
lpPEP
À v
PK
ð3Þ
d½G6P
dt
¼

of v
HK
with v
storage
+v
PFK
. The resulting model is the
20E6D model defined by Table 3 and the rate expres-
sions in Table 4. The corresponding ODEs are con-
structed from the tables according to
dc
s
dt
¼
P
r
m
sr
v
r
.
Elimination of variables allowed us to lump a number
of parameters as indicated in Table 4. All the underly-
ing parameter values are the same as in the 20D model,
i.e. no parameter optimization was done in the elimin-
ation of variables approach. The model’s parameters
are listed in Table S2.
With a 20D model as the starting point it is feas-
ible to do a complete scan of all the possible reduced
models that can be constructed by elimination of vari-

x
]
0
¼ 24 mm. Of the reduced
models with complex eigenvalues, those of lowest
dimensionality are 3D; we find one with positive
real parts of the complex eigenvalues and one with
Table 3. The model structure of the 20E6D model. The stoichio-
metric constraint A
tot
¼ [ATP] + [ADP] + [AMP] reduces the dimen-
sion of the model to six.
Reaction r
HK
a
: ATP fi ADP
PFK: ATP fi ADP + FBP
ALD: FBP Ð GAP þ DHAP
TIM: DHAP Ð GAP
GAPDH: GAP Ð BPG
lpPEP ADP þ BPG Ð ATP
PK
b
: ADP Ð ATP
glycerol: DHAP fi
storage: ATP fi ADP
ATPase: ATP fi ADP
AK: ATP þ AMP Ð 2ADP
a
The rate expression of the hexokinase reaction has been substi-

As was the case with the construction of the 20E6D
model, a search within the subset of reduced models
defined by the ranking of the species according to their
decreasing importance, successfully identifies the oscil-
latory, reduced model of lowest dimension.
Model properties
We judge the effects of the model reductions by compar-
ing the dynamic and biochemical properties of the ori-
ginal 20D model to those of the three reduced models.
Table 4. Rate expressions of the 20E6D model. The reaction names r refer to Table 3. The model reduction allowed us to lump some of the
parameters (indicated by ~), but the underlying parameters are as in the parent 20D model. A list of the parameter values is given in Table
S2.
r Rate expression v
r
HK:
~
V
5m
~
K
5
þ
½ATP
½AMP

2
þ
~
k
22

6
½GAPK
6DHAP
þ½DHAPK
6GAP
þ½GAP½DHAPðÞ
TIM:
V
7m
½DHAPÀ
½GAP
K
7eq

K
7DHAP
þ½DHAPþ
K
7DHAP
½GAP
K
7GAP
GAPDH:
~
V
8m
½GAPÀ½BPG=
~
K
8eq

~
K
15
þ½DHAP
storage:
~
k
22
½ATP
ATPase: k
23
½ATP
AK: k
24f
½AMP½ATPÀk
24r
½ADP
2
Table 5. Model structure of the 20LE3D model. The stoichiometric
constraint A
tot
¼ [ATP] + [ADP] reduces the dimension of the
model to three. The corresponding set of ODEs are constructed
according to
dc
s
dt
¼
P
r

6
[ATP]
storage: 2ATP fi 2ADP
~
k
7
½ATP
glycerol: trioseP fi
~
k
8
½trioseP
Fig. 5. Reaction network of the 20LE3D model.
S. Danø et al. A case study in model reduction
FEBS Journal 273 (2006) 4862–4877 ª 2006 The Authors Journal compilation ª 2006 FEBS 4869
Stuart–Landau parameters and location of Hopf
bifurcations
The dynamic properties of the models are reflected by
their Stuart–Landau parameters (Eqn 1) at the Hopf
bifurcations found with [Glc
x
]
0
as bifurcation param-
eter. In cases where the [Glc
x
]
0
parameter has been
eliminated, we used the parameter that corresponds

Table 6. Stuart–Landau parameters of the models. The table lists
the Stuart–Landau parameters (Eqn 1) of the Hopf bifurcations
found on the borders of the oscillatory region. Re(r) and Im(r)
determine the rates at which the linear instability of the stationary
state and the frequency of oscillations at the stationary state,
respectively, increase with the bifurcation parameter l. The addi-
tional frequency change caused by amplitude changes is deter-
mined by Im(g) ⁄ Re(g). The bifurcation parameter l is ([Glc
x
]
0
)
[Glc
x
]
0
,
bif
) ⁄ [Glc
x
]
0
,
bif
in the 20D and the 8D models, but because
this parameter has been eliminated from the 20E6D and 20LE3D
models, we used l ¼ (C
G6P
) C
G6P,bif

)
Im(g) ⁄
Re(g)
20D [Glc
x
]
0,bif
¼ 18.5 mM 10 1.1 )2.1 1.4
20L8D [Glc
x
]
0,bif
¼ 18.5 mM 17 1.0 1.8 1.5
20E6D C
G6P,bif
¼ 11.4 lM 15 (0.18) (0.028) 2.4
20E6D C
G6P,bif
¼ 5.61 mM 8.9 ()27) ()6.3) 15
20LE3D C
Glc,bif
¼ 6.12 mM 18 (4.1) (8.1) 1.0
FB P
AT P
DHAP
F6P
ADP
GAP
G6 P
20D

the lower part of glycolysis via allosteric activation of
phosphofructokinase (PFK), and substrate for the
lower part of glycolysis inhibits low energy charge by
increasing ATP production in the phosphoglycerate
kinase (PGK) and pyruvate kinase (PK) reactions [31].
In conclusion, Fig. 6 shows that the reduced models
can be seen as depicting the oscillatory core of the full-
scale model. Because this is one of the main objectives,
this is a strong indication that the reduction proce-
dures are good choices for reduction to an oscillating
core.
Interaction analysis
The nature of the regulatory mechanisms underlying
the oscillations can be assessed by ranking the species
according to the importance of the feedback loops they
are involved in. We do this by employing a slightly
modified version of one of the methods presented in
Schmidt & Jacobsen [42]. The basis of the method is
the fact that the appearance of complex behaviour,
such as bistability and oscillations, can be traced back
to changes in the local stability properties of the sys-
tem’s steady state. In the case of autonomous oscilla-
tions, the underlying steady state is an unstable focus.
This is reflected by the Jacobian matrix of the steady
state, which has at least one pair of unstable conjugate
complex eigenvalues. For each species, there exists a
feedback loop conveying its effect on the other species
of the system. For each of these feedback loops, the
original method determines the minimal, real valued,
relative perturbation required for stabilization of the

the central part of glycolysis. The importance ranking
matches the polar phase plane plot analysis above, and
the conclusions of Madsen et al. [31]. This is partic-
ularly so when it is noted that the high importance
ranking of BPG is due to its very low average concen-
tration, which results in a high relative amplitude. In
broad terms, the ranking is conserved in the model
reductions. In the process of lumping and fitting which
leads to the 20L8D model, the relative importance of
BPG is increased, whereas the importance ranking of
the pooled species Glc and trioseP is in good agree-
ment with that of the corresponding species in the 20D
model (Fig. 7). Further reduction of the 20L8D model
to the 20LE3D model preserves the ranking. Reduc-
tion of the 20D model to the 20E6D model preserves
the relative ranking of ATP, ADP, BPG and FBP,
whereas GAP is now ranked more important than
DHAP (Fig. 7, lower left).
It is interesting to note that some of the variables in
the 20E6D and 20LE3D models do not have very high
importance values, even though these models were
constructed by the elimination of nonessential varia-
bles. We suggest that such variables are essential for
the connectivity of the network, rather than for the
generation of the oscillations per se.
Flux control
The flux-control pattern is an important biochemical
property. We determine this pattern using metabolic
control analysis as described previously [45]. The ana-
lyses are carried out at the (lower) Hopf bifurcation

from supply to demand control can readily be under-
stood as a consequence of the model reduction,
because the mechanical flow of the reactor, k
0
, has
Fig. 7. Comparison of species’ importance
rankings. The heights of the bars indicate
the importance of the feedback loops asso-
ciated with each of the species, as deter-
mined by interaction analysis. e
s
is the
smallest scalar perturbation of the linear
feedback of species s which causes the
unstable complex conjugate eigenvalues of
the Jacobian to disappear [42]. Hence, 1/|e|
is a measure of importance. A large value of
the importance measure indicates that the
stability of the system is very sensitive to
the feedback of that particular species.
Fig. 8. Comparison of flux-control patterns.
The heights of the bars indicate the magni-
tude of the flux control coefficients, and the
colour coding indicates the signs: black is
positive (reactions increasing the flux) and
white is negative (reactions decreasing the
flux). All flux-control patterns are calculated
at the (lower) Hopf bifurcations of the mod-
els (Table 6).
A case study in model reduction S. Danø et al.

) ⁄ p
0
. The scaled
(C) and unscaled (G) control coefficients are then
given by
C
x
lc
p
¼
d ln x
d ln p
¼
1
x
0
Imðr
p
ÞÀReðr
p
Þ
ImðgÞ
ReðgÞ

ð5Þ
for the frequency-control coefficient, and by
C
a
2
p

reactions HK, storage and ATPase are the major sta-
bilizing reactions. As expected [43], the frequency-
control pattern (Eqn 5) of the 20D model (Fig. 10,
upper left) involves more reactions than the stability-
control pattern, most notably the redox reactions
(GAPDH, ADH and the glycerol branch) and the spe-
cific flow rate of the reactor, k
0
. (The results for the
20D model have been published previously [31].) In
comparison, the stability-control pattern of the 20L8D
model (Fig. 9, upper right) reveals a less important
role for PFK in the destabilization of the stationary
state; the glucose transporter is now equally important.
Among the stabilizing reactions, GAPDH is now more
Fig. 9. Comparison of stability-control pat-
terns. The heights of the bars indicate the
magnitude of the unscaled stability control
coefficients C
ReðkÞ
p
¼ ReðrÞ (Eqn 6). Black
bars represent positive values (destabilizing
reactions), and white bars represent negat-
ive (stabilizing reactions). The sensitivity
analyses are made at the (lower) Hopf bifur-
cations of the models (Table 6).
S. Danø et al. A case study in model reduction
FEBS Journal 273 (2006) 4862–4877 ª 2006 The Authors Journal compilation ª 2006 FEBS 4873
important than in the 20D model. As in the 20D

The first method, LASCO, is based on lumping and
subsequent optimization. We have shown that this
optimization can be carried out in a very efficient man-
ner, by constraining the reduced model by the oper-
ating point (concentrations, fluxes, etc.) of the parent
model. This is an application of the direct method of
optimization [30,35,36]. For the chosen test case we
found that variations in the intrinsic parameters, which
are the only remaining free parameters, did not affect
the frequency and relative phases. We do not know for
how big a class of models this holds, but we want to
point out that this feature might be useful when sol-
ving the general problem of model rejection.
LASCO can be applied at any stationary state, sta-
ble or unstable. Comparison between the parent and
the reduced model can be based on any property of a
stationary state, e.g. control coefficients as determined
in metabolic control analysis [45] or elements of the
Jacobian matrix [46,47]. A bifurcation is not needed,
but, if present, it can be exploited for efficient compar-
ison of dynamic properties.
The second method, ENVA, performs a complete
search among all possible combinations of eliminated
variables. From these, the reduced model is picked as
the most highly reduced model, which retains the basic
dynamic features of the original (here, oscillations).
The search among candidate models constitutes a
potential combinatorial problem. We have shown how
this problem can be overcome by restricting the search
Fig. 10. Comparison of frequency-control

be advantageous.
We evaluated the two reduction methods by com-
paring the properties of the three reduced models and
the parent model. In short, we found that the dynamic
structures of the models are similar, but that their bio-
chemical properties are different.
That the dynamic structures are similar can be seen
from Fig. 6, which shows that the oscillatory modes of
the four models have the same biochemical composi-
tions. The central feedback mechanisms between these
modes are the allosteric regulation of PFK and posi-
tive stoichiometric feedback from the lower ATP-pro-
ducing steps to the upper ATP-consuming steps; it is
thus the same mechanism as reported in Madsen et al.
[31]. Because the reduction methods did not ensure the
preservation of these modes, but only the preservation
of oscillations, this study provides further support for
this mechanism. That this feedback can, in principle,
give rise to oscillations was shown several decades ago,
using minimal modelling [24–27,29]. The models pre-
sented here are a verification of these results, but this
study has the additional strength that the models were
obtained through the reduction of a realistic full-scale
model. Consequently, they have, for example, more
realistic parameter values, fluxes, and steady-state
concentrations.
That the biochemical properties of the models are
different can be seen from the sensitivity analyses of
frequency, stability and flux. For example, the flux-
control analysis reveals that flux through the 20E6D

In conclusion, we can thus say that both methods
have been successful in the fulfilment of the given
objective: to produce reduced, but biochemically mean-
ingful, models that reproduce the basic dynamic prop-
erties. The strength of ENVA is that it is algorithmic,
and that it does not require any input in the form
of biochemical knowledge. A major advantage of
LASCO, however, seems to be that it results in models
with more well-preserved biochemical properties. We
have shown these statements through extensive analy-
sis of the resulting models.
Acknowledgements
The work was supported by the Swedish Foundation
for Strategic Research and the European Commission
through the BioSim project (Contract LSHB-CT-2004–
005137), which are gratefully acknowledged.
References
1 Isidori I (1995) Nonlinear Control Systems, 3rd edn.
Springer, London.
2 Cedersund G, Danø S, Sørensen PG & Jirstrand M
(2005) From in vitro biochemistry to in vivo
understanding of glycolytic oscillations in Saccharo-
myces cerevisiae. Proceedings of the Conference on
Modeling and Simulation in Biology, Medicine and Bio-
medical Engineering. BioMedSim, 2005, Linko
¨
ping,
Sweden.
S. Danø et al. A case study in model reduction
FEBS Journal 273 (2006) 4862–4877 ª 2006 The Authors Journal compilation ª 2006 FEBS 4875

of low-dimensional manifolds in systems approaching
equilibrium. J Chem Phys 111, 859–847.
12 Glad T & Ljung L (2000) Control Theory: Multivariable
and Nonlinear Methods. Taylor & Francis, London.
13 Hahn J & Edgar TF (2002) An improved method for
nonlinear model reduction using balancing of empirical
gramians. Comput Chem Eng 26, 1379–1397.
14 Liebermeister W (2005) Dimension reduction by bal-
anced truncation applied to a model of glycolysis. Pro-
ceedings of the 4th Workshop on Computation of
Biochemical Pathways and Genetic Networks, pp. 21–28.
Logos-Verlag, Berlin
15 Liebermeister W, Baur U & Klipp E (2005) Biochemical
network models simplified by balanced truncation.
FEBS J 272, 4034–4043.
16 Betz A & Chance B (1965) Phase relationship of glyco-
lytic intermediates in yeast cells with oscillatory meta-
bolic control. Arch Biochem Biophys 109, 585–594.
17 Winfree AT (1972) Oscillatory glycolysis in yeast: the
pattern of phase resetting by oxygen. Arch Biochem
Biophys 149, 388–401.
18 Kreuzberg KH & Betz A (1979) Amplitude and period
length of yeast NADH oscillations fermenting on differ-
ent sugars in dependence of growth phase, starvation
and hexose concentration. J Interdis Cycle Res 10, 41–
50.
19 Richard P, Teusink B, Hemker MB, van Dam K &
Westerhoff HV (1996) Sustained oscillations in free-
energy state and hexose phosphates in yeast. Yeast 12,
731–740.

30 Hynne F, Danø S & Sørensen PG (2001) Full-scale
model of glycolysis in Saccharomyces cerevisiae. Biophys
Chem 94, 121–163.
31 Madsen MF, Danø S & Sørensen PG (2005) On the
mechanisms of glycolytic oscillations in yeast. FEBS J
272, 2648–2660.
32 Kuramoto Y (1984) Chemical Oscillations Waves and
Turbulence. Springer, New York.
33 Danø S, Hynne F, De Monte S, d’Ovidio F, Sørensen
PG & Westerhoff H (2001) Synchronization of glycoly-
tic oscillations in a yeast cell population. Faraday Dis-
cuss 120, 261–276.
34 Ipsen M, Hynne F & Sørensen PG (1998) Systematic
derivation of amplitude equations and normal forms for
dynamical systems. Chaos 8, 834–852.
35 Hynne F, Sørensen PG & Møller T (1993) Current and
eigenvector analyses of chemical reaction networks at
Hopf bifurcations. J Chem Phys 98, 211–218.
36 Hynne F, Sørensen PG & Møller T (1993) Complete
optimization of models of the Belousov–Zhabotinsky
A case study in model reduction S. Danø et al.
4876 FEBS Journal 273 (2006) 4862–4877 ª 2006 The Authors Journal compilation ª 2006 FEBS
reaction at a Hopf bifurcation. J Chem Phys 98, 219–
230.
37 Brusch L, Cuniberti G & Bertau M (2004) Model eva-
luation for glycolytic oscillations in yeast biotransforma-
tions of xenobiotics. Biophys Chem 109, 413–426.
38 Kohout M, Schreiber I & Kubı
´
e

47 Mihaliuk E, Skødt H, Hynne F, Sørensen PG & Sho-
walter K (1999) Normal modes for chemical reactions
from time series analysis. J Phys Chem 103, 8246–8251.
48 Henson MA, Mu
¨
ller D & Reuss M (2002) Cell popula-
tion modelling of yeast glycolytic oscillations. Biochem J
368, 433–446.
Supplementary material
The following supplementary material is available
online:
Doc. S1. Parameter optimization of the constrained
20L8D model.
Table S1. Parameters of the 20L8D model.
Table S2. Parameters of the 20E6D model.
Table S3. Parameters of the 20LE3D model.
Table S4. Unstable stationary state of the 20L8D
model provided as a check of model implementations.
Table S5. Unstable stationary state of the 20E6D
model provided as a check of model implementations.
Table S6. Unstable stationary state of the 20LE3D
model provided as a check of model implementa-
tions.
This material is available as part of the online article
from
S. Danø et al. A case study in model reduction
FEBS Journal 273 (2006) 4862–4877 ª 2006 The Authors Journal compilation ª 2006 FEBS 4877