Computer-implemented method for creating a fermentation model

Abstract

The application relates to a computer-implemented method for creating a model of a bioreactionfermentation process or whole-cell catalysis processwith an organism on the basis of measured data.

Claims

1. A computer-implemented method for process control of a bioreaction with an organism comprising: creating a model of a bioreaction with an organism, said model including a matrix L, kinetics of macroreactions and model parameters of calculated values, which comprises the following steps: a. defining selected metabolic pathways of the organism, their properties of stoichiometry and reversibility as background knowledge and calculating elementary modes from this input; b. combining the elementary modes from a) in a matrix K, wherein the elementary modes combine the metabolic pathways from a) in macroreactions and the matrix K contains the stoichiometry and reversibility properties of all macroreactions; c. entering measured data for the bioreaction with the organism; d. calculating, using an interpolation method, specific rates for the organism of the metabolic pathways based on the measured data entered from c); e. selecting relevant macroreactions as a subset r(t) of the elementary modes from a) by i. data-independent and/or data-dependent prereduction of the elementary modes from a); ii. selection of the subset from the prereduction from e) i. using the measured data from c) and/or one or more rates from d) by means of an algorithm according to a mathematical quality criterion and combination of the subset in a matrix L; iii. optionally, the subset is shown graphically, f. calculating, using an interpolation method, reaction rates of the macroreactions of the subset r(t) on the basis of the measured data from c) and/or the rates from d); g. devising kinetics of the macroreactions of the subset from e) ii. using the following intermediate steps: i. devising generic kinetics from the stoichiometry of the macroreactions from e); ii. determining factors influencing the macroreactions from e) from the reaction rates from f); iii. expanding the generic kinetics from g) i. by model parameter values which quantify the factors determined in g) ii; h. performing separately for each macroreaction, a first adjustment of the model parameter values of the model parameters of g) iii to the calculated reaction rates from f); i. optionally repeating steps g) and h) until a predefined quality of adjustment is reached; j. adjusting the model parameter values of the model parameters of g) iii, h) or i) to the measured data from c); k. forming the matrix L from e) ii, the kinetics from g) iii and the model parameter values of the model parameters from j) as an output for transferring to a process control module or process development module; and l. closed-loop controlling the bioreaction using the output of step k).

2. The computer-implemented method according to claim 1, wherein in step d), growth rates of the organism, and optionally also death rates of the organism, are calculated.

3. The computer-implemented method according to claim 1, wherein in step g), an individual adjustment of the kinetics takes place, based on an analysis of the reaction rates from f).

4. The computer-implemented method according to claim 1, wherein in step h), the adjustment of the parameter values of the kinetics from g) takes place by combining several methods of adjustment.

5. The computer-implemented method according to claim 1, wherein in step e) ii., for selecting the subset of the macroreactions, linear estimation of reaction rates of selected macroreactions is carried out.

6. The computer-implemented method according to claim 1, wherein in step e) ii., for selecting the subset of macroreactions, linear estimation of reaction rates of selected macroreactions is carried out in combination with an evolutionary algorithm.

7. The computer-implemented method according to claim 1, wherein the measured data are shifted before using the interpolation method in step d), in order to achieve the description of constant consumption without feed peaks.

8. The computer-implemented method according to claim 1, wherein in step f), linear estimation of reaction rates of selected macroreactions is carried out.

9. The computer-implemented method according to claim 1, wherein in step e) i., a data-dependent prereduction is carried out and the method of linear estimation of reaction rates of selected macroreactions with NNLS is used for this.

10. The computer-implemented method according to claim 1, wherein in step e) iii., validity of selection of the subset of macroreactions is tested by means of a flux map.

11. The computer-implemented method according to claim 1, wherein in step e) ii., the selection from the prereduction from e) i. is carried out using the measured data from c).

12. The computer-implemented method according to claim 1, wherein the method is executed by a computer program.

13. The computer-implemented method according to claim 1 wherein the method is executed by a software program.

Description

EXAMPLE

Modeling of a Hybridoma Cell Culture

(1) 1 Step a)

(2) The background knowledge in the form of a metabolic network was taken from the work of Niu et al. (Metabolic pathway analysis and reduction for mammalian cell culturesTowards macroscopic modeling. Chemical Engineering Science (2013) 102, pp. 461-473. DOI: 10.1016/j.ces.2013.07.034.). The metabolic network of an animal cell described here contains 35 reactions, which link together 37 internal and external metabolites (see FIG. 5; see Table 1).

(3) TABLE-US-00001 TABLE 1 Reactions of the metabolic network accordingto Niu et al. (Metabolic pathway analysis and reduction for mammalian cell cultures-Towards macroscopic modeling. Chemical Engineering Science (2013) 102, pp. 461-473. DOI:10.1016/j.ces.2013.07.034.) 1 Glucose .fwdarw. 1 G6P 1 G6P + 2 NAD .fwdarw. 2 Pyruvate 1 Pyruvate .fwdarw. 1 Lactate + 1 NAD 1 Pyruvate .fwdarw. 1 Pyruvate_m 1 NADm + 1 Pyruvate_m .fwdarw. 1 Acetyl coA_m 1 Acetyl coA_m + 1 NADm + 1 Oxaloacetate_m .fwdarw. 1 -ketoglutarate_m 1-ketoglutarate_m + 1 NADm .fwdarw. 1 Succinyl CoA_m 1 FADm + 1 Succinyl CoA_m .fwdarw. 1 Fumarate 1 Fumarate .fwdarw. 1 Malate_m 1 Malate_m + 1 NADm .fwdarw. 1 Oxaloacetate_m 1 Glutamine .fwdarw. 1 Glutamate + 1 NH3 1 Glutamate + 1 NADm .fwdarw. 1 -ketoglutarate_m + 1 NH3 1 Malate_m .fwdarw. 1 Malate 1 Malate + 1 NAD .fwdarw. 1 Pyruvate 1 Glutamate + 1 Pyruvate .fwdarw. 1 -ketoglutarate_m + 1 Alanine 1 Glutamate + 1 Oxaloacetate_m .fwdarw. 1 -ketoglutarate_m + 1 Aspartate 1 Arginine + 2 NADm .fwdarw. 1 Glutamate + 3 NH3 1 Asparagine .fwdarw. 1 Aspartate + 1 NH3 2 Glycine + 1 NADm .fwdarw. 1 NH3 1 Histidine + 1 NADm .fwdarw.1 Glutamate + 2 NH3 1 Isoleucine + 2 NADm .fwdarw. 1 Acetyl coA_m + 1 NH3 + 1 Succinyl CoA_m 1 Leucine + 3 NADm .fwdarw. 3 Acetyl coA_m 1 Lysine + 6 NADmv 2 Acetyl coA_m 1 Methionine + 4 NADm .fwdarw. 1 NH3 + 1 Succinyl CoA_m 1 NADm + 1 Phenalanine .fwdarw. 1 Tyrosine 1 Serine .fwdarw. 1 NH3 + 1 Pyruvate 1 NADm + 1 Threonine .fwdarw.1 NH3 + 1 Succinyl CoA_m 19 NADm + 1 TRP .fwdarw. 3 Acetyl coA_m 5 NADm + 1 Tyrosine .fwdarw.2 Acetyl coA_m + 1 Fumarate 5 NADm + 1 Valine .fwdarw. 1 NH3 + 1 Succinyl CoA_m 1 NADm .fwdarw. 1 NAD 0.5 Oxygen(O2) .fwdarw. 1 NADm 1 NADm .fwdarw. 1 FADm 0.0156 Alanine + 0.0082 Arginine + 0.0287 Aspartate + 0.0167 G6P + 0.0245 Glutamine + 0.0039 Glutamate + 0.0196 Glycine + 0.0038 Histidine + 0.0099 Isoleucine + 0.0156 Leucine + 0.0119 Lysine + 0.0039 Methionine + 0.0065 Phenylalanine + 0.016 Serine + 0.0094 Threonine + 0.0047 Tyrosine + 0.0113 Valine .fwdarw.1 X (Biomass) + 0.0981 NAD 0.01101 Alanine + 0.005033 Arginine + 0.007235 Asparagine + 0.0081787 Aspartate + 0.010381 Glutamine + 0.010695 Glutamate + 0.01447 Glycine + 0.0034602 Histidine + 0.005033 Isoleucine + 0.014155 Leucine + 0.01447 Lysine + 0.0028311 Methionine + 0.007235 Phenylalanine + 0.026738 Serine + 0.016043 Threonine + 0.0084932 Tyrosine + 0.018874 Valine .fwdarw.1 IgG (Antibody)

(4) Reversibility of the reactions is not explicitly stated in the published work. Instead, the data on metabolic flux analysis from the same publication were evaluated and were used for identifying the irreversible reactions.

(5) With the stoichiometric matrix N, which contains the stoichiometry, i.e. the stoichiometric coefficients, of the internal metabolites, and the information about the reversibility of the reactions, all elementary modes (EMs) of the network were calculated using METATOOL 5.1 (Pfeiffer et al. METATOOL: for studying metabolic networks, Bioinformatics 199915 (3), pp. 251-257). The number of EMs is in this case over 300 000.

(6) 2 Step b)

(7) The matrix with the calculated EMs E was obtained in step a). Similarly to matrix N, the matrix N.sub.p contains the stoichiometry, i.e. the stoichiometric coefficients, of the external metabolites. Possible macroreactions of the stoichiometric network were combined in matrix K with Formula 21:
K=N.sub.p.Math.E(Formula 21)
3 Step c)

(8) The measured data of the process was taken from Baughman et al., who gives various measured quantities of a fermentation of hybridoma cells over the course of a batch process (cf. FIG. 6) [On the dynamic modeling of mammalian cell metabolism and mAb production. In: Computers & Chemical Engineering (2010) 34 (2), pp. 210-222]. The measured data was entered into the method.

(9) 4 Step d)

(10) Using spline-interpolated measured values from c) (C.sup.int), the growth and death rates and the specific uptake and secretion rates were calculated (cf. FIG. 7). Lysis was incorporated with a predefined lysis factor K.sub.l=0.1, which was entered into the method and was constant over the process time. Shifting of the measured data was not necessary, because in this case the data comes from a batch process without further additions. Accordingly, the data shows a steady course, because all concentration changes are caused by the cells, and not by additions.

(11) Additional information is employed for calculating the rates {tilde over (q)}. Thus, with the aid of the total biomass (BM in

(12) $\left[\frac{C-\mathrm{mol}}{l}\right]),$
also given in the data set from Baughman et al., and the total cell count, an average C-mol content of

(13) $f_{C-\mathrm{mol},X_v}=18.41\left[\frac{C-\mathrm{mol}}{{10}^9\mathrm{cells}}\right]$
could be calculated. The C-mol-based growth rate could now be calculated from Formula 22:

(14) $\begin{array}{cc}\stackrel{~}{}\left[\frac{C-\mathrm{mol}}{h.Math.{10}^9\mathrm{cells}}\right]=\left[\frac{1}{h}\right].Math.f_{C-\mathrm{mol},X_v}\left[\frac{C-\mathrm{mol}}{{10}^9\mathrm{cells}}\right]& \left(\mathrm{Formula}22\right)\end{array}$

(15) The C-mol based rate of formation of the antibody can be estimated similarly. For this, the molar composition of the antibody was estimated as CH.sub.1.58O.sub.0.31N.sub.0.26S.sub.0.004 with a formal molar mass of

(16) $M_{\mathrm{mAb},C-\mathrm{mol}}=22.45\frac{g}{\mathrm{mol}}.$
Here, it is assumed that the molar composition corresponds to an average molar composition of proteins as indicated by Villadsen et al. [Bioreaction engineering principles (2011), Chapter 3, Elemental and Redox Balances, p. 73, Springer Verlag, ISBN: 978-1-4419-9687-9]. The molar mass of the whole antibody was estimated at

(17) $M_{\mathrm{mAb}}=150000\frac{g}{\mathrm{mol}}.$
The rate of formation of the antibody was then obtained from the formula:

(18) $\begin{array}{cc}{\tilde{q}}_{\mathrm{mAb}}\left[\frac{C-\mathrm{mol}}{h.Math.{10}^9\mathrm{cells}}\right]=q_{\mathrm{mAb}}\left[\frac{{10}^{-4}\mathrm{mol}}{h.Math.{10}^9\mathrm{cells}}\right].Math.\frac{M_{\mathrm{mAb},C-\mathrm{mol}}}{M_{\mathrm{mAb}}}.Math.{10}^4& \left(\mathrm{Formula}23\right)\end{array}$

(19) The course of the {tilde over (q)}(t) over time could then be employed for selecting the macroreactions.

(20) 5 Step e)

(21) In step e), an EM-subset of macroreactions was set up, with which the data set was reproduced as well as possible. This required the matrix K from step b). As the number of over 300 000 macroreactions would have led to an excessively large number of possible combinations, a data-dependent prereduction was carried out first.

(22) To this end, the rates {tilde over (q)}(t) determined in step d) were used for calculating the yield coefficients Y.sup.m for all combinations of two external metabolites. The lower limit of a yield coefficient Y.sub.i,j was selected such that 99% of the determined yield coefficients Y.sub.i,j.sup.m(t) are above this value. The upper limit was selected such that 99% of the determined yield coefficients Y.sub.i,j.sup.m(t) are below this value. By way of example, Table 2 shows some determined limits and also the proportion of EMs having yield coefficients Y.sub.i,j.sup.EM within these limits. Overall, the number of EMs could thus be reduced to approx. 3000.

(23) TABLE-US-00002 TABLE 2 External metabolites, their maximum and minimum yield coefficients Y.sub.i,j, and the proportion of EMs having yield coefficients within the specified limits External Lower Upper Proportion of components limit limit EMs within the limits Ala:Asn 14.6811 0.1752 99.9134% Asn:Glc 0.0293 16.1603 64.0488% Asp:Ala 1.1130 2.1011 74.0839%

(24) After the data-dependent reduction, a data-independent reduction was subsequently additionally carried out. In this case, a maximum value for the cosine similarity of two EMs of 0.995 was defined. Beginning with the first reaction, all macroreactions that exceeded this value were thus removed from matrix K. There remained approx. 500 macroreactions from matrix K (also called reduced matrix K), which still cover more than 95% of the volume of the solution space spanned by the approx. 3000 EMs.

(25) Before the selection process, a reconciliation of the components indicated in the metabolic network according to Niu et al. (which correspond to the external metabolites of the metabolic network from a)) with the measured concentrations of the components from c) was additionally carried out. Apart from proline, all concentrations measured by Baughman et al. are also taken into account in the metabolic network according to Niu et al. So as to be able to employ the measurement of the proline concentration, it would be possible either to use another simplified network that contains proline as an external metabolite, or it would be possible to expand the existing metabolic network.

(26) Components that did in fact occur in the calculated macroreactions, but for which no data were available, were also ignored in the following. The corresponding rows of matrix K were accordingly deleted from the matrix. Deletion of the corresponding rows does not mean that these inputs or outputs are not used by the cell. They still exist in the metabolic network, but no measurements are available with which they can be reconciled. In this example, the inputs or outputs of arginine, glutamate, glycine, histidine, leucine, lysine, methionine, ammonium, oxygen, phenylalanine, serine, threonine, tryptophan, tyrosine and valine were ignored.

(27) In the next steps of the method, the reduced matrix Kwhich represents the background knowledgeand the rates {tilde over (q)}(t) from d) and the measured data from c)which form the process know-howare then used for obtaining a smallest possible subset L of the macroreactions from K.

(28) The inventive linear estimation of reaction rates of selected macroreactions was used as quality criterion.

(29) As with the rates {tilde over (q)}(t), the measured values of the cell count and of the antibody were normalized here to C-mol. This is necessary so that the dimension of the macroreactions agrees with those of the measured values.

(30) The subset was selected with a genetic algorithm. In the calculation of the target function of this genetic algorithm, the linear optimization problem addressed in the linear estimation of reaction rates of selected macroreactions was solved. The final sum of the least error squares of the linear optimization problem calculated here was at the same time the value of the target function for the particular selection of the macroreactions.

(31) For selecting the size of the subset L from K, optimization was carried out repeatedly with a different number of macroreactions in L. The number represents a compromise between the complexity of the model and the accuracy of reproduction. To determine how many reactions are sufficient for reproduction, either selection of the subset L may be repeated for a varying number of macroreactions, or a penalty term for the number of reactions can be added directly to the target function of the genetic algorithm. In this case several optimizations were carried out with a predefined number of macroreactions (10, 7, 5, 4 and 3). The minimum error found with the genetic algorithm is plotted in FIG. 9 against the number of macroreactions. It was found that in this case fewer than seven macroreactions are too few for representing the process course sufficiently well. The selected macroreactions are given in Table 3.

(32) TABLE-US-00003 TABLE 3 Selected subset of the macroreactions (L). Components that are not underlined are not taken into account in the model, as no measurements are available for these. 0.474 Alanine + 0.474 Methionine .fwdarw. 0.158 Asparagine + 0.316 Aspartate + 0.632 Glycine + 0.158 Tryptophan 0.015 Alanine + 0.00789 Arginine + 0.0304 Asparagine + 0.0161Glucose + 0.0236 Glutamine + 0.00375 Glutamate + 0.00366 Histidine + 0.00953 Isoleucine + 0.015 Leucine + 0.112 Methionine + 0.00626 Phenalanine + 0.0154 Serine + 0.0109 Valine .fwdarw. 0.963 X(Biomass) + 0.00276Aspartate + 0.24 Glycine + 0.0208 Tryptophan 0.295Asparagine + 0.147 Glutamate .fwdarw. 0.295 Aspartate + 0.885 Glycine + 0.147Lactate 0.00753 Arginine + 0.113Asparagine + 0.0603 Glucose + 0.0225 Glutamine + 0.0824 Histidine + 0.00909 Isoleucine + 0.00597 Phenalanine + 0.0216 Tryptophan + 0.00431 Tyrosine + 0.0104 Valine .fwdarw. 0.918X (Biomass) + 0.061Alanine + 0.0865 Aspartate + 0.343 Glycine + 0.0631 Methionine 0.0654 Arginine + 0.412Aspartate + 0.00991 Glucose + 0.0145 Glutamine + 0.554 Glycine + 0.00226 Histidine + 0.00588 Isoleucine + 0.00926 Leucine + 0.00706 Lysine + 0.0649 Phenalanine + 0.0095 Serine + 0.00671 Valine .fwdarw. 0.594X(Biomass) + 0.049Alanine + 0.395 Asparagine + 0.0503 Threonine + 0.0388 Tryptophan 0.0077 Arginine + 0.179 Aspartate + 0.0157 Glucose + 0.104 Glutamine + 0.216 Glycine + 0.00357 Histidine + 0.00929 Isoleucine + 0.0146 Leucine + 0.0112 Lysine + 0.038 Tyrosine + 0.0106 Valine .fwdarw. 0.939 X(Biomass) + 0.0624 Alanine + 0.152 Asparagine + 0.0183 Tryptophan 0.0342 Arginine + 0.211Aspartate + 0.00762Glucose + 0.0195Glutamine + 0.244 Glycine + 0.00452 Histidine + 0.0546 Isoleucine + 0.0185 Leucine + 0.0171 Lysine + 0.00406 Methionine + 0.0178 Tyrosine + 0.0203 Valine .fwdarw. 0.457X(Biomass) + 0.804 IgG (Antibody) + 0.185 Asparagine + 0.0153 Tryptophan

(33) In the macroreactions shown, all external metabolites of the metabolic network from a) are indicated. However, only the underlined external metabolites form part of the model, as measured data from c) are only available for these.

(34) 6 Step f)

(35) For the selected set of macroreactions, the reaction rates over time were determined. In this example, using the inventive method linear estimation of reaction rates of selected macroreactions, the measured values shown in FIG. 10 were approximated by an estimation of the reaction rates r(t). The result of the method is a piecewise linear course of the individual (volumetric) reaction rates. By dividing by the interpolated course of the viable cell count X.sub.(t), the cell-specific reaction rates r(t) of the macroreactions shown in Table 3 were obtained. The reaction rates r(t) thus obtained are shown in FIG. 10.

(36) 7 Step g)

(37) For all macroreactions shown in Table 3, generic kinetics according to Formula 24 were assumed:

(38) $\begin{array}{cc}{\hat{r}}_k\left(t\right)=r_{k,\max }.Math.\underset{i=1}{\overset{N_l}{.Math.}}{\tilde{r}}_i(\underset{\_}{C}\left(t\right),\underset{\_}{p},\mathrm{.Math.})& \left(\mathrm{Formula}24\right)\end{array}$

(39) In this case, they were realized by Monod kinetics, i.e. for each reaction k for each substrate i, a limitation according to Formula 25:

(40) $\begin{array}{cc}{\tilde{r}}_i\left(t\right)={\left(\frac{C_i\left(t\right)}{K_{m,k,i}+C_i\left(t\right)}\right)}^{n_i}& \left(\mathrm{Formula}25\right)\end{array}$

(41) was introduced. Here, r.sub.k,max is the maximum reaction rate, N.sub.l is the number of limitations taken into account, C.sub.i is the concentration of the component i, K.sub.m,k,i are the associated Monod constants and n.sub.i is the Hill parameter for the reaction order. Their values are adjusted in steps h) and j). Further terms are found from the analysis of the reaction rates r(t) from f). In this example, in addition to substrate limitations, inhibitions according to Formula 26 were also taken into account.

(42) 0$\begin{array}{cc}{\tilde{r}}_i={\left(\frac{K_{I,k,i}}{K_{I,k,i}+C_i}\right)}^{n_i}& \left(\mathrm{Formula}26\right)\end{array}$

(43) For this limitation too, it was necessary to adjust the values of the parameters K.sub.l,k,i, and n.sub.i. The kinetic terms used for the reactions are given in Table 4.

(44) TABLE-US-00004 TABLE 4 Kinetic terms of the selected macroreactions from L ${\hat{r}}_1\left(t\right)=r_{1,\max }.Math.\left(\frac{\left[\mathrm{Ala}\right]\left(t\right)}{K_{m,\mathrm{Ala},1}+\left[\mathrm{Ala}\right]\left(t\right)}\right).Math.{\left(\frac{\left[\mathrm{Glc}\right]\left(t\right)}{K_{m,\mathrm{Glc},1}+\left[G\mathrm{lc}\right]\left(t\right)}\right)}^2.Math.{\left(\frac{K_{I,\mathrm{Asn},1}}{K_{I,\mathrm{Asn},1}+\left[\mathrm{Asn}\right]\left(t\right)}\right)}^2$ ${\hat{r}}_2\left(t\right)=r_{2,\max }.Math.\left(\frac{\left[\mathrm{Glc}\right]\left(t\right)}{K_{m,\mathrm{Glc},2}+\left[\mathrm{Glc}\right]\left(t\right)}\right).Math.\left(\frac{\left[\mathrm{Gln}\right]\left(t\right)}{K_{m,\mathrm{Gln},2}+\left[G\ln \right]\left(t\right)}\right).Math.\left(\frac{\left[\mathrm{Asn}\right]\left(t\right)}{K_{m,\mathrm{Asn},2}+\left[\mathrm{Asn}\right]\left(t\right)}\right).Math.{\left(\frac{\left[\mathrm{Ala}\right]\left(t\right)}{K_{m,\mathrm{Ala},2}+\left[\mathrm{Ala}\right]\left(t\right)}\right)}^2$ ${\hat{r}}_3\left(t\right)=r_{3,\max }.Math.\left(\frac{\left[\mathrm{Glc}\right]\left(t\right)}{K_{m,\mathrm{Glc},3}+\left[\mathrm{Glc}\right]\left(t\right)}\right).Math.\left(\frac{\left[\mathrm{Asn}\right]\left(t\right)}{K_{m,\mathrm{Gln},3}+\left[\mathrm{Asn}\right]\left(t\right)}\right).Math.{\left(\frac{K_{I,\mathrm{Lac},3}}{K_{I,\mathrm{Lac},3}+\left[\mathrm{Lac}\right]\left(t\right)}\right)}^2$ ${\hat{r}}_4\left(t\right)=r_{4,\max }.Math.\left(\frac{\left[\mathrm{Glc}\right]\left(t\right)}{K_{m,\mathrm{Glc},4}+\left[\mathrm{Glc}\right]\left(t\right)}\right).Math.\left(\frac{\left[\mathrm{Gln}\right]\left(t\right)}{K_{m,\mathrm{Gln},4}+\left[G\ln \right]\left(t\right)}\right).Math.\left(\frac{\left[\mathrm{Asn}\right]\left(t\right)}{K_{m,\mathrm{Asn},4}+\left[\mathrm{Asn}\right]\left(t\right)}\right).Math.\left(\frac{\left[X_t\right]\left(t\right)}{K_{m,\mathrm{Xt},4}+\left[X_t\right]\left(t\right)}\right)$ ${\hat{r}}_5\left(t\right)=r_{5,\max }.Math.\left(\frac{\left[\mathrm{Glc}\right]\left(t\right)}{K_{m,\mathrm{Glc},5}+\left[\mathrm{Glc}\right]\left(t\right)}\right).Math.\left(\frac{\left[\mathrm{Gln}\right]\left(t\right)}{K_{m,\mathrm{Gln},5}+\left[G\ln \right]\left(t\right)}\right).Math.\left(\frac{\left[\mathrm{Asp}\right]\left(t\right)}{K_{m,\mathrm{Asp},5}+\left[\mathrm{Asp}\right]\left(t\right)}\right).Math.\left(\frac{K_{I,\mathrm{Asp},5}}{K_{I,\mathrm{Asp},5}+\left[\mathrm{Asp}\right]\left(t\right)}\right).Math.\left(\frac{K_{I,\mathrm{Asn},5}}{K_{I,\mathrm{Asn},5}+\left[\mathrm{Asn}\right]\left(t\right)}\right)$ ${\hat{r}}_6\left(t\right)=r_{6,\max }.Math.\left(\frac{\left[\mathrm{Glc}\right]\left(t\right)}{K_{m,\mathrm{Glc},6}+\left[\mathrm{Glc}\right]\left(t\right)}\right).Math.\left(\frac{\left[\mathrm{Gln}\right]\left(t\right)}{K_{m,\mathrm{Gln},6}+\left[\mathrm{Gln}\right]\left(t\right)}\right).Math.{\left(\frac{\left[\mathrm{Asp}\right]\left(t\right)}{K_{m,\mathrm{Asp},6}+\left[\mathrm{Asp}\right]\left(t\right)}\right)}^2$ ${\hat{r}}_7\left(t\right)=r_{7,\max }.Math.\left(\frac{\left[\mathrm{Glc}\right]\left(t\right)}{K_{m,\mathrm{Glc},7}+\left[\mathrm{Glc}\right]\left(t\right)}\right).Math.\left(\frac{\left[\mathrm{Gln}\right]\left(t\right)}{K_{m,\mathrm{Gln},7}+\left[\mathrm{Gln}\right]\left(t\right)}\right).Math.\left(\frac{\left[\mathrm{Asp}\right]\left(t\right)}{K_{m,\mathrm{Asp},7}+\left[\mathrm{Asp}\right]\left(t\right)}\right)$

(45) 8 Step h)

(46) For each reaction rate, the course of the reaction rate {circumflex over (r)}.sub.i(p, C.sup.int(t)) could be calculated algebraically with the kinetics given in Table 4 and the interpolated values of the concentrations C.sup.int(t) taken into account in the kinetics.

(47) The parameters of these kinetics were adjusted separately for each reaction i to the reaction rate r.sub.i(t) determined in step f). The target function for optimization of the parameters occurring in reaction i was in this example:

(48) $\begin{array}{cc}\underset{{\underset{\_}{p}}_k}{\min }\left(\underset{l=0}{\overset{T}{.Math.}}{\left({\hat{r}}_k\left({\underset{\_}{p}}_k,{\underset{\_}{C}}^{\mathrm{int}}\left(t_l\right)\right)-r_k\left(t_l\right)\right)}^2\right)& \left(\mathrm{Formula}27\right)\end{array}$

(49) The courses of all calculated {circumflex over (r)}.sub.k(p.sub.k, C.sup.int(t)) adjusted in this way are shown together with the corresponding r.sub.k(t) in FIG. 11. The courses of the former are shown with dashes, and those of the latter are shown with solid lines. It can be seen that the course agrees qualitatively. This means that with the selected kinetics, the dynamics of the process can also be reproduced satisfactorily. This information is very useful in this modeling step, because if reproduction is unsatisfactory, the quick steps g) (selection of other kinetics) and h) (estimation of parameter values) can be repeated, until the desired degree of adjustment is reached. Thus, step i) was not necessary here.

(50) 9 Step j)

(51) Further adjustment of the model parameter values p was carried out with the measured data from c). For this, all parameters were optimized at the same time. Moreover, the processes of apoptosis and lysis, which have not been examined previously, were also included. These are required in the differential equations that describe the development of the viable cell count and total cell count:

(52) $\begin{array}{cc}\frac{d{X}_v}{dt}=\left({}_x-{}_d\right).Math.X_v& \left(\mathrm{Formula}28\right)\\ \frac{d{X}_t}{dt}={}_x.Math.X_v-K_l.Math.\left(X_t-X_v\right)& \left(\mathrm{Formula}29\right)\end{array}$

(53) The selected kinetics for describing apoptosis was:

(54) 0$\begin{array}{cc}{}_d\left(t\right)={}_{d,\max }.Math.\frac{\left(\left[\mathrm{Lac}\right]\left(t\right)-C_{\mathrm{Lac},\mathrm{cr}}\right)}{K_{d,\mathrm{Lac}}+\left(\left[\mathrm{Lac}\right]\left(t\right)-C_{\mathrm{Lac},\mathrm{cr}}\right)},\left[\mathrm{Lac}\right]C_{\mathrm{Lac},\mathrm{cr}}& \left(\mathrm{Formula}30\right)\\ {}_d\left(t\right)=0,\left[\mathrm{Lac}\right]<C_{\mathrm{Lac},\mathrm{cr}}& \left(\mathrm{Formula}31\right)\end{array}$

(55) The lysis rate K.sub.l was assumed to be constant over the process. In addition to the parameters of the reaction rates, the parameters C.sub.Lac,cr (critical lactate concentration), .sub.d,max (maximum death rate), K.sub.d,Lac (Monod parameter for describing the influence of the lactate concentration on the death rate) and K.sub.l (lysis rate) introduced by apoptosis and lysis were determined in this step. In the example, the course of the estimated concentrations (t) was determined from the starting values of the data set, by numerical solution of the ODE system. The difference between the measured concentrations C.sup.m(t) and the estimated concentrations (t) was minimized by usual methods with the following target function:

(56) $\begin{array}{cc}\underset{\underset{\_}{p}}{\min }\left(\underset{i=1}{\overset{n_{\mathrm{comp}}}{.Math.}}\left(\underset{l=0}{\overset{T}{.Math.}}{\left({\hat{C}}_i\left(\underset{\_}{p},t_l\right)-C_l^m\left(t_l\right)\right)}^2\right)\right)& \left(\mathrm{Formula}32\right)\end{array}$

(57) With a total of 33 parameters p, as a rule this optimization is difficult to perform, as the target function has many local optima. If a deterministic optimization algorithm is started, for example the Levenberg-Marquardt algorithm on the starting values of the parameters known from step h), there is a much greater chance of success. The adjusted process course is shown in FIG. 12. The adjusted parameters are shown in Table 5.

(58) TABLE-US-00005 TABLE 5 Parameters of the kinetics and of apoptosis and lysis K.sub.m,Glc,1 14.6 K.sub.m,Gln,7 0.0187 K.sub.m,Ala,1 3.41 K.sub.m,Asp,7 0.872 K.sub.m,Glc,2 0.0508 r.sub.1,max 9.47 K.sub.m,Gln,2 0.00881 r.sub.2,max 9.91 K.sub.m,Asn,2 1.38 r.sub.3,max 57.6 K.sub.m,Ala,2 2.19 r.sub.4,max 21.7 K.sub.m,Glc,3 7.13 r.sub.5,max 0.345 K.sub.m,Asn,3 6.84 r.sub.6,max 49.4 K.sub.m,Xt,4 0.0315 r.sub.7,max 3.03 K.sub.m,Glc,4 1.29 K.sub.IAsn,1 16.1 K.sub.m,Gln,4 2.19 K.sub.I,Lac,3 0.681 K.sub.m,Asn,4 1.68 K.sub.I,Asp,5 9.74 K.sub.m,Glc,5 100 K.sub.I,Asn,5 1.10 K.sub.m,Gln,5 28.2 .sub.d,max 0.125 K.sub.m,Asp,5 102 K.sub.d,Lac 1.01 K.sub.m,Glc,6 0.0451 C.sub.Lac,cr 1.22 K.sub.m,Gln,6 0.791 K.sub.l 0.00843 K.sub.m,Asp,6 1.06 K.sub.m,Glc,7 0.0145
10 Step k)

(59) The model, consisting of the matrix L, the kinetics from Table 4 and the kinetics of apoptosis with the associated parameter values from Table 5, was produced as the output.

LIST OF SYMBOLS

(60) TABLE-US-00006 _ (underline) denotes a vector .sub.i (subscript i) denotes the i-th element of a vector .sub.k (subscript k) denotes the k-th element of a vector [ ] denotes the concentration of the component in brackets C concentration C concentration difference C.sup.Int interpolated concentration estimated concentration (e.g. by solving a differential equation) C.sub.s shifted concentration C.sub.cr critical concentration C.sup.m measured concentration D dilution rate q determined cell-specific secretion and uptake rate {tilde over (q)} determined cell-specific secretion and uptake rate that has been converted from any unit to $\left[\frac{\mathrm{Substance}\mathrm{amount}}{\mathrm{Time}.Math.\mathrm{Cell}\mathrm{count}}\right]$ r determined reaction rate {circumflex over (r)} estimated reaction rate (e.g. by calculating reaction kinetics) {tilde over (r)} limitation of kinetics r.sub.max parameter of reaction kinetics N stoichiometric matrix N.sub.p external stoichiometric matrix K matrix that contains macroreactions E matrix that contains all elementary modes X.sub.t total cell count X.sub.v viable cell count growth rate .sub.d death rate {tilde over ()} growth rate that has been converted from any unit to $\left[\frac{\mathrm{Substance}\mathrm{amount}}{\mathrm{Time}.Math.\mathrm{Cell}\mathrm{count}}\right]$ K.sub.d lysis rate K.sub.I parameter of an inhibition limitation K.sub.M parameter of a substrate limitation n Hill parameter of an inhibition or substrate limitation L subset of the macroreactions that is used for the model p model parameter S substrate SSR.sub.q sum of squared residuals of the specific uptake or secretion rates SSR.sub.C sum of squared residuals of the concentration SSR.sub.r sum of squared residuals of the reaction rates

BRIEF DESCRIPTION OF THE DRAWINGS

(61) FIG. 1 shows the shift of measured data: It shows the actual course of a measured quantity (C.sub.i(t)), which changes suddenly when there are changes of the dilution rate (D(t)). The course of the shifted concentration (C.sub.i,s(t)) only comes from changes caused by the cell.

(62) FIG. 2 shows the flux map of two specific rates q.sub.1 and q.sub.2. The contour lines indicate the frequency with which the particular combination of the rates occurs in the measured data.

(63) FIG. 3 shows a three-dimensional representation of the solution space, which is spanned by a positive linear combination of EMs. The solution space of the complete set is shown in black, and that of a subset is shown in grey.

(64) FIG. 4 shows the flux map of two specific rates q.sub.1 and q.sub.2. The 2-dimensional projections of the macroreactions of a set L are shown as vectors.

(65) FIG. 5 shows a schematic representation of the metabolic network from Niu et al. Here, the boundary of the cell is shown as a box. The intracellular border of the mitochondrium is shown with a dashed line. External components are marked with the subscript xt. The arrows and dotted arrows denote reactions.

(66) FIG. 6 shows the measured data of a fermentation with hybridoma cells from Baughman et al. The total cell count (total cells) is calculated here from the total of viable cells and dead cells. The abbreviations GLC, GLN, ASP, ASN, LAC, ALA and PRO denote the substrate glucose and the amino acids glutamine, aspartic acid, asparagine, alanine and proline and the metabolic product lactate. The abbreviation MAB denotes the product of monoclonal antibodies and BM denotes biomass.

(67) FIG. 7 shows the growth and death rates and the cell-specific uptake and secretion rates. All cell-specific rates except q.sub.MAB are given in

(68) $\left[\frac{\mathrm{mM}}{h.Math.{10}^9\mathrm{Cells}}\right].$
The rate q.sub.MAB is given in

(69) $\left[\frac{{10}^{-4}\mathrm{mM}}{h.Math.{10}^9\mathrm{Cells}}\right].$

(70) FIG. 8 shows the concentrations approximated with the linear estimation of reaction rates of selected macroreactions with the selected reaction set. The total cell count (X.sub.t) and the antibody concentration (MAB) were converted for this to C-mol.

(71) FIG. 9 shows the minimum error plotted against the number of macroreactions in the subset (n.sub.R).

(72) FIG. 10 shows the reaction rates of the macroreactions r(t) determined with the inventive method linear estimation of reaction rates of selected macroreactions.

(73) FIG. 11 shows the reaction rates of the macroreactions r(t) determined with the inventive method linear estimation of reaction rates of selected macroreactions (solid line) together with the algebraically calculated reaction rates {circumflex over (r)}(p, C.sup.int(t)) (dashed line).

(74) FIG. 12 shows a comparison of the measured concentrations Cm (t) (points) and the simulated process course (t) (solid line). The concentrations are given in [mM]. The viable cell count and total cell count (X.sub.v/X.sub.t in [10.sup.9 cells/l]) and the concentration of the antibody (mAb in [10.sup.4 mM]) are exceptions.