## Abstract

Explicit modelling of metabolic networks relies on well-known mathematical tools and specialized computer programs. However, identifying and estimating the values of the very numerous enzyme parameters inherent to the models remain a tedious and difficult task, and the rate equations of the reactions are usually not known in sufficient detail. A way to circumvent this problem is to use ‘non-mechanistic’ models, which may account for the behaviour of the systems with a limited number of parameters. Working on the first part of glycolysis reconstituted *in vitro*, we showed how to derive, from titration experiments, values of effective enzyme activity parameters that do not include explicitly any of the classical kinetic constants. With a maximum of only two parameters per enzyme, this approach produced very good estimates for the flux values, and enabled us to determine the optimization conditions of the system, i.e. to calculate the set of enzyme concentrations that maximizes the flux. This fast and easy method should be valuable in the context of integrative biology or for metabolic engineering, where the challenge is to deal with the dramatic increase in the number of parameters when the systems become complex.

- biochemical pathway
- flux computation
- glycolysis
- hyperbolic fitting
- optimal enzyme concentration
- systemic enzyme parameter

## INTRODUCTION

The ‘omics’ approaches are challenging biochemical practice. The complete genome sequence data and microarray technologies reveal frequent occurrences of gene families with complex patterns of gene regulations, quantitative proteomics give expression data for hundreds of enzymes and isoenzymes, and metabolomics techniques give access to thousands of major metabolites. The way all these systemic variables are affected by development, environmental conditions and/or genetic background gives invaluable information to investigate the nature, organization and control of complex metabolic networks. Modelling metabolic systems relies on classical theoretical methods, such as differential equations, and on powerful computer programs. With these tools, extensive efforts have been made to achieve a comprehensive dynamic description of some essential pathways, such as glycolysis, especially in yeast [1,2], *Trypanosoma brucei* [3], erythrocytes [4–6] and muscle [7]. However, such approaches, which require a long-term effort of numerous laboratories for simple pathways, have severe limitations for more complex systems. The main obstacle for biochemical modelling is usually not the availability of mathematical and informatics tools, but the determination of the parameter values. Single rate equations may include more than ten parameters {as arises in the equation for PFK (6-phosphofructokinase; EC 2.7.1.11) [2]}, and the number of parameters increases dramatically with the size and complexity of the systems. Except for the concentration of the most abundant enzymes, which can be estimated with quantitative proteomics (e.g. [8]), there is no simple or direct way to measure essential parameters such as catalytic constants, Michaelis–Menten constants, or inhibition constants. For example, the robot-based platform recently developed to measure activities of approx. 20 enzymes in *Arabidopsis* is promising, but does not provide values for such essential parameters [9]. *In vivo*, it is questionable whether all the parameters can be estimated from perturbation experiments [10], in particular because the range of variation of metabolite levels is limited by homoeostasis constraints [11]. For some well-described pathways, computer simulations can be performed to estimate the parameter values by fitting the model with *in vivo* observations [12]. In any case, estimating the parameter values remains a tedious and time-consuming task, with multiple factors of inaccuracy, and no assurance that all the key parameters are identified.

Consequently, modelling efforts based on conceptual short-cuts are essential to simulate complex cellular behaviours from a limited amount of biological data. Flux balance analysis (see [13] for a review) represents a radical approach in this direction, since no kinetic information is included in the models. Flux balance analysis does not attempt to predict the exact behaviour of metabolic networks, but uses known constraints on the integrated function of multiple enzymes to separate the state that a system can reach from those it cannot. Constraints considered are the steady-state mass balance of metabolites, the effective irreversibility of reactions due to thermodynamic constraints, or the maximum flux capacity of enzymes or transport proteins. In the MCA (metabolic control analysis) [14,15] or in the BST (biochemical systems theory) [16], rate laws are not explicitly written, but kinetic parameters are not ignored. Both theoretical frameworks use systemic coefficients (control coefficients or logarithmic gains) and local coefficients (elasticity or kinetic order), which can be estimated experimentally, to quantify the behaviour of the systems. Derived from MCA or BST, various so-called ‘non-mechanistic’ models have been proposed that aim at using approximate equations to describe the behaviour of the systems with a reduced number of parameters. The power law representations, such as the S-systems, are obtained by specifying a rate constant (the multiplicative term) and writing one power term for each variable that directly influences the rate law (metabolites, enzymes and effectors) [17]. The (log)linear kinetic models rely on knowledge about the strength of interaction among the rates of enzyme-catalysed reactions and the concentrations of the various metabolites and regulators of the network [18]. The linlog kinetics proposes general expressions giving steady-state fluxes and metabolite levels as a function of enzyme concentration and response coefficients [11]. Finally, the ‘tendency modelling’ approach results in rate equations containing a minimal number of parameters, because it assumes that the behaviour of the metabolites is determined more by the stoichiometric structure of the network than by the exact enzyme mechanisms [19].

The present study is in the line of the modelling approaches using a limited number of parameters. We show how to use a simplified formalism based on MCA to derive global parameters that account for the kinetic behaviour of enzymes in their metabolic context. According to the structure of the pathway and to the position of the enzyme in the pathway, one or two parameters per enzyme are sufficient. We evaluated the validity of the method on a pathway reconstituted *in vitro*, an approach that just proved to be very efficient to study the systemic properties of the pathways, such as elasticities, flux control coefficients and transition time control coefficients [20–22], dynamic properties [23,24] and allosteric controls [25]. The upstream part of glycolysis, from HK (hexokinase; EC 2.7.1.1) to G3*P*DH (glycerol-3-phosphate dehydrogenase; EC 1.1.1.8), allowed us to measure the flux for a large number of different distributions of enzyme concentrations, and to compare the values with the theoretical predictions. We showed that the predicted flux values were tightly correlated to the measured flux values. In addition we derived equations to calculate the optimal distribution of enzyme concentrations for a given total enzyme content allocated to the network. Again the predicted values were fully consistent with observations. As the method does not require the estimation of any classical kinetic parameter such as *K*_{m}, *V*_{max} or *K*_{i}, the experimental effort for modelling and determining the conditions for optimization is thus significantly reduced. Application of the approach to metabolic engineering and to *in vivo* measurements is discussed.

## THEORETICAL DEVELOPMENTS

### Linear pathway of Michaelis–Menten enzymes

The flux *J* through a linear pathway of *n* unimolecular reversible reactions catalysed by Michaelis–Menten enzymes far from saturation is [14,15]:
(1)
where *S*_{1} and *S*_{n+1} are respectively the concentrations of the initial and final substrates of the pathway, *K*_{1,n+1} is the product of the equilibrium constants of all the reactions, *K*_{1,j} is the product of the equilibrium constants from reaction 1 to reaction *S*_{j−1}→*S*_{j}, *K*_{mj} is the Michaelis–Menten constant of enzyme *j* and *V*_{j} is the maximum velocity of enzyme *j*. Since the maximum velocity is given as *V*_{j}=*E*_{j}*k*_{cat j}, with *E*_{j} the concentration of enzyme *j* and *k*_{catj} its catalytic constant, we have:
(2)
where *S*=[*S*_{1}]−[*S*_{n+1}]/*K*_{1,n+1}.

For a given total enzyme concentration,
there is an optimal set of enzyme concentrations {*E*_{1}*, *E*_{2}*,…, *E*_{n}*}, that is, a distribution that maximizes the flux. At maximum flux, we have for all enzyme *i* [26,27]:
(3)
where *A*_{j}=*k*_{catj}*K*_{1,j}/*K*_{mj}.

### Estimating systemic kinetic parameters

Equations (2) and (3) do not apply to pathways such as glycolysis that contain non-Michaelian–Menten and regulated enzymes. Nevertheless, we wanted to evaluate the extent to which a similar formalism could help in computing estimates of flux values and determining the conditions for optimizing metabolic systems of any complexity.

Consider a system converting a metabolite *S*_{1} into *S*_{m+1} through a series of reactions catalysed by *n* reversible enzymes far from saturation. As branching is possible, *n*≥*m*. With a single input (*S*_{1}) and output (*S*_{m+1}), the flux through the whole pathway at steady state is the rate of consumption of *S*_{1}, which is equal to the rate of formation of *S*_{m+1}. On increasing the concentration *E*_{i} of enzyme *i*, the concentrations of other enzymes being fixed, the response of the flux is expected to display a quasi-hyperbolic ascending curve (saturation curve). This intuitive prediction is corroborated by an abundant theoretical and experimental literature (see the Discussion section). However, depending on the topology, the regulations and the constraints in the network, various mechanisms could result – at least theoretically – in a decrease in the flux when the concentration of certain enzymes increases (some examples are presented in the Discussion section). Moreover, non-linear behaviour, such as oscillations, can prevent a stable steady state being observed.

In the following developments, we considered only stable steady states, but for more generality we considered that the asymptotic response of the flux upon enzyme increase may be ascending or descending.

For fitting to experimental curves, we chose an expression of the flux *J* derived from eqn (2):
(4)
where *S* and *E*_{j} have the same definition as above. The parameter *A*_{j} is a complex systemic kinetic parameter that depends on the type of enzyme and regulation, and takes the above value (*A*_{j}=*k*_{catj}*K*_{1},_{j}/*K*_{mj}) only for Michaelis–Menten non-regulated enzymes in a linear pathway. Note that *A*_{j} may be negative, for example if the enzyme belongs to the backward branch of a substrate cycle. The dimension of *A*_{j} is μM^{−1}·s^{−1}. The parameter *p*_{j} accounts for the ‘dispensability’ of enzyme *j*: if *p*_{j}=0, then *J*=0 when *E*_{j}=0, and if |*p*_{j}|≠0, then *J*≠0 when *E*_{j}=0. For example, in the system considered (see Scheme 1 in the Experimental section), removing the TPI (triose-phosphate isomerase; EC 5.3.1.1) does not bring the flux to zero because of the branching of the pathway at the reaction catalysed by the FBA (aldolase; EC 4.1.2.13). So the *p*_{TPI} value is not null, unlike those of the other enzymes. The dimension of *p*_{j} is s^{−1}.

All those parameters can be estimated with titration experiments, by varying in turn the concentration of each enzyme, the other concentrations being kept constant. For enzyme *i*, we have the following expression of the flux:
(5)
When *E*_{i} tends to infinity, the flux tends to the value:
(6)
Thus the flux may be given as:
(7)
The three parameters *SA*_{i}, *Sp*_{i} and *J*_{∞ i} of this function of *E*_{i} can be estimated by hyperbolic fitting of the titration curves [we used the *nls* package of the version 2.0.1 of the R software (http://www.r-project.org/)]. Note that the ordinate of the intercept for *E*_{i}=0, when *Sp*_{i}≠0, is:
(8)
Hence:
(9)
For indispensable enzymes, we have *Sp*_{i}=0, and the flux value tends to 0 when *E*_{i} tends to 0. With all the *SA*_{i} and *Sp*_{i} estimates, the flux can be computed from eqn (4) for any set of *E*_{i} values.

### Determination of the optimal distribution of enzyme concentrations

For a given total enzyme concentration, the optimal distribution of the concentrations {*E*_{1}*, *E*_{2}*,…, *E*_{n}*} for a flux expressed by eqn (4) can be calculated using the Lagrange multiplier method (see the Appendix). Writing *e*_{i}*=*E*_{i}*/*E*_{tot}, we get:
(10)

This equation is valid only if *E*_{tot} corresponds to a fixed total molarity. If *E*_{tot} is a fixed total weight of enzymes (e.g. in g/l), as considered in the present study, the molecular masses of enzymes have to be included in the optimization equation, as shown in the Appendix. We get:
(11)
where *A*′_{j}=*A*_{j}/*M*_{j}, with *M*_{j} the molecular mass of enzyme *j*, and where the left subscript ‘g’ means that the enzyme concentrations are in g/l.

## EXPERIMENTAL

### Chemicals

D-Glucose, Pipes, ATP, NADH, phosphocreatine and all enzymes were obtained from Sigma. HK, PGI (phosphoglucoisomerase; EC 5.3.1.9) and TPI were purified from baker's yeast. FBA, G3*P*DH and CK (creatine kinase; EC 2.7.3.2) were purified from rabbit muscle, and PFK was from *Bacillus stearothermophilus*. All the purchased commercial enzymes were crystalline suspensions in ammonium sulphate. These proteins were desalted by using PD-10 desalting columns according to the manufacturer's recommendations (Amersham Biosciences). The protein concentrations were determined at λ=205 nm [28]. The molecular masses of the enzymes are: PGI, 61300 Da; PFK, 34000 Da; FBA, 39212 Da; and TPI, 26795 Da.

### Measurements of kinetic parameters of isolated enzymes

All enzyme assays were carried out in a thermoregulated quartz cuvette (25 °C), in a buffer containing 50 mM Pipes (pH 7.5), 100 mM KCl and 5 mM magnesium acetate. HK activity was measured with 1 μM lactate dehydrogenase, 1 μM pyruvate kinase, 1 mM phosphoenolpyruvate and various concentrations of glucose. PGI activity was measured in the forward direction with 0.5 μM PFK, 5 μM FBA, 3 μM G3*P*DH and glucose 6-phosphate as substrate. PFK activity was measured with 0.25 μM FBA, 0.41 μM TPI, 0.9 μM G3*P*DH and various concentrations of fructose 6-phosphate. FBA activity was measured with 0.5 mM TPI and 1 mM G3*P*DH and fructose 1,6-bisphosphate as substrate. TPI activity was assayed with 1 μM G3*P*DH and glyceraldehyde 3-phosphate as substrate.

For estimating the kinetic parameters of isolated enzymes, we used the classical Eadie–Hofstee representation [29] where the maximal velocity *V* and the Michaelis–Menten constant *K*_{m} were estimated by fitting the linear regression model
where *Y*_{i}=*v*_{i}/*S*_{i} are experimental data points for the reaction rate *v*_{i} at substrate concentration *S*_{i}; *X*_{i}=*v*_{i}; and ϵ_{i} is a random error term. In this representation, *K*_{m}=1/β and *V*=α/β. Unbiased estimators *K*_{m} and *V*, and their corresponding sampling variance, can be found by a Taylor series expansion of the ratios around the true values α and β (see for example [30]):
where

and

are estimated from the linear regression of *v*_{i}/*S*_{i} on *v*_{i}.

*In vitro* pathway and flux measurement

The system constructed *in vitro* was operated following Scheme 1.

The ATP pool is kept constant via the regeneration reaction catalysed by CK. HK is saturated due to the high glucose concentration (100 mM). The HK concentration was fixed in all experiments at 0.1 μM so that the glucose 6-phosphate input is the same in all the experiments. G3*P*DH, catalysing the oxidation of NADH, was kept constant, at high concentration (1 μM), in order to use the slope of the linear NADH decay as a measurement of the steady-state flux through the whole pathway. NADH consumption was monitored every 2 s with an Uvikon 850 spectrophotometer at 390 nm from 60 to 120 s. The assays were carried out at 25 °C in a buffer containing 50 mM Pipes (pH 7.5), 100 mM KCl, 20 mM phosphocreatine, 3 mM NADH and 5 mM magnesium acetate. Each reaction was started by addition of 1 mM ATP.

### Titration curves

We varied in turn the concentrations of each enzyme, the concentration of the others being maintained at reference values computed on the basis of their physiological concentration estimated in the yeast strain S288C [8]. The reference values are 0.15 μM for PGI, 0.29 μM for PFK, 1.54 μM for FBA and 0.84 μM for TPI (respectively 9.1, 10.4, 60.1 and 22.3 mg/l). Hyperbolic fitting of the curves gave the relevant enzyme parameters that can be used for computing the predicted flux values (see the Theoretical developments subsection).

### Selecting distributions of enzyme concentrations for flux measurement

To assess the reliability of the predicted flux values, and the validity of the model giving the optimal enzyme concentrations, we measured the flux values for a series of 121 distributions of enzyme concentrations, the total amount being kept constant. The sum of PGI, PFK, FBA and TPI concentrations was fixed at *E*_{tot}=101.9 mg/l, chosen as the sum of the reference concentrations. In order to cover a large range of variation, we varied the FBA concentration (the most abundant enzyme) from 0 to *E*_{tot}, taking ten values evenly distributed, and the proportions of the remaining enzymes (PGI, PFK and TPI) were drawn in β distributions with shape parameter α=1 and scale parameter
with *e*^{φ}_{i} the reference proportion of enzyme *i*. In the 121 distributions obtained, enzyme concentrations ranged from 0 to 70 mg/l for PGI and PFK and from 1.66 to 66.9 mg/l for TPI.

Each flux measurement was repeated three times in three independent tubes. As the fluxes should be expressed in μM/s, the enzyme concentrations (in mg/l) were converted into μM using the molecular masses given by the manufacturer.

## RESULTS

### Estimation of the kinetic parameters of isolated enzymes

The kinetics of the commercial enzymes in the forward direction was studied under the same conditions as those used for the system constructed *in vitro*. The values of the enzyme kinetic parameters estimated using the Eadie–Hofstee representation were of the same order of magnitude as those indicated by the manufacturer (Table 1).

### Titration experiments in the pathway constructed *in vitro*

Relationships between the flux and the concentration of each enzyme of the upstream part of glycolysis constructed *in vitro*, and the curves of the hyperbolic fittings, are shown in Figure 1. For each enzyme, the parameters *J*_{∞ i}, *SA*_{i} and *Sp*_{i} were estimated as indicated in the Theoretical developments subsection (Table 2). As expected, *SA*_{i} values were poorly correlated with the ratio *k*_{cat i}*K*_{1,i}/*K*_{m i} that can be estimated from isolated enzymes (results not shown). TPI alone is ‘dispensable’, since the flux with no TPI is estimated to be *J*_{0 TPI}=10.47 μM/s, whereas *J*_{0 i} for other enzymes is zero. So we have *Sp*_{TPI}=61.77 μM/s, and *Sp*_{i}=0 for PGI, PFK and FBA.

### Comparing predicted and measured flux values

With the estimates of the *SA*_{i} and *Sp*_{i}, we calculated a predicted flux value, *J*_{pred}, for any set of enzyme concentrations from eqn (4):
This predictor was used to compute the flux values for 121 distributions of enzyme concentrations, and these values were compared with the flux values measured *in vitro*, *J*_{obs}. Table 3 gives the *J*_{pred} and the *J*_{obs} values and the S.D. of the *J*_{obs}. The cv (coefficients of variation) never exceeded 0.14, and 84% of the cv were below 0.05. The linear correlation coefficient between the predicted and measured flux value was high (*r*=0.93) and the intercept of the regression line was not statistically different from zero (Figure 2). The slope of the regression line was statistically higher than 1 (*a*=1.38±0.18), which means that the predictor overestimates the observed flux by a constant factor. Thus *J*_{pred}/*a* appears to be a reliable predictor of the flux in the *in vitro* system.

### Optimizing the metabolic system

As the total enzyme concentration was fixed, the relationship between the concentration of a particular enzyme and the flux can no longer be hyperbolic, since there is a global negative correlation between the concentrations. High concentration of one enzyme reduces that of the others, which may decrease the flux. Consistent with this prediction, plotting the observed flux as a function of the concentration of any one enzyme showed that the experimental points were distributed below a hump-shaped curve (Figure 3). For a given concentration of one of the enzymes, the flux may vary between 0 and a local maximum, depending on the distribution of the concentrations among the other enzymes. The distributions were randomly drawn, which explains the high extent of scatter of the points (some enzyme concentrations may be very low by chance).

The optimal proportion of each enzyme in the system was computed using eqn (11), which gave the proportions 0.155 for PGI, 0.240 for PFK, 0.590 for FBA and 0.015 for TPI. As the total enzyme concentration, *E*_{tot}, was 101.9 mg/l, these values corresponded respectively to 15.8, 24.5, 60.1 and 1.5 mg/l. As expected, the flux value computed with these concentrations (*J*_{pred max}/*a*=14.62 μM/s) is higher than all the predicted or measured flux values, and the position of the optimal concentration for each enzyme corresponds to the top of the hump-shaped curve that ‘envelops’ all the points (Figure 3). For an enzyme concentration lower than the optimum, the local maximum flux increases with the enzyme concentration. For an enzyme concentration that exceeds the optimum, the local maximum flux decreases with the enzyme concentration.

## DISCUSSION

We used the reconstituted upstream part of glycolysis to evaluate a method for estimating *effective* enzyme parameters that account for the behaviour of the enzyme within its metabolic pathway. The flux values predicted from these parameters proved to be quite close to the measured flux values. The systemic parameter *A*_{i} that we defined includes all the kinetic and regulation properties of the enzyme, but its estimation does not require any measurement of parameters such as *K*_{m}, *k*_{cat} or *K*_{i}. Thus the kinetic type of the enzyme does not interfere with the modelling. For instance, PFK, which is a regulated enzyme, was treated like PGI, a Michaelian enzyme. It is worth noting that *A*_{i} also contains information on the position of the catalysed reaction in the pathway, via the equilibrium constants. In linear pathways of Michaelis–Menten enzymes far from saturation, *A*_{i} includes merely the product of the equilibrium constants from the initial reaction to the reaction catalysed by enzyme *i* [14], but in the general case of complex networks there is no simple relationship between the equilibrium constants. The other parameter that we defined is the ‘dispensability’ (*p*_{i}), which is different from zero whenever removing the enzyme from the system does not drive the flux to 0, a situation that occurs in the case of reticulation in the system. The parameters *A*_{i} and *p*_{i} could not be estimated independently of *S*, a parameter that depends also on the equilibrium constants of the reactions and on the initial and final substrates of the pathway, but this had no incidence on the calculus of the predicted flux.

The visual inspection of the titration curves used to estimate these parameters revealed that for every enzyme, the experimental points follow a slightly more convex curve than the fitted hyperbola, so that the estimated *J*_{∞ i} are higher than the apparent maximum experimental flux values. Thus the slope *a* of the regression line of *J*_{pred} on *J*_{obs} is higher than 1, so the right predictor of flux is *J*_{pred}/*a*. Other, more complex, functions could lead to better fittings. However, the quality of the prediction that we obtained with only two parameters does not justify making the model more complex. In addition, we think that the prediction could be even improved in a system in which the concentration of the first substrate could be maintained constant. Under our conditions, the first substrate of the system was glucose 6-phosphate; its input was constant because HK was saturated by glucose and did not vary, but its steady-state concentration depended on the concentration of the other enzymes. As a consequence, unlike what was assumed in the model, the parameter *S* could slightly vary, which inevitably increased the dispersion of the values.

The method for estimating global parameters can apply whether flux response towards enzyme variation is positive or negative, i.e. whether enzymes have positive or negative flux control coefficients. Here, the flux displayed a saturation curve for all enzymes, which means that the enzymes have a positive control coefficient for any concentration. This situation has been described for diverse pathways in many classical papers (e.g. [31–33]), and is theoretically expected in case of linear pathways [14] even if there are feedback or feedforward loops [34]. In this connection, it would be interesting to test the robustness of the method on the complete glycolysis, which would introduce strong feedback loops. In branched pathways, there are some particular situations where increasing an enzyme concentration could result in a decrease in the flux. (i) When the enzyme under study is located inside a feedback loop, negative control coefficients may be observed, depending on the parameter values of the network [34]. This is consistent with observations showing that increasing the trehalose-6-phosphate synthase may result in a decrease in trehalose synthesis, via the inhibition of HK by trehalose 6-phosphate [35]. (ii) When there is a substrate cycle (or ‘futile’ cycle), i.e. when a reaction, or a set of reactions, that converts metabolite *S*_{i} into *S*_{j}, is opposed by a second set that reconverts *S*_{j} into *S*_{i}, as observed for example with the cycle glucose–glucose 6-phosphate. Overexpression of glucose-6-phosphatase in rat hepatocytes decreases glucose 6-phosphate levels, with a negative control coefficient of this enzyme on glycogen synthesis [36]. In such cycles, increasing concentration of a backward enzyme is expected to decrease the whole flux to a minimal level, whereas a maximum flux value will be observed if the concentration of any one backward enzyme is null. (iii) When there is regulation by covalent modifications of enzymes, negative control coefficient may be observed. In the classical example of reversible phosphorylation of enzymes by protein kinases, the effects on the target enzymes can be activatory or inhibitory. Interestingly, when antagonistic enzymes in cycles are phosphorylated, usually one is activated and the other is inhibited [37]. Analogous, though more complex, mechanisms account for the high negative control coefficient of glycogen phosphorylase on glycogen synthesis [38]. All these situations necessarily result in negative values for *SA*_{i} and *Sp*_{i}, which is allowed in our modelling since we made no assumption on the sign of those parameters. (It is worth noting that when the enzyme under study does not belong to the branch the flux of which is measured, negative flux control coefficients are expected [39,40], but this situation is not relevant to our approach, because we considered the flux through the whole pathway, with input flux equals output flux.)

The equations that we derived also allowed us to find the optimal distribution of enzyme concentrations in case of constraint on total enzyme. Despite the fact that the amount of enzyme the cell allocates to a given metabolic pathway is necessarily limited, biochemical modelling implicitly considers in most cases that enzyme concentrations are independent of each other, i.e. there is no limitation of resources (matter and energy) in the cell. However, some results clearly indicate that there are constraints on the variation of enzyme concentrations, particularly for the most abundant ones, due to competition for space, energy and cellular resources. Thus global negative correlations may occur between enzyme concentrations, resulting in the so-called ‘protein burden effect’ [41]. For example, it has been shown that increasing the expression of various glycolytic enzymes in *Zymomonas mobilis* causes a significant decrease in the glycolytic flux and growth rate [42]. In terms of economy of resources, and hence from the evolutionary point of view, it is important to achieve a given flux value with the minimum amount of total enzyme. Formally, this question is equivalent to the one that we addressed, namely how to achieve the maximum flux value with a given amount of total enzyme (for a general treatment of this question, see [26,27,43]). The equations that we derived allow this optimal distribution to be found.

Another noticeable interest of predicting optimal distribution of enzyme concentrations is for *in vitro* metabolic engineering. For instance, a process for the production of 2-oxo-L-gulonic acid, a precursor to vitamin C, implies two enzymatic reactions exterior to the non-living permeabilized *Pantoea citrea* cells used in the bioreactor [44]. The optimal concentrations of those extra enzymes could be determined following the method described in the present paper. Similarly, the approach could be applied to functional protein chips based on mRNA–enzyme fusion molecules hybridized on a DNA microarray, where immobilized enzymes retain activity proportional to the amount of capture DNA, allowing modulation of the relative activity of enzymes; this novel and elegant technique was used for the optimization of the trehalose synthesis pathway [35], but on an empirical base, while our theoretical treatment would make it rational.

The approach is not restricted to systems reconstituted *in vitro*, and can be applied to two experimental systems closer to *in vivo* conditions. With cell-free extracts, titration experiments using external enzymes are possible, as has been shown for characterizing control properties in glycolysis of mouse muscle [45] and rat liver [46]. Interestingly, it was shown that the kinetic parameters of commercial enzymes did not produce a reliable model of control properties of HK and PFK as compared with kinetic parameters determined in muscle extract [45]. In plants, isolated mesophyll cells from leaves of *Digitaria sanguinalis* keep functional plasmodesmata that allow the free exchange of low-molecular-mass compounds with the culture medium. This system has been used successfully to analyse the kinetic properties of the phosphoenolpyruvate carboxylase, with external malate dehydrogenase added to the medium [47]. Varying external enzymes in such experimental systems would allow the parameters we defined to be estimated ‘in cell’.

In conclusion, we developed a general method that allows us both to predict flux values from global enzyme parameters and to estimate optimal enzyme concentrations. This method is promising for integrative biology and could find applications in metabolic engineering.

## Acknowledgments

We thank Professor Hans V. Westerhoff (Faculty of Earth and Life Sciences, Department of Molecular Cell Physiology, Vrije Universiteit Amsterdam, Amsterdam, The Netherlands) for fruitful discussions, Professor David Fell (School of Biological and Molecular Sciences, Oxford Brookes University, Oxford, U.K.) for useful suggestions and Professor Olivier C. Martin (Laboratoire de Physique Théorique et Modèles Statistiques, Université Paris-Sud, Orsay Cedex, France) and Dr Delphine Sicard (UMR de Génétique Végétale, INRA/UPS/CNRS/INAPG, Gif-sur-Yvette, France) for careful reading of this paper and relevant remarks. J.B.F. was supported by a Thesis grant from the French Ministère de la Jeunesse, de l'Education Nationale et de la Recherche.

## APPENDIX

### Calculation of the optimal concentrations of the enzymes

We used the Lagrange multiplier method to find the distribution of enzyme concentrations that maximizes the flux for a given amount of total enzyme. At maximum flux, we have: (A1)

Taking the expression of the flux of the main text (eqn 4), we have:
(A2)
where *E**_{i} is the optimal concentration of enzyme *i*.

From eqn (A1), we see that all the derivatives are equal to λ, so we have: (A3) which is equivalent to: (A4) Hence: (A5)

Introducing *E*_{tot}, it yields:
(A6)

After rearrangements, we get: (A7) This equation leads to eqn (10) of the main text, with

This general treatment, where implicitly the concentrations are expressed in mol/l, has to be adapted if the constraint on *E*_{tot} is a fixed total weight of enzymes, as we did in the present study. In this case, the enzyme concentration, _{g}*E*_{i} (left subscript ‘g’ for g/l), can be written as:
(A8)
where *M*_{i} is the molecular mass of enzyme *i*. Thus the constraint is:
(A9)
and eqn (A1) becomes:
(A10)
Here the derivatives of the flux are not equal, since we have:
(A11)
The flux equation can be rewritten to make the molecular masses apparent:
(A12)
with A′_{j}=*A*_{j}/*M*_{j}.
Hence eqn (A11) can be written as:
(A13)
with _{g}*E**_{i} the optimal concentration of enzyme *i*.

Thus replacing *A*_{i} (respectively *A*_{j}) by *A*′_{i} (respectively *A*′_{j}), and *E**_{i} (respectively *E**_{j}) by _{g}*E**_{i} (respectively _{g}*E**_{j}) in eqns (3)–(7), we obtain directly:
(A14)
This equation leads to eqn (11) of the main text, with

**Abbreviations:**
BST, biochemical systems theory;
CK, creatine kinase;
cv, coefficients of variation;
FBA, aldolase;
G3PDH, glycerol-3-phosphate dehydrogenase;
HK, hexokinase;
MCA, metabolic control analysis;
PFK, 6-phosphofructokinase;
PGI, phosphoglucoisomerase;
TPI, triose-phosphate isomerase

- The Biochemical Society, London