Dynamic modeling of complex biological systems: a link between metabolic and macroscopic description

Mathematical Biosciences 193 (2005) 25–49 www.elsevier.com/locate/mbs

Dynamic modeling of complex biological systems: a link between metabolic and macroscopic description Jens E. Haag a, Alain Vande Wouwer

a,*

, Philippe Bogaerts

b,1

a b

Service d’Automatique, Faculte´ Polytechnique de Mons, 31 Boulevard Dolez, 7000 Mons, Belgium Service de Chimie ge´ne´rale et Biosyste`mes, Universite´ Libre de Bruxelles, 50 Avenue, F.D. Roosevelt, CP 165/61, Brussels 1050, Belgium Available online 28 December 2004

Abstract In this study, a class of dynamic models based on metabolic reaction pathways is analyzed, showing that systems with complex intracellular reaction networks can be represented by macroscopic reactions relating extracellular components only. Based on rigorous assumptions, the model reduction procedure is systematic and allows an equivalent Ôinput–outputÕ representation of the system to be derived. The procedure is illustrated with a few examples.
2004 Elsevier Inc. All rights reserved. Keywords: Biotechnology; Mathematical modeling; Model reduction; Parameter estimation; Non-linear systems

1. Introduction In the last few decades, much effort has been put toward the mathematical description of mammalian cell cultures. Early attempts have revealed the cells as biocatalysts that reproduce themselves, excrete metabolites and often overexpress proteins of interest, e.g. for pharmaceutical purposes. Simple models *

Corresponding author. Tel.: +32 65 37 4141; fax: +32 65 37 4136. E-mail address: [email protected] (A.V. Wouwer). 1 Tel.: +32 2 650 4076; fax: +32 2 650 3575.

0025-5564/$ - see front matter
2004 Elsevier Inc. All rights reserved. doi:10.1016/j.mbs.2004.11.007

26

J.E. Haag et al. / Mathematical Biosciences 193 (2005) 25–49

have been developed, which give a phenomenological representation of the cell activities. Their validity for describing eucaryotic cells is, however, limited to a relatively small range of culture conditions. Nevertheless, simple models have been successfully used for process design [1], and a survey on unstructured models for hybridoma cell growth and metabolite production is given by Po¨rtner and Scha¨fer [2]. The structured kinetic model of Batt and Kompala [3] distinguishes four compartments inside the cell and considers the major substrates (glucose and amino acids) and metabolites (lactate and ammonia). Glutamine is distinguished from the other amino acids to emphasize its important role in the metabolism as a carbon and nitrogen source. In recent years, mass balancing techniques have been used to calculate the specific reaction rates in metabolic pathway networks in steady-state conditions [4,5]. In [6], the problem of underdetermined networks and the choice of additional constraints is discussed. In addition, Metabolic Flux Analysis has also been applied to determine the average network fluxes in batch and fed-batch cultures, assuming that the fluxes remain constant during the experiments [7]. More recently, Biener [8] has formulated a dynamic model based on the major metabolic pathways including amino acid metabolism. This model allows several limiting influences to be described, e.g. energy limitation, and accounts for irreversible reaction pathways. With the increasing knowledge about the cell metabolism, more detailed models are now being developed, which describe a number of metabolic features. Such modeling approaches are mostly intended for pathway analysis in metabolic engineering. On the other hand, simple macroscopic models, i.e. phenomenological models describing the dynamic evolution of a few macroscopic components (e.g. main substrates, biomass, products of interest), are now routinely used in bioprocess engineering for optimization, monitoring and control purposes [9]. As macroscopic models involve less components, and therefore less parameters, they have better identifiability properties, and offer the possibility to study the dynamics of a cell culture, which is useful for process optimization and control. Until recently, both branches of research – microscopic modeling for metabolic engineering and macroscopic modeling for process engineering – have developed nearly independently. It would however be interesting to establish a link between these two model classes. On the one hand, an existing macroscopic model can be used to check the consistency of a more sophisticated metabolic model under development. On the other hand, if a detailed metabolic model exists, it could be interesting to reduce this model, in a systematic and consistent way, to get macroscopic model candidates, and in some sense, to ÔlegitimateÕ the use of such models for the considered cell culture. This study considers these issues, and aims at providing a link between both approaches and at introducing a systematic model reduction procedure based on a few rigorous assumptions. Starting from a complex reaction network with many intracellular metabolites, an equivalent macroscopic Ôinput–outputÕ representation is derived, which involves only extracellular metabolites. This new dynamic model can be useful in situations where the full details of the intracellular metabolism are not needed, or cannot be obtained. This is particularly true in bioprocess engineering applications where dynamic models are derived from the measurements of a few extracellular metabolites and are used to design simulators, software sensors and controllers.

J.E. Haag et al. / Mathematical Biosciences 193 (2005) 25–49

27

The paper is organized as follows. In Section 2, some general ideas about mathematical modeling of industrial biochemical processes are given. In Section 3, dynamic metabolic modeling of bioprocesses is introduced. A model reduction procedure is developed in Section 4, which allows complex intracellular reaction networks to be represented by macroscopic reactions. In Section 5, a systematic search for the constraint matrix defining the macroscopic reaction scheme is presented. Section 6 summarizes and concludes the theoretical aspects and the properties of the systematic model reduction procedure.

2. Modeling of bioprocesses In biological systems of cultured biomass in suspension, at least two characteristic entities are distinguished: the culture growth medium and the biomass. The properties of the individual cells are significantly determined by the transport mechanisms across the cell membranes separating the extracellular medium from the inside of the cells, but also across the intracellular membranes separating, e.g., the cytosol from the mitochondria in eucaryotic systems, and of course by the bioreaction system, also called metabolism. Cell metabolism is characterized by a huge number of individual reaction steps regulated by modulation effects through the presence of enzymes. The complete reaction pathway is still not fully investigated with respect to stoichiometry and regulatory mechanisms, even for well studied systems like Escherichia coli and Saccharomyces cerevisiae, despite intensive research. In contrast to the high complexity of the metabolism, simple mathematical models have been developed, mostly based on experimental observations describing phenomena like limitation, activation, inhibition, saturation, multiple substrate uptake, bottlenecks and multiplicity of metabolic steady states. These models are usually valid for a limited, but often sufficiently large range of operating conditions. Of course, the level of complexity of the mathematical description depends on the application. For instance, for the purpose of controlling a system with small perturbations around a constant set-point, it is usually sufficient to use a linearized model, which is only valid in the close vicinity of the set-point. However, changing the set-point would necessitate the derivation of a new model. In other applications, for instance for process observation along a dynamic trajectory crossing a large range of operational conditions, a more general non-linear formulation of the biological system is required. Very often, the system properties (e.g. the Michaelis constants or the reaction rates) are neither known quantitatively nor even qualitatively. In these cases, a systematic and flexible modeling approach is required to find out the characteristics of the biological system. On the one hand, the model must be sufficiently general, with a sufficient number of parameters in order to represent the observed phenomena. On the other hand, the model must not be overparametrized in order to ensure identifiability from experimental data. A modeling strategy must therefore provide means to simplify the model, if it is too complex, e.g. by fixing or canceling unnecessary parameters, but also to complicate it, if the model agreement is not satisfactory. In the following section, a reduced-order reaction scheme is derived from a dynamic model containing a complex metabolic pathway. The model reduction, which is based on a few rigorous assumptions, provides a systematic working procedure.

28

J.E. Haag et al. / Mathematical Biosciences 193 (2005) 25–49

3. Dynamic metabolic modeling Consider a biological system described by a set of state variables, i.e. the concentrations of substrates, metabolic products and biomass in the extracellular solution and the intracellular intermediate metabolites. A reaction scheme with nr reactions describes the conversion of the individual components through uptake and excretion (exchange of extra- and intracellular components), intracellular metabolism and, possibly, extracellular reactions. Conversion is related by a stoichiometry, which imposes the amounts of converted mass or moles for each reaction (see for instance [8] for a complete introduction to the subject) m n X X jSi Si ! j Pi P i ; ð1Þ i¼1

i¼1

where generic notations S and P are used to denote the components consumed in the reactions and those produced in the reaction, respectively. This reaction can be rewritten in a more detailed way as 2 3 2 3 S1 P1

 6 . 7 6 . 7  T S T T T 7 7 6 6 ¼ kT n ¼ 0: ½jS1    jSm 4 .. 5 þ ½jP1    jPn 4 .. 5 ¼ kS S þ kP P ¼ kS kP ð2Þ P Sn Pn The stoichiometry of a reaction scheme with nr reactions can be written analogously 2 T3 k1 6 . 7 6 . 7n ¼ KT n ¼ 0; 4 . 5 kTnr

ð3Þ

where the array n contains all the components occurring in the reaction scheme. In an ideally mixed continuous stirred tank bioreactor (CSTR), the mass (or molar) balances for each component in the extracellular solution expressed in terms of the concentration vector ce 2 Rne and the balance for the reactor volume V 2 R are written in the form of a differential equation system: dðce V Þ ð4Þ ¼ Ke rV þ cin Qin  ce Qout ; dt dV ¼ Qin  Qout dt with terms on the right-hand side of (4) respectively accounting for

ð5Þ

• the biological reactions with the stoichiometry Ke 2 Rne nr and the reaction rate r 2 Rnr , • the reactor inflow rate Qin 2 R and a concentration cin 2 Rne , • the outflow rate Qout 2 R and the reactor concentration c 2 Rne . Eq. (4) can be rewritten as dce dV V ¼ Ke rV þ cin Qin  ce Qout  ce dt dt

ð6Þ

J.E. Haag et al. / Mathematical Biosciences 193 (2005) 25–49

29

and with Eq. (5) dce Q Q Q  Qout ð7Þ ¼ Ke r þ cin in  ce out  ce in ¼ Ke r þ ðcin  ce ÞD; dt V V V where the dilution rate is defined as D ¼ QVin . If the biomass is segregated into active (viable) and inactive (i.e. non-viable, dead) cells, their respective mass balances in (7) become: dcXv ¼ kTXv r  DcXv ¼ ðl  k d ÞcXv  DcXv ; dt

ð8Þ

dcXd ¼ kTXd r  DcXd ¼ k d cXv  DcXd ð9Þ dt assuming that there is no biomass in the inflow (cXv ;in ¼ cXd ;in ¼ 0) and dead biomass is uniquely formed from viable biomass (Xv ! Xd). The specific growth rate l 2 R is thus the sum of all specific reaction rates with the exception of the rate of mortality k d 2 R, which is considered as a morphological reaction transforming one form of biomass into another one. The total specific growth rate ltot is defined as the net rate of formation of the biomass, specific to the viable biomass: nX dcXtot X dcXi dcXv dcXd ð10Þ ¼ ¼ þ ¼ ltot cXv  DcXtot ; dt dt dt dt i¼1 which is equivalent to the viable cell growth rate, i.e. ltot = l, since there is no other biomass production pathway. Metabolic reactions in biological systems usually take place inside the biomass, and their reaction rates depend on the intracellular concentrations ci 2 Rni , which are specific to the viable cell mass (or number). A global balance of the components present in the viable biomass is written dðci cXv V Þ ð11Þ ¼ Ki qcXv V  ci k d cXv V  ci cXv Qout dt with the terms on the right-hand side representing reaction and loss due to cell death and reactor outflow. The vector q 2 Rnr contains the cell-specific reaction rates, also called metabolic fluxes in the metabolic pathway. The matrix Ki 2 Rni nr represents the corresponding stoichiometry. The application of the differentiation rule to (11) results in dci dcXv dV : cXv V ¼ Ki qcXv V  ci k d cXv V  ci cXv Qout  ci V  ci cXv dt dt dt Together with (8) and (5), the balance equation for ci is written:

ð12Þ

dci Q Q  Qout ð13Þ ¼ Ki q  ci k d  ci out  ci ðl  k d  DÞ  ci in ¼ Ki q  ci l: dt V V The balances of the extracellular and intracellular concentrations can now be formulated as the following differential equation system: c_ e ðtÞ ¼ Ke rðce ðtÞ; ci ðtÞÞ þ ge ðce ðtÞ; tÞ;

ð14Þ

c_ i ðtÞ ¼ Ki qðce ðtÞ; ci ðtÞÞ  ci ðtÞlðce ðtÞ; ci ðtÞÞ

ð15Þ

30

J.E. Haag et al. / Mathematical Biosciences 193 (2005) 25–49

with the function vector ge 2 Rne containing all the terms representing the exchange with the external world (upstream and downstream process units), in the case of a CSTR without gaseous inputs and outputs: ge ðce ðtÞ; tÞ ¼ ðcin ðtÞ  ce ðtÞÞDðtÞ:

ð16Þ

If the reaction scheme contains metabolic fluxes only, i.e. reactions catalyzed by the biomass, such as substrate uptake, product excretion, cell growth, metabolism as well as cell death, the reaction rate r is linearly related to the viable biomass concentration, rðce ; ci Þ ¼ qðce ; ci ÞcXv ;

ð17Þ

since a multiplication of catalyzing cells will lead a multiplication of the absolute metabolic reaction rates. A counterexample is given by the glutamine or protein degradation, which takes place without the presence of biomass and which therefore does not belong to the biocatalyzed or metabolic reactions. System {(14),(15)} defines completely the dynamics of a bioreaction system and takes intracellular metabolites and the corresponding metabolic pathways into account. In practice, however, the major drawback of such a modeling approach is its complexity with respect to the intracellular fluxes, and especially the dynamics, which are usually not well identified from measurements of the intracellular concentrations. It is therefore difficult to set up and use such a system representation for parameter and/or state estimation due to the lack of a priori knowledge and experimental information. Therefore, the need for a simplified representation of the reaction scheme arises. The following section deals with the systematic reduction of the number of reactions leading to a differential equation system of reduced order.

4. Model reduction The first fundamental assumption for the simplification of the model {(14),(15)} is a quasi-steady state hypothesis for the intracellular pools or at least parts of them. This assumption is verified, if the dynamics of these components are so fast compared to the remaining components that they are dominated by the slower dynamics of these latter components. This means that they are quasi statically related to the slower dynamical processes in the bioreactor. Mathematically, this is expressed by a vanishing time derivative: c_ i ðtÞ 0:

ð18Þ

This leads with (15) to the non-linear algebraic equation system Ki qðce ðtÞ; ci ðtÞÞ  ci ðtÞlðce ðtÞ; ci ðtÞÞ 0:

ð19Þ

If it is further assumed that the second term in Eq. (19) representing the dilution of the intracellular pools due to the growth of biomass is negligible with respect to the metabolic fluxes in the first term, then the algebraic equation becomes linear with respect to the flux vector q: Ki qðce ðtÞ; ci ðtÞÞ 0:

ð20Þ

J.E. Haag et al. / Mathematical Biosciences 193 (2005) 25–49

31

Eq. (20) is therefore a linear equation system with ni equations for nr fluxes. Zupke and Stephanopoulos [4] have shown that the above-mentioned assumptions are verified in most cases for mammalian cell cultures. The number of intracellular components, ni, is necessarily smaller than the number of metabolic reactions, nr. This is proved by the fact that each intracellular metabolite has at least one in- and one outgoing flux and the whole pathway is characterized by at least one flux transporting a substrate into the metabolic system inside the biomass and at least one flux leaving the biomass toward a product or the biomass itself, which is also considered as extracellular component. It is therefore impossible that the number of intracellular metabolites exceeds the number of metabolic fluxes, i.e. ni < nr. This means that the flux vector is only partly determined by the algebraic equation system (20) and still has nf = nr  ni degrees of freedom. Some additional conditions are therefore needed to fully determine the flux vector. Such conditions are usually imposed by defining one flux explicitly, e.g. by a Michealis-Menten expression, which was derived for single enzymatically catalyzed reaction steps: qj ðci Þ ¼

qmax;j ci : ci þ K M;i;j

ð21Þ

Another common condition is given by relating two or more fluxes, e.g. in a growth reaction based on essential precursors: X tj qj ¼ l ð22Þ j

with the respective yields tj. These two common approaches to formulate further conditions are linear with respect to the flux vector q and can therefore be written in the form Kf qðce ; ci Þ ¼ fðce ; ci Þ

ð23Þ

with the constraint matrix Kf 2 Rnf nr and a possibly non-linear function vector f 2 Rnf . From (19) and (23), the complete linear algebraic equation system defining the flux vector q is:



 Ki 0 fðce ; ci Þ ð24Þ qðce ; ci Þ ¼ I Kf with the (ni · nf) zero matrix 0 and the nf-dimensional identity matrix I. The flux vector itself is calculated as

1  0 Ki ð25Þ fðce ; ci Þ; qðce ; ci Þ ¼ I Kf  provided that the combined matrix KTi KTf is invertible, i.e.  T ð26Þ det KTi KTf 6¼ 0: The flux vector q defined by Eq. (25) can be substituted into the extracellular balance (14) leading, with (17), to the following new balance equation:

32

J.E. Haag et al. / Mathematical Biosciences 193 (2005) 25–49

c_ e ðtÞ ¼ Ke ½Ki Kf 1 ½0I fðce ; ci ÞcXv þ ge ðce ðtÞ; tÞ ¼ Ke

fðce ; ci ÞcXv þ ge ðce ðtÞ; tÞ Ke

ð27Þ

ne nf

with a new stoichiometry matrix 2R representing nf < nr reactions. Although the intracellular balances have been eliminated thanks to the quasi-steady state assumption and additional flux constraints, the unknown intracellular concentrations could still be present in the function vector f. In order to derive a ÔlumpedÕ network of only extracellular concentrations, it is necessary to introduce an additional assumption, i.e. the function vector f is exclusively a function of the extracellular concentrations ce: fðce ; ci Þ ¼ fðce Þ:

ð28Þ

Consequently, the original reaction scheme involving nr reactions and ni metabolites is reduced to a reaction scheme involving nf = nr  ni reactions only. These latter reactions relate extracellular components only, and allow an equivalent input–output representation of the system. It is important to stress that the potential benefit of this ÔlumpedÕ network is that extracellular metabolites can be more easily measured, and that a model of reduced complexity can be more easily identified and validated. In turn, this model can be used for bioprocess analysis, optimisation and control, all such tasks requiring a valid (but not necessarily sophisticated) dynamic model. For a given set of independent constraints f (i.e., the determinant condition (26) holds true), the solution of Eq. (25) is unique. Some alternatives might however be considered when defining the constraint set, depending on the information at hand (or the lack of information. . .). Alternatives can be checked and compared by repeatedly applying a procedure in which the model is successively reduced, identified, and validated (see also the examples in the following section). If assumption (28) does not hold, then as many additional conditions have to be defined as intracellular concentrations occurring in the function vector f. In the extreme case, ni additional constraints have to be defined in order to determine the complete vector ci, which is of course a frustrating situation in which the proposed procedure does not bring anything! In the study of complex biological systems, the lack of information about intracellular concentrations and their rate of change is certainly a limitation. However, our procedure suggests that, if additional constraints (i.e. relation between metabolic fluxes) can be expressed in terms of external concentrations only, i.e. in terms of uptake, formation, or known ratio, then significant simplifications are legitimate. The basic ideas (and assumptions) behind the reduction of a metabolic pathway into a macroscopic reaction scheme are illustrated in the following simple example. Then, more elaborate examples are developed in order to show the usefulness of the concepts and their systematic application. Example 1. Consider a biological reaction system with four components: a substrate S, a product P, the biomass X and an intracellular metabolite M. The reaction scheme consists of three reactions: q1

S ! 2M; q2

M ! P; q3

cM ! X:

ð29Þ ð30Þ ð31Þ

J.E. Haag et al. / Mathematical Biosciences 193 (2005) 25–49

33

The mortality of the biomass is neglected here for simplicity. The complete biomass is thus considered as active biocatalyst. The respective external and internal concentration vectors are defined as follows: cTe ¼ ½ s p

x ;

cTi ¼ ½m:

ð32Þ

Mass balances are given by {(14),(15)} c_ e ðtÞ ¼ Ke rðce ðtÞ; ci ðtÞÞ þ ge ðce ðtÞ; tÞ;

ð33Þ

c_ i ðtÞ ¼ Ki qðce ðtÞ; ci ðtÞÞ  ci ðtÞlðce ðtÞ; ci ðtÞÞ

ð34Þ

with the stoichiometric matrices 3 2 1 0 0 7 6 Ke ¼ 4 0 1 0 5; Ki ¼ ½ 2 0 0 1

1

c :

ð35Þ

The growth rate in (34) is calculated as l ¼ kTX q

ð36Þ

with kTX ¼ ½ 0

0 1 :

ð37Þ

Consider a kinetic structure for the fluxes, which follows a first order law with respect to the reactants of the individual reactions: q1 ¼ k 1 s;

ð38Þ

q2 ¼ k 2 m;

ð39Þ

q3 ¼ k 3 m;

ð40Þ

which can be written in the form of Eq. (23): 3 2 3 2 k1 0 0 0  k1s 7 ce 6 7 6 Kf q ¼ fðce ; ci Þ ¼ 4 k 2 m 5 ¼ 4 0 0 0 k 2 5 ; ci k3m 0 0 0 k3

ð41Þ

where Kf = I. At this stage, the constraint matrix Kf defining the kinetic structure of the reaction scheme is square (nr · nr) = (3 · 3), because no assumption has been made yet. In the following, the role of each assumption is highlighted. From the quasi-steady-state assumption (19), the intracellular concentrations in (34) are isolated as non-linear functions of the fluxes ci ðqÞ ¼ Ki q½kTX q1 ;

ð42Þ

34

J.E. Haag et al. / Mathematical Biosciences 193 (2005) 25–49

which for this example leads to 2 32 q1 6 76 m ¼ ½ 2 1 c 4 q2 54½ 0 q3

2 0

q1

331

6 77 1 4 q2 55 q3

¼

2q1  q2  cq3 : q3

ð43Þ

According to (41), the flux vector is then the solution of the following non-linear equation system: fðce ; ci ðqÞÞ  q ¼ 0:

ð44Þ

In the considered example, the solution for the fluxes and the intracellular metabolite can be computed analytically: q1 ðce Þ ¼ k 1 s;

ð45Þ

q2 ðci ðce ÞÞ ¼ k 2 mðce Þ;

ð46Þ

q3 ðci ðce ÞÞ ¼ k 3 mðce Þ; 0sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1

2

 1 k2 q ðce Þ k2 A cþ  cþ þ8 1 mðce Þ ¼ @ : 2 k3 k3 k3

ð47Þ ð48Þ

If the growth rate l is small, k3  k2, i.e. if the reaction, which leads to the formation of biomass, is slow, then the intracellular metabolite concentration cannot be calculated from (42), since the determinant of ½kTX q vanishes and the inverse cannot be built. Rewriting the algebraic equation resulting from the quasi-steady state assumption (43) as 2q1  q2  cq3  mq3 ¼ 0

ð49Þ

leads with (39) and (40) to a quadratic equation in m: k 3 m2 þ ðk 2 þ ck 3 Þm  2q1 ¼ 0;

ð50Þ

where the quadratic term vanishes if k3 0. The metabolite concentration is therefore calculated as mðq1 Þ ¼

2q1 ; k 2 þ ck 3

ð51Þ

and the respective fluxes for product formation and cell growth are q2 ðq1 Þ ¼

k2  2q1 ; k 2 þ ck 3

cq3 ðq1 Þ ¼

ck 3  2q1 : k 2 þ ck 3

ð52Þ ð53Þ

The ratio of the fluxes leaving the metabolite M is therefore determined by the respective kinetic constants k2 and k3. In this case, the ratio

J.E. Haag et al. / Mathematical Biosciences 193 (2005) 25–49

q2 k 2 ¼ ) k 3 q2  k 2 q3 ¼ 0 q3 k 3

35

ð54Þ

is invariant and could be directly used in the constraint matrix Kf in (23): Kf q ¼ fðce Þ;



 1 0 0 k1s : q¼ 0 k 3 k 2 0

ð55Þ ð56Þ

The explicit dependence of the kinetics on the intracellular metabolite concentrations ci is therefore replaced by a function of ce only. The global stoichiometry – without intracellular metabolites – is therefore defined by the matrix Ke in Eq. (27):

1  Ki 0  ð57Þ Ke ¼ Ke I Kf 3 2 1 0 3 2 32 31 2 2 1 c 0 0 1 0 0 c 7 6 2k 2 7 7 6 6 76 7 6 6 k k þ ck þ ck ð58Þ 0 5 41 05 ¼ 6 2 ¼4 0 1 0 54 1 0 3 2 37 7: 5 4 2k 3 1 0 k 3 k 2 0 1 0 0 1 k 2 þ ck 3 k 2 þ ck 3 The first column of Ke contains the stoichiometric coefficients for the reaction rate defining the uptake of the substrate S, f1 = k1s, while the second column characterizes the stoichiometry of a reaction, which does not take place, since the function f2 = 0 according to the constraint definition (56). It remains one macroscopic reaction for the extracellular components ! ! 2 kk 32 q1 2 Pþ X ð59Þ S! 1 þ c kk32 1 þ c kk32 with the reaction rate q1 in (38). Thus, the transformation (24) maps the metabolic pathway network into a reduced-order global reaction scheme, which gives a black-box Ôinput–outputÕ representation of the biomass. The systematic application of the procedure is now demonstrated in the next examples. Example 2. Consider a simple metabolic reaction scheme with nr = 4 reactions involving two substrates (S), one product (P), ni = 2 metabolites (M) and the biomass (X): (1) (2) (3) (4)

S1 ! 2M1, S2 ! M2, M1 ! 2M2 + P, M2 ! 0.1X.

For a graphical representation of the metabolic reaction scheme (without stoichiometry), see Fig. 1.

36

J.E. Haag et al. / Mathematical Biosciences 193 (2005) 25–49

S1

M1

S2

M2

P

X

Fig. 1. Metabolic reaction scheme for Example 2.

With the state vectors cTe ¼ ½ s1 s2 p x  and cTi ¼ ½ m1 matrices are defined as 2 3 1 0 0 0

 6 0 1 0 0 7 2 0 1 0 6 7 : Ke ¼ 6 7; Ki ¼ 4 0 0 1 0 5 0 1 2 1 0 0 0 0:1

m2 , the respective stoichiometric

ð60Þ

Since nr  ni = 2, two additional equations are required to complete the definition of the reaction system, which are chosen as follows in this example: q1 ¼ f1 ðcÞ;

ð61Þ

q2 ¼ f2 ðcÞ:

ð62Þ

The constraint matrix is then written:

 1 0 0 0 : Kf ¼ 0 1 0 0 From Eq. (27), the following global stoichiometric matrix results: 3 2 1 0 6 0 1 7 7 6 Ke ¼ 6 7; 4 2 0 5 0:4

0:1

which is equivalent to the global reaction scheme (1) S1 ! 2P + 0.4X, (2) S2 ! 0.1X (see Fig. 2).

ð63Þ

ð64Þ

J.E. Haag et al. / Mathematical Biosciences 193 (2005) 25–49

S1

37

P

S2

X

Fig. 2. Simplified reaction scheme for Example 2.

It is to be noted that no a priori knowledge on the kinetic structure is needed to simplify the metabolic reaction scheme. The only restriction is its exclusive dependence on extracellular concentrations. The influence of the different types of constraints is illustrated in different cases of the following, more complex example: Example 3. Consider the biological system consisting of two substrates, S1 and S2, three metabolic products, P1, P2 and P3, and two intracellular intermediates, M1 and M2. The following reaction scheme is assumed (Fig. 3): (1) (2) (3) (4)

q1

S1 ! 2M1 þ 2M2 þ 2P3 , q2 M1 þ M2 ! P1 , q3 M1 ! 4M2 þ 3P2 þ P3 , q4 S2 þ 2M2 ! 6P3 .

S1 is directly converted into P3 via reaction 1 and indirectly via reactions 3 and 4. P2 is indirectly produced from S1 via reactions 1 and 3, and P1 is indirectly produced from S1 via reactions 1 and 2. The only product of S2 is P3. The concentration vectors are defined by cTe ¼ ½ s1

s2

p1

p2

p3 ;

cTe ¼ ½ m1

m2 ;

ð65Þ

which leads to the following stoichiometric matrices:

S1

M1

P1 P2

S2

M2

P3

Fig. 3. Simple reaction scheme with four reactions. Stoichiometry is not shown.

38

J.E. Haag et al. / Mathematical Biosciences 193 (2005) 25–49

2

1

6 6 0 6 Ke ¼ 6 6 0 6 4 0 2

0

0

0

0

1

0

0

3

0

1

0

3

7 1 7 7 0 7 7; 7 0 5 6

Ki ¼

2

1

1

0

2

1

4

2

 ð66Þ

:

Since nr = 4 reactions are formulated with ni = 2 balancing metabolites, two additional algebraic constraints have to be introduced. Three cases are considered in the sequel: Case 1 (Constraint on substrate consumption). If the substrate uptake rates are given by the following expressions: qS1 ¼ q1 ¼ f1 ðcS1 Þ;

ð67Þ

qS2 ¼ q4 ¼ f2 ðcS2 Þ;

ð68Þ

then the constraint arrays are written:



 1 0 0 0 f1 Kf ¼ ; f¼ : 0 0 0 1 f2 Consequently, the flux vector q can be written as a linear transformation of the specified reaction rate vector f(s1, s2) according to Eq. (25): 2 3 2 3 q1 1 0 6 q 7 6 2 0:4 7 f ðs Þ  6 27 6 7 1 1 ; ð69Þ 6 7¼6 7 4 q3 5 4 0 0:4 5 f2 ðs2 Þ q4

0

1

and the following global reaction scheme results according to the matrix Ke in Eq. (27): f1

(1) S1 ! 2P1 þ 2P3 , f2 (2) S2 þ 0:4P1 ! 1:2P2 þ 6:4P3 , which is visualized in Fig. 4.

S1

P1 P2

S2

P3

Fig. 4. Global reaction scheme for first case consisting of two reactions (stoichiometry not shown).

J.E. Haag et al. / Mathematical Biosciences 193 (2005) 25–49

39

The explanation of the first reaction with rate f1 is straightforward, as it can easily be obtained by linear combination of reactions 1 and 2, and the yield from substrate S1 to product P3 remains the same as in reaction 1. The second reaction is the consequence of the need for the intermediate metabolite M2 in reaction 4, provided by reactions 2 and 3, which themselves produce another 0.4 P3 per S2 consumed. The net yield from S2 to P3 increases therefore compared to the yield of the single reaction 4, and the metabolite is equivalently replaced by the other products involved in the intermediate reactions 2 and 3. Thus, with a model reduction from four to two reactions, the Ôcomplex systemÕ of metabolic pathways is lumped in a global reaction scheme for the extracellular components Technically, the Ôwhite-boxÕ model is reformulated as an equivalent Ôinput–outputÕ representation. Case 2 (Constraint on product formation). The effect of a small modification in the constraint formulation is investigated next: the second constraint (68) on the uptake rate of S2 is replaced by the specification of the production rate of P3, qP3 ¼ 2q1 þ q3 þ 6q4 ¼ f2 ðce Þ:

ð70Þ

In this case, the constraint matrix becomes

 1 0 0 0 Kf ¼ ; 2 0 1 6 and the resulting global stoichiometry is written as f1

(1) S1 þ 38 P2 ! (2)

5 S 32 2

5 S 16 2 f2 3

þ 161 P1 !

16

þ 2 18 P1 , P2 þ P3 ,

which is schematically shown in Fig. 5. The production of P3 is here, in contrast to the previous case, decoupled from the uptake of S1. This is due to the explicit definition of the kinetics of these two components. The substrate S2 and the product P1 are treated as intermediates in the pathway from S1 to P3, while P2 plays the role of a catalyst in this reaction chain. The combination of the two resulting reactions by eliminating the ÔintermediateÕ S2 reveals that for each molecule of S1 consumed in the first reaction, two molecules of P3 are produced with two molecules of P1 as a side product, which corresponds exactly to the first reaction of Case 1. However, the kinetics of this overall reaction is different.

S1

P1 P2

S2

P3

Fig. 5. Global reaction scheme for second case consisting of two reactions (stoichiometry not shown).

40

J.E. Haag et al. / Mathematical Biosciences 193 (2005) 25–49

Case 3 (Linearly dependent fluxes). Another kind of constraint that is often used to define a metabolic system kinetically is to linearly relate – i.e. to associate – fluxes to each other. In this example, such a constraint could be to fix the ratio between the consumption of S1 and the production of P1, q4 ¼ a ) aq2  q4 ¼ f2 ¼ 0: ð71Þ q2 Using (71) with a = 0.1 as the second constraint, while the first remains the same as in the first case, the global reaction scheme becomes f1

(1) S1 þ 0:19S2 ! 1:92P1 þ 0:23P2 þ 3:23P3 , (2) 1.15P2 + 6.15P3 9 0.96S2 + 0.38P1, see also Fig. 6. Since the right-hand side of the second constraint (71) is equal to zero, the global reaction system is reduced to one equation only, and the second equation in the system becomes meaningless. In this case, the complete set of reactions 1–4 is reformulated in one global reaction driven by the uptake rate of substrate S1. This global reaction scheme represents a stronger model reduction compared to the previous cases, since there is only one global transformation from the substrates to the products, i.e. the dynamics of all component concentrations are linearly dependent. The missing degree of freedom is taken by fixing the value of a, which could also be considered as an unknown parameter in the model identification procedure. One problem becomes apparent in this example, which was not evident in the detailed reaction system: The (global) kinetics is formulated as a function of the concentration of S1 and does not depend on the second substrate, which could however be limiting under certain circumstances, since it appears as a reactant in the global reaction equation. Clearly, there is no precaution that prevents the concentration cS2 from becoming negative in the case where s2 = 0 ^ s1 > 0. To avoid this problem, the kinetics would have to be rewritten such that f1 ðsi ¼ 0Þ  0;

i ¼ 1; 2:

ð72Þ

Otherwise the model would be restricted to cases with non-limiting concentrations of S2, which would limit the range of validity of the reaction model.

S1

P1 P2

S2

P3

Fig. 6. Global reaction for third case (stoichiometry not shown).

J.E. Haag et al. / Mathematical Biosciences 193 (2005) 25–49

41

Example 4. As a last illustrative example, the complex metabolic network of hybridoma cells by Zupke and Stephanopoulos [4] is considered. It consists of nr = 20 reactions including glycolysis (EMP), the pentose phosphate pathway (PPP), the tricarboxylic acid (TCA) cycle, the oxidative phosphorylation as well as cell maintenance, cell growth and antibody production. 23 components are balanced, and among them ni = 14 intracellular metabolites and co-metabolites are assumed in a quasi-steady state (see Table 1 and Fig. 7). Hence, nr  ni = 6 additional constraints are necessary to fully define the specific reaction rate vector q. The following assumptions are therefore introduced: (1) The glucose uptake depends on its extracellular concentration, e.g. Michaelis-Menten kinetics: qGlc ¼ q1 ¼ f1 ðcGlc Þ: (2) The glutamine uptake is also described in function of its extracellular concentration: qGln ¼ q13 þ 0:036q19 þ 0:012q20 ¼ f2 ðcGln Þ: (3) The specific maintenance rate is constant: q18 ¼ f3 ¼ mATP : Table 1 Reaction scheme of hybridoma metabolism considered by Zupke and Stephanopoulos [4] q1 q2 q3 q4 q5 q6 q7 q8 q9 q10 q11 q12 q13 q14 q15 q16 q17 q18 q19 q20

Glc + ATP ! G6P G6P ! F6P F6P + ATP ! 2GAP GAP ! Pyr + NADH + 2ATP Pyr + NADH ! Lac G6P ! Ru5P + 2 NADPH + CO2 Ru5P ! X5P Ru5P + X5P ! S7P + GAP S7P + GAP ! F6P + E4P X5P + E4P ! F6P + GAP Pyr + Mal ! aKG + 3 NADH + CO2 aKG ! 2 NADH + CO2 + ATP + Mal Gln ! Glu + NH3 Pyr + Glu ! Ala + aKG Glu ! aKG + NADH + NH3 2 NADH + O2 ! 6ATP Mal ! Pyr + CO2 + NADPH ATP ! 0.036Gln + 0.062Glu + 0.0031G6P + 0.013Ru5P + 0.042GAP + 0.123Pyr + 0.019Mal + 0.74ATP + 0.188NADPH ! X + 0.073aKG + 0.23NADH + 0.12CO2 0.012Gln + 0.07Glu + 0.04GAP + 0.011NADPH + 0.012Mal ! Ab + 0.082 NADH + 0.057aKG

All bioreactions are catalyzed by the viable cells, assumed to be homogeneous in their properties. Reaction rates are specific to the viable cell density. The stoichiometric coefficients in reactions 19 and 20 are in pmol/cell and pmol, respectively. Extracellular components are printed in bold. Quasi-steady state assumption is applied to the intracellular components.

42

J.E. Haag et al. / Mathematical Biosciences 193 (2005) 25–49 Glc 1

ATP

Ab

NADPH CO2

G6P

X

Gln

X

6

Ru5P 7

2

X5P 10

ATP

9

3

X Ab

X Ab

S7P

GAP

NH3

4

NADH

Lac

13

E4P

F6P

5

ATP

Pyr

CO2 NADPH

X Ab

3 NADH 15

17

CO2

Mal

14

α KG

CO2 2 NADH

X Ab

Glu 11

NH3

Ala X Ab

12

2 NADH ATP

18

(ADP)

X Ab

NADH

16

6 ATP

O2

Fig. 7. Major metabolic reaction pathway of hybridoma cells.

(4) The antibody production follows Luedeking–Piret kinetics [10]: qAb ¼ q19 ¼ a þ bl ()

 bq19 þ q20 ¼ f4 ¼ a:

(5) The reaction rate ratio c of the two alternative pathways from glutamine to a-ketoglutarate is constant: q15 ) cq14 þ ðc  1Þq15 : c¼ q14 þ q15 (6) The specific growth rate depends on the concentrations of glucose, glutamine, lactate and ammonia as modulators: l ¼ q19 ¼ f6 ðcGlc ; cGln ; cLac ; cNH3 Þ: Respecting these assumptions, the detailed modeling of the most important metabolic pathways is simplified to a global reaction scheme with five reactions only, since six additional constraints with five non-zero elements in f are introduced. The resulting reaction scheme is shown in Table 2. Table 2 Simplified global reaction scheme of hybridoma metabolism f1 f2 f3 f4 f6

Glc + 0.21CO2 ! 2.11Lac + 0.32O2 Gln ! (0.54 + c)Lac + (1c)Ala + (1 + c)NH3 + 0.13O2 + 0.75CO2 0.053Lac + 0.158O2 ! 0.105CO2 (0.082 + 0.011c)Lac + 0.082 (1c)Ala + (0.012 + 0.082c)NH3 + 0.018 O2 + 0.022 CO2 ! Ab (0.187 + 0.011b + 0.098c + 0.082bc)Lac + (1c) (0.098 + 0.082b)Ala + (0.154 + 0.018b)O2 + (0.036 + 0.012b + 0.098c + 0.082bc)NH3 ! X + bAb + (0.186  0.022b)CO2

J.E. Haag et al. / Mathematical Biosciences 193 (2005) 25–49

43

Flux f1 characterizes the uptake of glucose and its way through the glycolysis with the formation of the waste metabolite lactate. Flux f2 represents the glutamine uptake and its metabolism via glutamic acid, where alanine and ammonia are produced as side products. As a strange side effect, oxygen appears as a product in both reactions, which is physically impossible, since mammalian cells do certainly not produce and excrete oxygen. This result originates from the definition of the constraints, where the glucose and glutamine consumption as well as the growth rates are fixed by kinetic expressions. On the other hand, freedom is left to the consumption rate of oxygen, which occurs in flux q16, the oxidative phosphorylation. It is therefore closely related to the formation of energy in the cell and might be seen as an Ôextracellular energy equivalentÕ (although this makes biologically no sense). For instance, maintenance (f3) is globally represented by the consumption of oxygen, whereas, in reality, there is a consumption of the intracellular energy equivalent ATP in the metabolic pathway. The global reactions producing oxygen might therefore be seen as reducing the total oxygen consumption. This makes sense, as long as the total oxygen uptake rate OUR ¼ qO2 ¼ 0:32f 1  0:13f 2 þ 0:158f 3 þ 0:018f 4 þ ð0:154 þ 0:018bÞf6

ð73Þ

remains positive. Care has therefore to be taken in the formulation of the kinetic expressions fi(ce). In a similar way, the consumption of lactate with the flux f3 can be explained as a reduction of the overall lactate production. Also in this case, the kinetics have to be defined appropriately in order to prevent the unrealistic case of lactate consumption. The antibody is globally produced in two independent steps due to the assumption of partly growth-associated production kinetics. It occurs therefore independently of growth in the global reaction f4 and as a side product of growth f6. The biosynthesis globally leads to the consumption of all essential precursors. Glucose and glutamine do not appear in the growth reaction, because their uptake rates are formulated independently of the growth rate. Metabolic stoichiometry is the subject of intense research led by biologists, and can be considered as relatively well known a priori. In contrast, there is usually little information about the choice of the constraint matrix Kf. However, the preceding examples have demonstrated the key influence of this latter matrix on the resulting global macroscopic reaction schemes. Efforts have therefore to be put in the selection of appropriate constraints on the metabolic fluxes, and the next section deals with a systematic procedure for the definition of a constraint matrix. 5. Systematic construction of the constraint matrix Kf Since the metabolic stoichiometry is supposed to be known, the search for a global, macroscopic reaction scheme reduces to the search for an appropriate constraint matrix Kf. The full matrix Kf contains nf · nr = (nr  ni) · nr elements. With a known function fi in the ith constraint kTf;i q ¼ fi ðce Þ;

ð74Þ

the nr elements in kf,i would be linearly independent. However, in the case of a yet unknown function fi (which is the case in practice), the corresponding row of the constraint matrix could be normalized with respect to one element in kf,i or to one factor in fi.

44

J.E. Haag et al. / Mathematical Biosciences 193 (2005) 25–49

A simple example is given by the direct definition of a single reaction rate as constraint, e.g. by the classical Michaelis-Menten function c ; ð75Þ q ¼ f ðcÞ ¼ qmax c þ KM which could be written in the following two forms: c 1  q ¼ qmax ; ð76Þ c þ KM c : ð77Þ c þ KM In (76), the constraint is normalized with respect to (in this case sole) one element in the constraint matrix, while in (77), the constraint is normalized with respect to the maximum reaction rate qmax. Hence, a full matrix Kf with a priori unknown kinetics leads to the problem of estimating nf(nr  1) linearly independent unknown parameters directly related to the stoichiometry. In order to avoid additional stoichiometric parameters, the simplest form of constraint matrix is considered in the following, which contains a single non-zero element in each row. This element, which is normalized to one, corresponds therefore to the definition of a specific reaction rate. In this case, the maximum number of possible global reaction schemes is the number of permutations

 nr : ð78Þ N¼ nf q1 max  q ¼

In order to evaluate the quality of each individual choice for Kf, the kinetic functions f have to be structurally defined and parametrized. This definition has to be done in an as general as possible way in order to avoid shortcuts in the evaluation. In this study, the following general structure is proposed: ne Y aij ðcij Þ; i ¼ 1 . . . nf ð79Þ fi ðce Þ ¼ fmax;i j¼1

with the multiplicative modular functions 8 cj sgnfK ij g if K > 0 < cj þK 2 cj ij aij ðcij Þ ¼ sgnfK ij g ¼ 1 otherwise 2 c þ K2 : j

ij

ð80Þ

1þcj K ij

featuring two fundamental modulation effects: • limitation according to the Michaelis-Menten law for Kij > 0: fi is monotonically increasing with cj; the Michaelis constant is K M;ij ¼ K 2ij ; • non-competitive inhibition for Kij < 0: fi is monotonically decreasing with cj; the inhibition constant is calculated as K I;ij ¼ K 2 ij . At Kij = 0, the influence of cj vanishes, i.e. component j has no impact on the reaction i. Note that the unit of Kij changes in a reciprocal way at this point. Pure first order kinetics are a special case f . of Kij > 0 with Kij ! 1 and a proportionality factor k i ¼ max;i K ij

J.E. Haag et al. / Mathematical Biosciences 193 (2005) 25–49

45

Hence, the unknown model parameters are • nf maximum reaction rates qmax,i and • nfnc modulation constants Kij, i.e. the model parameter vector p contains nf(nc + 1) elements. The structural parameters have to be estimated for each permutation k, with 1 6 k 6 N, and the cost functions j(p*) at the optimum have to be evaluated and compared. The best choice for the constraint matrix is then the permutation corresponding to the lowest cost: k  ¼ arg minfjðKf;k ; pk Þg

ð81Þ

pk ¼ arg minfjðKf;k ; pk Þg:

ð82Þ

k

with p

Due to the structure of the kinetic functions fi, the parameter estimation problem is smooth and structurally identifiable. The quality of the estimates depends of course also on the experimental conditions chosen to identify all unknown parameters. It is intuitively clear that all concentrations must undergo a dynamic variation in order to allow a good estimation of their corresponding modulation constants, i.e. steady state experiments would be inappropriate. Example 5. To illustrate the effectiveness of the constraint search procedure, consider again Example 2 in Section 4. The reference, i.e. the ÔtrueÕ system is characterized by the following kinetics: cS1 ; ð83Þ f1 ¼ q1;max cS1 þ K M1 f2 ¼ q2;max

cS2 KI cS2 þ K M2 cP þ K I

ð84Þ

with the parameter values given in Table 3. Note that Eqs. (83) and (84) match the general kinetics structure (79). The bioreactor is operated in perfusion mode, i.e. the biomass is retained in the bioreactor. The function ge(t) in the extracellular balance (14) then takes the following form: ge ðcðtÞ; tÞ ¼ ½cin  cout ðcðtÞÞDðtÞ:

ð85Þ

Table 3 Parameter values for ÔtrueÕ reference system Parameter

Value

q1,max q2,max KM1

1 = f1,max 2 = f2,max 1 = K 211

KM2

2 = K 222

KI

5 = K 2 23

46

J.E. Haag et al. / Mathematical Biosciences 193 (2005) 25–49

Table 4 Experimental conditions for model identification Experimental parameter

Value

Experiment duration Sampling interval Perfusion rate

tend = 30 Dt = 1  2 2t DðtÞ ¼ 1  tend 1

Initial state

cT0 ¼ ½ 1 10 0:1 0:1 

Feed concentrations Relative measurement error

cTin ¼ ½ 10 1 0 rrel = 1%

0

The objective is to assess the correct constraint matrix Kf defined in (63) among the N = 6 candidate schemes from the experiment with the reference scheme. To this end, the reference system is simulated numerically with the initial value solver routine ode15s in MATLAB with outputs every time step Dt. The experimental parameters and functions are listed in Table 4. Measurement noise with a relative standard deviation rrel is added to the results. This data set is used to compute a maximum-likelihood estimate of the unknown model parameters. Note that the initial conditions cT0 of the differential equation system (14) are also unknown because of the measurement error and are therefore subject to estimation. Thus, 14 unknowns are determined from 124 experimental data values, i.e. the redundancy factor of the estimation problem is 8.9. The initial estimates for the parameter estimation are chosen as follows: • ^c0 ¼ y0 , • ^k i ¼ 1,  b ij ¼ 1 • K 0:1

if fKgij < 0; otherwise:

The evaluation of the minimum cost for all candidate schemes leads to the ranking given in Table 5. Moreover, Table 6 computes the estimated kinetic parameter values corresponding to the best scheme with those of the reference model. From these tables and Fig. 8, it is apparent that the systematic search procedure provides good estimates of the exact pseudo-stoichiometry and kinetics. An important underlying assumption in this illustrative example is the knowledge of the kinetic structure (79). A careful analysis of the kinetic structure is required in real case studies. Table 5 Ranking of candidate schemes rank

Flux combination

ÔtrueÕ

{1, 2}

1 2 3 4 5 6

{1, 2} {2, 3} {2, 4} {3, 4} {1, 4} {1, 3}

Cost 259.7 264.1 578.6 1276.4 83865.4 126889.5  T Ki KTf singular

J.E. Haag et al. / Mathematical Biosciences 193 (2005) 25–49

47

Table 6 Kinetic parameter values for the ÔbestÕ reaction scheme (flux combination {1,2}) Parameter

Identified value

qmax,1 qmax,2 K11 K12 K13 K14 K21 K22 K23 K24

ÔtrueÕ value

1.00 1.92

1.00 2.00

0.96 0.05 0.04 0.03 0.00 1.36 0.45 0.00

1.00 0.00 0.00 0.00 0.00 1.41 0.45 0.00

c

S

1

10 5 0

0

10

20

30

0

10

20

30

0

10

20

30

0

10

20

30

10

c

S2

15

5 0

c

P

20 10 0

cX

100 50 0

t

Fig. 8. ÔExperimentalÕ data (+) from reference model in comparison with identified general model ().

6. Discussion Under the conditions that (i) the biomass growth is slow and (ii) the intracellular kinetics are so fast that the dynamics of the intracellular concentrations can be neglected compared to

48

J.E. Haag et al. / Mathematical Biosciences 193 (2005) 25–49

extracellular dynamics, the metabolic reaction system can be reduced to an unstructured representation of the cell behavior. The dimension of the reduced reaction system depends on the degree of freedom left to the system and is equal to the difference between the dimension of the metabolic reaction system and the number of balancing metabolites assumed in quasi-steady state. Zupke and Stephanopoulos [4] have shown that the above-mentioned assumptions are fulfilled with good approximation for hybridoma cell lines, which are characterized by rather complex metabolic pathways. As complex metabolic reaction schemes can be reduced to equivalent macroscopic reaction schemes relating extracellular components only, the use of unstructured models is henceforth legitimate. However, depending on the introduction of specific kinetic constraints, different macroscopic reaction schemes can be generated, among which some of them may be biologically acceptable, some others may not. In an attempt to overcome this potential difficulty, a systematic procedure for evaluating these different schemes is introduced. At this stage, a general formulation of the specific reaction rates featuring modulation effects (from limitation to inhibition) is proposed. This formulation uses only one parameter per component and reaction, and one maximum specific reaction rate per reaction. This keeps the number of structural model parameters low and offers good identifiability properties. The potential benefit of the resulting ÔlumpedÕ network, which appears in the form of an Ôinput– outputÕ representation, is that extracellular metabolites can be more easily measured, and that a model of reduced complexity can be more easily identified and validated. In turn, this model can be used for bioprocess analysis, optimisation and control, all such tasks requiring a valid (but not necessarily sophisticated) dynamic model.

Acknowledgment The financial support of the company 4C Biotech s.a., Seneffe, Belgium is gratefully acknowledged.

References [1] S. Dhir, J. Morrow Jr., R.R. Rhinehart, T. Wiesner, Dynamic optimization of hybridoma growth in a fed-batch bioreactor, Biotechnol. Bioeng. 67 (2) (2000) 197. [2] R. Po¨rtner, T. Scha¨fer, Modelling hybridoma cell growth and metabolism – a comparison of selected models and data, J. Biotechnol. 49 (1996) 119. [3] B.C. Batt, D.S. Kompala, A structured kinetic modeling framework for the dynamics of hybridoma growth and monoclonal antibody production in continuous suspension cultures, Biotechnol. Bioeng. 34 (1989) 515. [4] C. Zupke, G. Stephanopoulos, Intracellular flux analysis in hybridomas using mass balances and in vitro 13C NMR, Biotechnol. Bioeng. 45 (4) (1995) 292. [5] H.P.J. Bonarius, V. Hatzimanikatis, K.P.H. Meesters, C.D. de Gooijer, G. Schmid, J. Tramper, Metabolic flux analysis of hybridoma cells in different culture media using mass balances, Biotechnol. Bioeng. 50 (3) (1996) 299. [6] H.P.J. Bonarius, G. Schmid, J. Tramper, Flux analysis of underdetermined metabolic networks: the quest for the missing constraints, TIBTECH 15 (1997) 308. [7] L. Xie, D.I.C. Wang, Material balance studies on animal cell metabolism using a stoichiometrically based reaction network, Biotechnol. Bioeng. 52 (5) (1996) 579.

J.E. Haag et al. / Mathematical Biosciences 193 (2005) 25–49

49

[8] R. Biener, Modellbildung und modellgestiitzte ProzeBfiihrung bei tierischen Zellkulturen, Fortschritt-Berichte VDI, Reihe 17, VDI Verlag, Dusseldorf, vol. 178, 1998. [9] G. Bastin, D. Dochain, On-Line Estimation and Adaptive Control of BioreactorsProcess Measurement and Control, vol. 1, Elsevier, Amsterdam, 1990. [10] R. Luedeking, E.L. Piret, A kinetic study of the lactic acid fermentation. Batch process at controlled pH, J. Biochem. Microbiol. Technol. Eng. 1 (4) (1959) 393.