CN107408161A

CN107408161A - 用于创建发酵模型的计算机实现的方法

Info

Publication number: CN107408161A
Application number: CN201680007741.XA
Authority: CN
Inventors: T.奈曼; L.赫宾; S.恩格尔
Original assignee: Bayer AG
Current assignee: Bayer AG
Priority date: 2015-01-29
Filing date: 2016-01-28
Publication date: 2017-11-28
Anticipated expiration: 2036-01-28
Also published as: IL253584A0; EP3051449A1; US10872680B2; EP3251039A1; EA201791659A1; TW201643744A; AU2016212059B2; US10296708B2; JP6816003B2; WO2016120361A1; BR112017016198A2; TWI690813B; CA2975012A1; SG10202006972VA; AR103564A1; CN107408161B; US20160224721A1; CA2975012C; HK1247340A1; JP2018511849A

Abstract

本申请涉及一种用于利用生物体基于测量数据来创建生物反应（发酵或者全细胞催化）的模型的计算机实现的方法。

Description

用于创建发酵模型的计算机实现的方法

技术领域

本发明涉及一种用于在使用生物体的情况下创建生物反应（尤其是发酵或者全细胞催化）的模型的计算机实现的方法。

就本申请而言的生物体是由植物或动物细胞（如哺乳动物细胞、酵母菌、细菌、藻类等等）构成的培养，所述培养被用在生物反应中。

对发酵过程的以传感器技术方式的监控以及例如借助于Bayer技术服务有限公司的质量源于设计分析自动化平台BaychroMAT®对来自过程的样品的分析实时地提供了关于在生物反应器中的过程的状态的不同的信息。通常，细胞数目、细胞活力、基质的浓度，如碳源（例如葡萄糖）、氨基酸或者O₂、产品和副产品（例如乳酸或CO₂）、过程参数（如温度和/或pH值）或者产品特征被确定。这些数据可以通过所计算的数据和/或例如来自现有技术的外插法来补充。这些数据一起形成就本申请而言的测量数据或者过程知识（Prozesswissen）。

背景技术

就本申请而言，关于生物体的背景知识指的是关于生物体的生物化学反应（特定的和不特定的反应）以及尤其是细胞内部的反应或用于描述生物体特定的新陈代谢网（也称作SN或metabolische Netzwerke）的宏观反应的知识，所述新陈代谢网由基质、代谢物（也称作新陈代谢网的节点）、产品以及在它们之间的生物化学反应组成。这些生物化学反应通过它们的：

（a）化学计量，

（b）可逆性（在生物条件下），

（c）到化学计量网中的结合

来限定。

到目前为止，所述测量数据主要被用于对过程的定性监控。在下文，示出了对如下技术提问的选择，对于所述技术提问的解决方案来说，动态过程模型是必要的。

就本申请而言对过程知识的技术应用提供了对在生物反应器中的过程的基于模型的状态估计。如“扩展卡尔曼滤波器（extended Kalman Filter）”那样的方法能够实现对过程参量的连续估计，关于所述过程参量不连续地存在测量[Welch G，Bishop G. 1995年的An Introduction to the Kalman Filter. Chapel Hill, NC, USA: University ofNorth Carolina at Chapel Hill]。也可以根据其它测量来计算不能测量的参量的变化过程。为此的前提是正确地描述所基于的过程的过程模型。

另一应用是基于模型的、最优的过程控制。在这种情况下，使用动态过程模型，以便在基于模型的、预测性的闭合调节回路中使过程控制在产量、产品特征或副产品的形成或其它目标参量方面优化。所述动态过程模型例如由Frahm等人针对杂交瘤细胞培养来示出[Frahm B, Lane P, Atzert H, Munack A, Hoffmann M, Hass VC, Portner R. 2002年的Adaptive, Model-Based Control by the Open-Loop-Feedback-Optimal (OLFO)Controller for the Effective Fed-Batch Cultivation of Hybridoma Cells.Biotechnol. Prog. 18(5):第1095–1103页]。

对于两种提到的技术应用来说重要的是，所创建的过程模型具有尽可能低的复杂性，即同时在对过程的重现的良好的精确度的情况下的有限数目的状态变量和/或方程。

除了所提到的用于过程控制的应用之外，也可以在过程研发期间使用动态过程模型，以便计划具有最优的信息收益的实验。所述行为被称作基于模型的实验计划[Franceschini G, Macchietto S. 2008年的Model-based design of experiments forparameter precision: State of the art. Chemical Engineering Science 63(19):第4846–4872页]。除了上面提到的针对模型的复杂性的前提之外，对于所述技术应用还必要的是，动态过程模型在研发阶段已经存在。所述动态过程模型应该能够尽可能快地根据已经存在的过程知识来创建，以便使过程研发的时间花费保持得低。

因而，已经存在如下需求：提供一种方法，所述方法能够实现在使用背景知识和过程知识的情况下创建动态过程模型。为此，为了能够将该模型例如用于状态估计、最优的过程控制或者用于基于模型的实验计划，模型的复杂性必须低。相关性、即过程参量或过程状态对过程表现的影响应该在设计空间（Design Space）之内足够精确地量化。为此，应该使用所有现存的关于过程状态的信息。对产品特征的基于模型的描述应该能在有需求时集成到所述模型中。在其中存在过程知识的区域被称作设计空间。该方法应该能被应用到上面列举的生物反应上并且这种动态模型的研发时间应该显著地被缩短。到目前为止的用于研发动态模型的方案需要几个月直至几年，直至将过程模型完成。按照本发明，本方案将研发时间减少到几周。

就本申请而言的典型的产品特征例如是蛋白质的糖基化模式或者蛋白质完整性，而对此不进行限制。被用在上面提到的上下文中的动态模型迄今为止不具有该特性。本方案能够实现对产品特征的简单的基于模型的集成。

对发酵的基于模型的过程控制由Frahm等人以杂交瘤细胞培养为例来示出（FrahmB, Lane P, Atzert H, Munack A, Hoffmann M, Hass VC, Portner R. 2002年的Adaptive, Model-Based Control by the Open-Loop-Feedback-Optimal (OLFO)Controller for the Effective Fed-Batch Cultivation of Hybridoma Cells.Biotechnol. Prog. 18(5):第1095–1103页）。基本的过程参量这里基于模型地来控制。这里，不进行对产品特征的集成。细胞的数学模型已经针对所述特定的过程来设计并且只能以高花费被转用到具有相同或不同生物体或者相同生物体的品系的过程上。以细胞内部的反应为形式的背景知识在该模型中没有明确地予以考虑。这里，只能以非常高的花费来实现将其它测量参量集成到该模型中并且因此完整地使用关于过程状态的信息。因此，该方案是单独的解决方案，所述单独的解决方案又能被转用到其它过程上，还能够实现对所确定的数据的完整使用。所提到的方法由于该模型的所要期望的研发时间以及该解决方案到其它具有相同生物体或具有不同生物体的过程上的花费高的可转用性而没有解决上面提到的技术问题。

在Kontoravdi等人的出版物中找到进一步的建模，所述建模也一并包括产品特征（如糖基化）。描述主新陈代谢的模型没有把以细胞内部的反应为形式的背景知识包括在内，而且也不能被转用到其它具有相同或不同生物体的过程上。这里，不能将其它测量参量集成到该模型中[Kontoravdi C, Asprey SP, Pistikopoulos EN, Mantalaris A. 2007年的Development of a dynamic model of monoclonal antibody production andglycosylation for product quality monitoring. Computers & ChemicalEngineering 31(5-6):第392–400页]。该方法也不能实现对关于过程状态的信息的完整使用，该方法需要对该模型的长的研发时间并且不能被转用到不同生物体或品系上。因此，该方法不是所述技术问题的解决方案。

Jedrzejewski等人和Jimenez等人的在结合核苷酸新陈代谢的情况下对糖基化的建模把以内部新陈代谢中间物的平衡方程为形式的背景知识包括在内 [JedrzejewskiPM, del Val, Ioscani Jimenez, Constantinou A, Dell A, Haslam SM, Polizzi KM,Kontoravdi C. 2014年的Towards Controlling the Glycoform: A Model FrameworkLinking Extracellular Metabolites to Antibody Glycosylation. Internationaljournal of molecular sciences 15(3):第4492–4522页; Jimenez del Val, Ioscani,Nagy JM, Kontoravdi C. 2011年的A dynamic mathematical model for monoclonalantibody N-linked glycosylation and nucleotide sugar donor transport within amaturing Golgi apparatus. Biotechnology progress 27(6):第1730–1743页]。但是，在将该模型用于过程控制时，整个模型的复杂性以及细胞内部的新陈代谢中间物的缺乏的可观察性都是不利的。此外，主新陈代谢的模型不容许转用到其它过程上或者不容许对关于过程状态的信息的完整使用。因此，该方法不是所述技术问题的解决方案。

由Leifheit等人从事针对生物过程的灵活的模型生成[Leifheit J, Heine T,Kawohl M, King R. 2007年的对生物技术过程的计算机辅助的半自动化建模（Semiautomatic Modeling of Biotechnical Processes）在自动化技术55(5)]。借助于过程知识、不过在没有背景知识的情况下实现模型生成。该步骤能被用于具有相同或者不同生物体的不同的过程。这里，基础是如下宏观反应，操作者自己预先给定所述宏观反应。所述宏观反应的精确的化学计量在该方法中被确定。该方法针对小数目的状态参量或测量参量来描述。其它状态参量或测量参量的集成会与该方法的复杂性的显著升高相关联。在使用广泛的数据基础时，如所述数据基础例如通过BaychroMAT®平台来提供的那样，该方法会不再能被执行。该方法不能实现产品特征的集成。因此，该方法不是上面提到的技术问题的解决方案。

由Provost来描述对以宏观反应为形式的背景知识的使用，所述宏观反应作为基本模式（EM）从生物体的已知的新陈代谢（化学计量）网获得[Provost A. 2006. Metabolicdesign of dynamic bioreaction models. Faculté des Sciences Appliquées,Université catholique de Louvain, Louvain-la-Neuve, Louvain-la-Neuve, 第 81页及其后，第107 页及其后，第118页及其后]。该方法能被用于不同生物体或者相同生物体的品系。在使用过程知识的情况下，选择用于过程模型的宏观反应。然而，被限定用于预先限定的数目的宏观反应的过程部分单独地随机地被选择。所描述的方法这里提供了基本模式的多种可能的组合。宏观反应的数目以及因此模型复杂性被规定并且不能被改变。该方法得到针对每个过程部分的单独的模型。在考虑所选择的宏观反应的化学计量的情况下，选择各个宏观反应的动力学。然而，所述动力学的参数（模型参数）没有与过程数据适配。作为替代，对单独的过程部分模型的使用描绘了在过程数据方面的改变。虽然也可以基于广泛的数据基础来随机地选择反应，但是通过所描述的用于选择动力学的方案以及所选择的动力学不能描绘该过程的变化过程或生物体在该过程中的表现。此外，多个过程部分模型的使用还导致过程模型的不必要的复杂性提高。相关性、即过程参量或过程状态对过程表现的影响利用该方法没有被量化。这里，也不进行产品特征的集成。因此，该方法不是上面提到的技术问题的解决方案。

发明内容

因而，已经存在如下需求：提供一种方法，所述方法能够实现在考虑背景知识的情况下基于过程知识和过程测量数据快速并且高效地提供模型并且优化产品周转和关键的产品特征，而且不具有上面提到的缺点。

该任务是通过一种如下面描述的那样的用于创建利用生物反应器中的生物体进行生物反应的模型的方法来解决的。

本申请的主题是一种用于创建利用生物体的生物反应（尤其是发酵或者全细胞催化）的模型的计算机实现的方法，所述方法包括下列步骤：

a. 生物体的所选择的新陈代谢路径、所述新陈代谢路径的化学计量以及可逆性特性作为背景知识被输入到该方法中。换句话说，所述生物体的一个或多个新陈代谢网被输入到该方法中。基本模式（EM）根据该输入来计算；

b. 所述EM在矩阵K中合并，其中所述EM将来自a)的新陈代谢路径合并成宏观反应。就此，该矩阵K包含来自背景知识的所有可能的宏观反应的化学计量和可逆性特性；

c. 用于利用生物体的生物反应的测量数据（也称作过程知识）被输入；

d. 借助于插值法，基于来自c)的所输入的测量数据来计算所输入的新陈代谢路径的对于生物体来说特定的速率（一个或多个输入参量和输出参量中的分泌和摄取速率）。优选地，生物体的生长速率也被计算，特别优选地，生物体的衰亡速率也被计算；

e. 重要的宏观反应以来自b)的基本模式的子集的形式通过如下方式来选择：

i. 与数据无关地和/或与数据有关地预先减少来自b)的EM的数目，

ii. 借助于根据数学质量评定标准的算法和子集在矩阵L中的合并，利用来自c)的测量数据和/或来自d)的一个或多个速率、优选地利用来自c)的测量数据，从来自e) i.的所述预先减少中选择子集，

iii. 可选地以图形方式呈现该子集；

f. 借助于插值法，基于来自c)的所输入的测量数据和/或来自d)的速率来计算子集r(t)的宏观反应的反应速率；

g. 来自e) ii.的子集的宏观反应的动力学以下面的中间步骤来设计；由此限定了模型参数：

i. 根据宏观反应的化学计量来设计通用动力学，

ii. 对宏观反应的影响参量根据来自f)的反应速率来确定，

iii. 来自g) i.的通用动力学被扩展了如下项，所述项使在g) ii.中确定的影响参量量化；

h. 可选地，针对来自g)的动力学，使每个宏观反应的模型参数值单独地与来自f)的所计算的反应速率进行首次适配，并且检查适配质量；

i. 可选地，重复步骤g)和h)，直至达到了预先限定的适配质量；

j. 模型参数值与来自c)的测量数据适配；

k. 矩阵L、来自g)的动力学以及来自j)的模型参数值形成该模型并且被输出和/或被转用到过程控制或过程研发模块中。

通常，过程控制模块与过程控制系统在线（on-line）进行通信，所述过程控制系统通常被用于控制生物反应器。

通常，过程研发模块被用于线下（off-line）优化过程或者被用于计划实验。

按照本发明的对生物反应的建模基本上建立在有代表性的宏观反应的假定上，所述有代表性的宏观反应简化地呈现了内部新陈代谢过程。对于反应的选择来说，不仅需要生物化学背景知识而且需要过程知识。

在该方法的第一步骤中，由用户通过用户界面或理想地自动通过选择生物体和该生物体的来自数据库模块的所存储的新陈代谢路径来输入新陈代谢网的反应，所述反应的化学计量和可逆性特性，关于生物体的背景知识存储在所述数据库模块中。新陈代谢网（在现有技术中也称作化学计量网）和所述新陈代谢网的各个反应的特性是该生物体的背景知识。优选地，新陈代谢网包含来自对于生物体来说重要的新陈代谢路径的反应、例如糖酵解的反应。特别优选地，该选择包含外部反应。就本申请而言的外部反应包括至少一个在细胞之外的成分，通常包含至少一个输入参量和/或至少一个输出参量（产品、副产品等等）。特别优选地，新陈代谢网包含如下反应，所述反应描述了细胞生长、例如以从内部代谢物到外部生物量的经简化的反应的形式的细胞生长。在该例子中的图5和表1描述了能应用的新陈代谢网，而不对此进行限制。

接着，根据该生物体的所输入的、在一个或多个化学计量网中合并的新陈代谢路径来计算基本模式。每个基本模式都是来自所述新陈代谢路径（也就是说新陈代谢网的内部和外部反应）的反应速率的线性组合，其中对于内部代谢物来说不仅满足“稳态（steadystate）”条件，而且考虑反应的可逆性或不可逆性。在反应的考虑内部代谢物的“稳态”条件的线性组合的情况下，内部代谢物不能积累。

就本申请而言的内部反应仅仅在细胞之外进行。通过内部的成分的外部化、也就是说通过将原来内部的成分归类为输入或输出参量，可能的是，将与经外部化的、内部的成分相关联的内部反应建模为外部反应并且因此在这种情况下绕开内部代谢物的“稳态”条件。就本申请而言的宏观反应将如下所有反应合并，所述所有反应从一个或多个输入参量导致一个或多个输出参量。因此，每个基本模式都描述了一个宏观反应。就本申请而言，与Leifheit等人的方法相比，基于所输入的背景知识来确定宏观反应。

基本模式（EM）在矩阵E中合并、优选地在用于矩阵构造的模块中合并，所述模块以相对应的算法来配置。可以将已知的算法用于构造基本模式矩阵（Elementary ModesMatrix）。提到METATOOL作为例子，而对此不进行限制：[Pfeiffer T, Montero F,Schuster S. 1999. METATOOL: for studying metabolic networks. Bioinformatics15(3):第251-257页]。利用METATOOL来生成第一矩阵E，所述第一矩阵E描述了所输入的内部和外部反应。

在步骤b)中，借助于（外部）化学计量矩阵N_p从矩阵（E）生成由可能的宏观反应K组成的矩阵：

K=N_p·E （公式1）。

矩阵E朝K的转换从Provost公知[Provost A. 2006年的Metabolic design ofdynamic bioreaction models. Faculté des Sciences Appliquées, Universitécatholique de Louvain, Louvain-la-Neuve, Louvain-la-Neuve, 第81页]。

矩阵K的列向量描述了宏观反应。行向量描述了宏观反应的成分（输入和输出参量）。在矩阵K中记入了宏观反应的化学计量。每个就新陈代谢网而言可能的反应速率都可以作为这些宏观反应的线性组合来示出。

在现有技术中公知将EM用作过程模型的基础（例如从Provost A. 2006年的Metabolic design of dynamic bioreaction models. Faculté des Sciences Appliquées, Université catholique de Louvain, Louvain-la-Neuve, Louvain-la-Neuve, 第87页, 第118页及其后，和 Gao, J.等人(2007年)的 Dynamic metabolic modeling for aMAB bioprocess. Biotechnology progress, 23(1), 第168-181页）。

在另一步骤c)中，输入用于利用生物体的生物反应的可支配的测量数据（过程知识）。通常，细胞数目、细胞活力、基质的浓度，如碳源（例如葡萄糖）、氨基酸或者O₂、产品和副产品（例如乳酸或CO₂）、过程参数（如温度和/或pH值）或者产品特征被确定。所述输入可以由用户手动地实现或者自动地实现，诸如通过从用于存储测量数据的数据库模块中选择并且将所选择的数据转换成数据分析模块来实现，所述数据分析模块与数据库模块连接。

在步骤d)中，根据这些测量数据，计算基质和（副）产品的细胞特定的分泌和摄取速率（一起称作特定速率q(t)），以及可选地计算生物体的生长和衰亡速率(μ(t)，μ_d(t))。所述计算的前提是借助于插值法对活细胞数目、总细胞数目以及介质浓度的插值。根据这些可以确定测量参量的随时间的变化。所计算的速率q(t)、μ(t)、μ_d(t)给出了关于该生物体的所要观察的在时间上的动态表现的情况。

针对上面提到的速率的计算可以组合地使用一个或多个不同的插值法。

示范性地，借助于经样条插值的测量数据，从测量数据中，Leifheit等人描述了对测量参量的随时间的变化（例如总细胞数目的随时间的变化、活细胞数目的随时间的变化或者其它介质浓度的随时间的变化）的确定[Leifheit, J., Heine, T., Kawohl, M., &King, R. (2007年), 对生物技术过程的计算机辅助的半自动化建模（SemiautomaticModeling of Biotechnical Processes）在自动化技术55(5), 第211-218页]。就此，该方法通过参考被集成到本申请中。

上面提到的速率q(t)、μ(t)、μ_d(t)根据这些随时间的变化来计算：

例如，生物体的生长速率μ(t)可以根据总细胞数目X_t(t)和活细胞数目X_v(t)的经样条插值的值以及总细胞数目的能从中计算的随时间的变化以公式2来计算：

（公式2），

其中，D(t)是稀释速率。

在μ(t)的已知的变化过程的情况下，衰亡速率μ_d(t)可以根据X_v(t)的变化过程以及活细胞数目的随时间的变化的变化过程以公式3来计算：

（公式3）。

对另一成分i的特定速率q_i(t)的计算可以根据活细胞数目的经样条插值的值X_v(t)和成分的浓度C_i(t)以及随时间的变化的变化过程以公式4来实现，所述变化过程可以根据C_i(t)的经样条插值的值来确定：

（公式4）。

在该方法的一个优选的实施方式中，来自步骤c)的测量数据在第一次插值之前按如下地来整理：为了考虑所有不是由细胞造成的浓度变化并且根据测量数据获得浓度变化的连续的变化过程，测量数据被移动（在本申请中称作shiften（变动））。浓度测量被移动了的数值ΔC_i(t)可以根据公式5来计算：

（公式5），

其中，D(τ)是稀释速率。接着，根据公式6得到被变动的浓度变化过程C_i,s(t)：

（公式6），

因此，根据公式4和6得到说明被变动的浓度C_i,s(t)的变化过程的微分方程：

（公式7），

对测量数据的所述整理（变动）防止了尤其是在补料分批培养过程中投放或停止饲料（饲料高峰（Feedpeak））时的所计算的特定速率的突变。

图1示出了就本申请而言的对测量数据的整理/变动。

在该方法的一个特别的实施方式中，所整理的数据接着被用于利用Leifheit等人的方法计算总细胞数目的斜度。该斜度按照微分方程8以样条插值来近似：

（公式8），

特别优选地，细胞溶解在计算中依据细胞溶解因数K_l而被包括在内（公式9）。所述细胞溶解因数K_l例如可以被假定为在过程变化上恒定，

（公式9），

这样，被变动的总细胞数目X_t,s(t)的降低可以通过细胞溶解来解释，由此可以避免生长速率μ(t)的负值。

优选地，所整理的数据也被用于计算衰亡速率μ_d(t)。

优选地，在如示例性地在图2中的那样的通量图（Flux-Map Diagramm）中示出了特定速率q(t)的可能的组合。该图示能够实现关于所计算的特定速率q(t)的良好的一览。这里，等高线说明了：哪些区域在生理上是重要的。

如果特定速率q(t)以及可选地其它速率μ(t)和μ_d(t)都具有不同的量级和单位，那么这些通常借助于简化来合并成具有相同单位的特定的速率向量。例如，以克为单位[g]来测量的宏观分子的特定速率将具有单位。如果该宏观分子的成分例如基于它的C-mol(摩尔)-含量来估计，那么它的特定速率可以从[g]改变到[C-mol]地来示出，使得所述特定速率具有单位。特定速率形成针对该方法的另一步骤e)的基础之一、即对重要的宏观反应的选择。

在步骤e)中，基于所述数据来选择EM的子集（L），利用所述子集（L），根据数学质量评定标准可以良好地描绘来自d)的特定速率和/或来自c)的测量数据。EM在子集（L）中的数目应该尽可能小，以便确保过程模型的尽可能低的复杂性。但是子集L应该保证对过程知识的良好的描述。与来自a)的原始EM集合（K）相比，EM的选择缩小了解空间，但是继续包含细胞的所确定的生理上重要的区域。

图3示出了解空间的图示，其中EM的原始集合（K）被减小到子集（L）。

对于步骤e)来说，所计算的特定速率以及来自c)的测量数据通常被转换成用于选择重要的宏观反应的模块，所述模块以相对应的算法来配置。

在步骤e) i)中，以任意的顺序与数据无关地和/或与数据有关地预先减少矩阵K：

优选地，通过几何减少来实现与数据无关的预先减少。在此，针对随机选择的EM计算与所有其它模式的所有余弦相似性。具有最高相似性的EM从集合中被除去。重复该行为方式，直至达到了EM的预先限定的数目。所希望的数目通常针对该方法在准备阶段被限定。可以将解空间的容积用作控制变量。令人感到意外地，已经确定的是，与原始容积相比，宏观反应的数目的明显减少在保持所展开的容积的90至98%、优选地92至95%的情况下是可能的。

与数据有关的预先减少可以通过将EM的收益系数（Y^EM）与根据来自d)的特定速率计算的收益系数（Y^m）进行比较来实现。第k个EM的收益系数根据公式10通过除以外部代谢物i和j的相对应的化学计量系数来确定。对于第k个EM来说，这是矩阵录入项K_i,k和K_j,k：

（公式10）。

如果化学计量系数K_j,k=0，那么收益系数不能被确定。收益系数根据公式11说明了两个所测量的或者根据d)计算的细胞特定速率彼此间的比例：

（公式11）。

根据收益系数Y^m可以针对两个外部成分i和j的每个可能的组合来确定上极限和下极限。例如，两个外部代谢物i和j的最小收益系数可以被用作下极限，而的最大值可以被用作上极限，然而其它极限也是可能的。从矩阵K除去如下EM，所述EM的收益系数在上极限之上或者在下极限之下。如果EM的收益系数不能被确定，那么该收益系数保留在矩阵K中。优选地，也可以将在15页上描述的、按照本发明的方法“利用NNLS对所选择的宏观反应的反应速率的线性估计”用于与数据有关的预先减少。在那里找到该方法在所述预先减少的上下文中的解释。在该方法中对在与数据有关的预先减少中的过程数据的使用具有如下优点：所述减少与过程结合地并且借此更有效地以及更聚焦地来实现。

在步骤e) ii.中，利用如下算法从宏观反应中选择子集：针对所述选择，需要如下质量评定标准以及用于选择子集的算法，利用所述质量评定标准可以量化的是，利用子集（L）可以将来自d)的特定速率和/或来自c)的测量数据描绘得多好。

根据Soons等人，可以根据公式12将特定速率的平方剩余的和（SSR_q）用作利用子集L描绘所计算的特定速率的质量评定标准[Soons, Z. I. T. A., Ferreira, E.C., Rocha, I. (2010年). Selection of Elementary Modes for Bioprocess Control.11th International Symposium on Computer Applications in Biotechnology,Leuven, Belgium, 2010年7月7-9日, 第156-161页]：

（公式12），

SSR_q的值应该尽可能小。

针对SSR_q的最小化，事先必须针对每个被考虑的时间点t_i借助于非负最小二乘法来确定向量r(t_i)，使得适用：

（公式13），

利用附加的边界条件：

（公式14）。

该方法的优点是，也可以针对具有多个EM的很大的子集来执行根据公式12-14的计算。显著的缺点是，对于该计算来说必要的是所计算的特定速率。因为所述特定速率从经插值的测量值获得，所以所述特定速率以在它们的真实值方面的高的不可靠性来存在。测量精确度也许可以非常强烈地对所计算的特定速率有影响。因此，质量评定标准SSR_q也可以只在高的不可靠性的情况下被确定。除了关于所述描绘的质量的信息之外，利用该方法，也获得了子集L的反应速率r(t)的所估计的变化过程，作为根据公式13和14的最小化的结果。

Leighty, R.等人描述了另一方法，在所述另一方法中，测量值（浓度测量）直接通过对在时间上的容积反应速率的线性估计来逼近。在假定反应的变化过程在采样点之间线性地变化的情况下，通过利用线性最小二乘求解器来解决线性优化问题，可以快速地估计所述反应的变化过程 [Leighty, R. W., & Antoniewicz, M. R. (2011), Dynamicmetabolic flux analysis (DMFA): a framework for determining fluxes atmetabolic non-steady state. Metabolic engineering, 13(6), 第745-755页]。该方法仅仅涉及可逆的宏观反应（如在该来源中提到的“自由通量（Free Fluxes）”），此外，稀释效果（即不是由细胞造成的浓度变化）不能被考虑。如果宏观反应和测量值的大小不一致，那么这些测量值不能被用于估计反应速率。例如，如果以形成外部生物量为形式的细胞生长是宏观反应的部分并且根据所述测量值只公知细胞干燥量，那么这是这种情况。因而，在该形式下，该方法有关不可逆的宏观反应和补料分批培养过程的应用是不适合的。

如果在按照本发明被整理的（被变动的）数据的情况下使用Leighty等人的就此通过参考被集成到本申请中的设计，那么该方法现在也可以被应用到补料分批培养过程上。此外，还可以通过添加针对宏观反应的反应速率的下极限作为线性优化问题的边界条件，将该方法也用于不可逆的反应（如基本模式）。如果宏观反应和测量值的大小不一致，那么通过适当的相关可以使测量值的大小与宏观反应的大小适配。根据Leighty等人的线性估计与来自该权利要求中的改善方案的所述组合在下文被称作“对所选择的宏观反应的反应速率的线性估计”。

这样可以检查，利用原始EM集合K的子集L的所选择的宏观反应是否可以充分地呈现所测量的数据。在利用该方法确定的被变动的浓度与被变动的浓度之间的根据公式15的平方剩余SSR_C的这里所计算出的最终的和说明了利用该子集可以将测量数据描绘得多好：

（公式15），

SSR_C的值越小，该子集L就越好。该方法相对于Soons等人的方法，尤其是对于补料分批培养过程的建模来说是优选的，因为对子集的质量的检查在没有对特定速率的事先的可能有错误的确定的情况下也是可能的。通过假定所估计的反应速率在采样点之间线性地变化，测量偏差对反应速率的估计有很小的影响。该方法的缺点是，所要检查的子集L的大小通过线性优化问题的解决方案被限制。在该子集中的反应的最大数目等于存在的测量的数目除以采样点的数目。

除了关于所述描绘的质量的信息之外，利用该方法也获得了子集r(t)的反应速率的所估计的变化过程。

在按照本发明的“对所选择的宏观反应的反应速率的线性估计”的一个优选的实施方案中，为了解决线性优化问题，替代线性最小二乘求解器，使用Lawson等人的非负最小二乘求解器（NNLS）[Lawson, C.L. and R.J. Hanson, Solving Least SquaresProblems, Prentice-Hall, 1974, 第23章, 第161页]。由此，可能的是，利用该方法也检查更大的子集的质量。在这种情况下，宏观反应的最大数目也可以显著大于存在的测量的数目除以采样点的数目。“对所选择的宏观反应的反应速率的线性估计”与非负最小二乘求解器的使用的所述组合在下文被称作“利用NNLS对所选择的宏观反应的反应速率的线性估计”。

按照本发明的“利用NNLS对所选择的宏观反应的反应速率的线性估计”的方法附加地可以被用作用于在步骤e) i)中预先减少EM的另一与数据有关的方法。为此，这里可以使用宏观反应的很大的集合K。该方法的结果一方面是SSR_C的值而另一方面是反应速率r(t)的变化过程。从矩阵K除去具有所属的速率r(t)的小的值的EM。重复该行为方式，直至达到了EM的预先限定的数目或者SSR_C的值超过了所规定的极限值。

例如从Provost等人和Soons等人公知用于选择该子集的算法[Provost A. 2006.Metabolic design of dynamic bioreaction models. Faculté des Sciences Appliquées, Université catholique de Louvain, Louvain-la-Neuve, Louvain-la-Neuve;Soons, Z. I. T. A., Ferreira, E. C., & Rocha, I. (2010). Selection ofelementary modes for bioprocess control. 11th International Symposium onComputer Applications in Biotechnology, Leuven, Belgium, 2010年7月7-9日, 第156-161页]。

Soons等人描述了以二级优化方法形成EM子集。如上面所描述的那样，针对不同的、随机选择的EM分别使SSR_q的值最小化。选择具有最小的被最小化的SSR_q值的集合。然而，在EM的大的数目的情况下，随机的选择时无效的，因为可能的组合的数目非常强烈地增长。例如，在从20000个EM的集合中选择10个反应时存在超过2.8·10³⁶个组合。这里，找到最优的组合的概率非常小。由于使用质量评定标准SSR_q，该方法是易受测量不可靠性和测量偏差影响的。

Provost描述了一种可替换的算法，其中针对的不同的特定值确定基本模式的所有可能的正线性组合，其中适用：SSR_q(t=t_i) = 0。紧接着，从所述大量可能的组合中随机地选择一个组合。该方法分别只使用向量并且不使用在时间上的整个变化过程。因此，针对这个过程选择EM是不可能的。在随机的选择的情况下，虽然可以示出向量（就此可以将所述过程的剩余呈现到所述向量的程度），然而不能确定所述向量。该方法的另一缺点是，不能使用不在所有EM的解空间中的向量。近似解不能被确定。这尤其在具有不可靠性的测量和特定速率的情况下是大的缺点。由于使用质量评定标准SSR_q，该方法同样是易受测量不可靠性和测量偏差影响的。

因而，在该方法的一个其它并且优选的实施方式中，为了选择重要的宏观反应、也就是说为了选择EM子集L，使用进化、尤其是基因算法。这种算法例如从Baker等人公知[Syed Murtuza Baker, Kai Schallau, Björn H. Junker. 2010年的 Comparison ofdifferent algorithms for simultaneous estimation of multiple parameters inkinetic metabolic models. J. Integrative Bioinformatics:-1-1.]。特别优选地，可以使用如下基因算法，在所述基因算法的目标函数中，针对EM的不同的组合，利用所述方法“对所选择的宏观反应的反应速率的线性估计”计算相应的值SSR_C。可替换地，也可以使用随机的选择。在步骤ii)结束之后，矩阵L包含必需的宏观反应（步骤iii）。

在可选的步骤iii)中，以图形方式检查EM子集L的可靠性。这里，可以将来自步骤d)的通量图用作EM子集L的投影。图4示出了具有对六个EM的子集的投影的通量图。此外，如果EM子集L可靠的，那么测量数据继续处在EM子集L之内。该图示能够实现对所述选择的可靠性的快速的以图形方式的检查。

在另一步骤f)中，利用来自d)的特定速率和/或来自c)的测量数据，计算子集L的宏观反应的反应速率。根据Soons等人，对r(t)的计算可依据如在e)中所描述的那样的特定速率来进行[Soons, Z. I. T. A., Ferreira, E. C., Rocha, I. (2010).Selection of Elementary Modes for Bioprocess Control. 11th InternationalSymposium on Computer Applications in Biotechnology, Leuven, Belgium, 2010年7月7-9日, 第156-161页]，优选地，对r(t)的计算依据来自c)的测量数据利用按照本发明的“对所选择的宏观反应的反应速率的线性估计”来进行。

在该方法的步骤g)中，设计宏观反应的动力学。所确定的动力学应该使过程状态对相应的反应速率的动态影响量化：

（公式16），

从所述动力学得到所要确定的模型参数。

在步骤g) i.中根据宏观反应的化学计量来设计通用动力学。针对宏观反应的基质i，假定Monod类型的极限。最大反应速率乘以不同的极限：

（公式17），

其中Monod常数K_m,k,i和Hill系数n_i是方程的参数，所述方程的第一值手动地被输入。通常，Monod常数K_m,k,i被设定到相应的被测量的最大浓度的十分之一，而Hill系数n_i被设定到值1。根据反应化学计量对通用动力学的确定由Provost或者由Gao等人描述[Provost A.2006. Metabolic design of dynamic bioreaction models. Faculté des SciencesAppliquées, Université catholique de Louvain, Louvain-la-Neuve, 第126页;[Gao, J., Gorenflo, V. M., Scharer, J. M., & Budman, H. M. (2007). Dynamicmetabolic modeling for a MAB bioprocess. Biotechnology progress, 23(1), 第168-181页]。就此，这些方法通过参考而被集成到本申请中。在这些方法中，将Monod类型的基质极限用于反应的相应的基质。尽管Provost或Gao没有描述这一点，但是利用该方法根据反应化学计量同样能推导出由于有毒产品而引起的抑制。

在步骤g) ii.中，确定对在f)中所确定的反应速率r(t)的影响参量。在此，考虑如下所有参量，所述参量描述了过程状态（即也包括生物反应条件，诸如pH值、反应器温度、分压力，所述分压力不能从宏观反应的化学计量中推导出来）。影响参量可以例如借助于统计学方法（如偏最小二乘法（Partial-Least-Square））手动地来确定。为此，在过程状态（所述过程状态被合并在矩阵中）与来自f)的反应速率r(t)之间的相关被确定。

接着，在步骤g) iii.中，使在g) ii.中确定的影响量化并且使来自i.的动力学扩展相对应的项。接着，过程状态的参量对反应速率的在g) ii.中找到的影响可以配备项。所述项是任意的函数，所述函数根据过程状态而得到在0到1之间的值。接着，所述反应的在g) i.中列出的通用动力学乘以该项。

例如，所找到的在成分i的浓度与反应k之间的负相关暗示由于浓度i(C_i)引起的对反应速率k的影响。根据Haldane，这例如可以配备抑制动力学：

（公式18），

其中K_l,k,i表示抑制常数并且是另一模型参数，所述模型参数的第一个值手动地被输入并且通常被设定到相应的被测量的最大浓度的十分之一。

在可选的步骤h)中，动力学的模型参数值p与宏观反应的在f)中确定的反应速率r(t)适配：

（公式19），

这在下文被称作模型参数值估计。可以省去在该步骤中根据公式2至4对一个或多个微分方程的数值求解；模型参数值可以在通常具有3至10个参数的无关的组中单独地针对每个宏观反应k来适配。所述适配通过常见的方法（如高斯牛顿（Gauss-Newton）法）[BatesDM, Watts DG. 1988. Nonlinear regression analysis and its applications. NewYork: Wiley. xiv, 365.]来实现。

所述对于每个宏观反应来说单独的模型参数值估计对于步骤i)和j)来说是特别有利的，因为所述对于每个宏观反应来说单独的模型参数值估计一方面能快速地被执行，而另一方面在步骤j)中提供了用于使模型参数值与来自c)的测量数据适配的经改善的初始值。

适配质量例如根据公式20利用平方剩余SSR_r的和来计算：

（公式20），

SSR_r的值越小，所述适配越好。可替换地，通过与的以图形方式的比较来检查适配质量。

在可选的步骤i)中，有关宏观反应的在g)中所选择的动力学的适配质量来检查宏观反应的在g)中所选择的动力学。基础是在步骤h)中所计算的值SSR_r，所述值SSR_r使模型参数值估计的适配质量量化。在不令人满意的适配质量的情况下，可以重复步骤g)和h)，直至达到预先限定的适配质量。

在另一步骤j)中，根据对于适配来说常见的方法，可以使来自g)的动力学的参数值与来自c)的测量数据适配。优选地，将来自步骤h)的初始值用于所述适配。在包括上面提到的差分方程（公式2至4）在内的情况下，例如借助于高斯牛顿法[Bates DM, Watts DG.1988. Nonlinear regression analysis and its applications. New York: Wiley.xiv, 365.]或者利用多重打靶算法[Peifer M, Timmer J. 2007. Parameter estimationin ordinary differential equations for biochemical processes using the methodof multiple shooting. Systems Biology, IET 1(2):78–88]来实现模型参数值适配。

优选地，产品特征可以被集成到该模型中。特别优选地，这可以被用于如下产品特征，所述产品特征与副产品或中间产品的浓度有关。副产品（所述副产品是在a)中输入的新陈代谢网的外部成分）的浓度已经集成到该模型中并且可以被计算。不过，在需要时，也可以将其它副产品或中间产品合并在一个或多个单独的新陈代谢网中。如果所期望的分泌或摄取速率在不同的量级下波动或者应该考虑在不同的详细度下的确定的新陈代谢过程，那么这是有利的。替换于集成的模型地，利用步骤a)至j)可以生成一个单独的用于计算产品特性的模型，所述模型同样利用具有自己的动力学的宏观反应的原理来描述单独的新陈代谢网的外部成分的过程变化。副产品或中间产品（所述副产品或中间产品不在生物体之外，但是它们的细胞内部的积累对一个或多个产品特征有影响）可以在步骤a)和b)中在计算EM以及表达宏观反应的情况下被外部化、即被归类为外部成分。接着，与细胞内部的或者细胞外部的浓度有关的产品特征的结合可以通过在浓度与产品特征之间的定量或定性的相互关系的附加的集成来实现。

其它的主题是一种用于执行按照本发明的方法的计算机程序或者软件。

利用按照本发明的方法提供的模型可以被用于过程控制或者过程控制的计划以及对在反应器中的方法的检查。

按照本发明的方法的一个特别的实施方式示范性地被描述，而不对此进行限制。同样，利用该方法已经示范性地创建了利用杂交瘤细胞的发酵的模型并且如所描述的那样检查了所述模型的可靠性。

具体实施方式

例子：杂交瘤细胞培养的建模

1 步骤a)

以新陈代谢网的形式的背景知识已经从Niu等人的出版物(Metabolic pathwayanalysis and reduction for mammalian cell cultures—Towards macroscopicmodeling. Chemical Engineering Science (2013年) 102, 第461–473页, DOI:10.1016/j.ces.2013.07.034.)得知。动物细胞的这里所描述的新陈代谢网包含35个反应，所述35个反应使37个内部和外部代谢物相联系（参见图5，参见表1）。

表1：根据Niu等人的新陈代谢网的反应(Metabolic pathway analysis andreduction for mammalian cell cultures—Towards macroscopic modeling. ChemicalEngineering Science (2013年) 102, 第461–473页, DOI: 10.1016/j.ces.2013.07.034.)：

1 葡萄糖 → 1 G6P

1 G6P + 2 NAD → 2 丙酮酸

1 丙酮酸 → 1 乳酸 + 1 NAD

1 丙酮酸 → 1 丙酮酸_m

1 NADm + 1 丙酮酸_m → 1 乙酰辅酶A_m

1 乙酰辅酶A_m + 1 NADm + 1 草酰乙酸_m → 1 α-酮戊二酸_m

1 -酮戊二酸_m + 1 NADm → 1 琥珀酰辅酶A_m

1 FADm + 1 琥珀酰辅酶A_m → 1 延胡索酸

1 延胡索酸 → 1 苹果酸_m

1 苹果酸_m + 1 NADm → 1 草酰乙酸_m

1 谷氨酰胺 → 1 谷氨酸 + 1 NH3

1 谷氨酸 + 1 NADm →1 α-酮戊二酸_m + 1 NH3

1 苹果酸_m → 1 苹果酸

1 苹果酸 + 1 NAD → 1 丙酮酸

1 谷氨酸 + 1 丙酮酸 → 1 α-酮戊二酸_m + 1 丙氨酸

1 谷氨酸 + 1 草酰乙酸_m → 1 α-酮戊二酸_m + 1 天冬氨酸

1 精氨酸 + 2 NADm → 1 谷氨酸 + 3 NH3

1 天冬酰胺 → 1 天冬氨酸 + 1 NH3

2 甘氨酸 + 1 NADm → 1 NH3

1 组氨酸 + 1 NADm →1 谷氨酸 + 2 NH3

1 异亮氨酸 + 2 NADm → 1 乙酰辅酶A_m + 1 NH3 + 1 琥珀酰辅酶A_m

1 亮氨酸 + 3 NADm → 3 乙酰辅酶A_m

1 赖氨酸 + 6 NADm → 2 乙酰辅酶A_m

1 蛋氨酸 + 4 NADm → 1 NH3 + 1 琥珀酰辅酶A_m

1 NADm + 1 苯丙氨酸 → 1 酪氨酸

1 丝氨酸 → 1 NH3 + 1 丙酮酸

1 NADm + 1 苏氨酸→1 NH3 + 1 琥珀酰辅酶A_m

19 NADm + 1 TRP → 3 乙酰辅酶A_m

5 NADm + 1 酪氨酸 →2 乙酰辅酶A_m + 1 延胡索酸

5 NADm + 1 缬氨酸 →1 NH3 + 1 琥珀酰辅酶A_m

1 NADm →1 NAD

0.5 氧气(O2) →1 NADm

1 NADm →1 FADm

0.0156 丙氨酸 + 0.0082 精氨酸 + 0.0287 天冬氨酸 + 0.0167 G6P + 0.0245 谷氨酰胺 + 0.0039 谷氨酸 + 0.0196 甘氨酸 + 0.0038 组氨酸 + 0.0099 异亮氨酸 +0.0156 亮氨酸 + 0.0119 赖氨酸 + 0.0039 蛋氨酸 + 0.0065 苯丙氨酸 + 0.016 丝氨酸 + 0.0094 苏氨酸 + 0.0047 酪氨酸 + 0.0113 缬氨酸 →1 X (生物量) + 0.0981NAD

0.01101 丙氨酸 + 0.005033 精氨酸 + 0.007235 天冬酰胺 + 0.0081787 天冬氨酸+ 0.010381 谷氨酰胺 + 0.010695 谷氨酸 + 0.01447 甘氨酸 + 0.0034602 组氨酸 +0.005033 异亮氨酸 + 0.014155 亮氨酸 + 0.01447 赖氨酸 + 0.0028311 蛋氨酸 +0.007235 苯丙氨酸 + 0.026738 丝氨酸 + 0.016043 苏氨酸 + 0.0084932 酪氨酸 +0.018874 缬氨酸 →1 IgG (抗体) 。

在所述出版物中，没有明确地说明反应的可逆性。作为替代，用于“新陈代谢流分析”的数据已经根据同一出版物被分析并且被用于标识不可逆的反应。

利用化学计量矩阵N（所述化学计量矩阵N包含内部代谢物的化学计量、也就是说化学计量系数）以及关于反应的可逆性的信息，已经借助于METATOOL 5.1(Pfeiffer etal. METATOOL: for studying metabolic networks, Bioinformatics 199915 (3), 第251–257页)来计算该网络的所有基本模式（EM）。这里，EM的数目为超过300000。

2 步骤b)

具有所计算的EM E的矩阵已经在步骤a)中获得。类似于矩阵N地，矩阵N_p包含外部代谢物的化学计量、也就是说化学计量系数。化学计量网的可能的宏观反应已经利用公式21合并在矩阵K中：

（公式21）。

3 步骤c)

过程的测量数据已经从Baughman等人得知，所述过程说明了杂交瘤细胞在分批培养过程（参见图6）的变化过程期间的发酵的不同的测量参量[On the dynamic modeling ofmammalian cell metabolism and mAb production. In: Computers & ChemicalEngineering (2010) 34 (2), 第210–222页]。所述测量数据已经被输入到该方法中。

4 步骤d)

借助于来自c)的经样条插值的测量值（C ^int），已经计算出生长和衰亡速率以及特定的分泌和摄取速率（参见图7）。细胞溶解是包含预先限定并且被输入到该方法中的、在过程时长内恒定的细胞溶解因数K_l=0.1的。测量数据的变动已经不是必需的，因为这里涉及在没有其它添加物的情况下的分批培养过程的数据。与此相应地，数据显示出升高的变化过程，因为所有浓度变化都是由细胞造成，而不是由添加物造成。

为了计算速率而使用附加信息。这样，借助于同样在数据组中由Baughman等人说明的总生物量（BM以为单位）以及总细胞数目已经可以计算出平均C-mol含量为。关于C-mol的生长速率现在已经可以利用公式22：

（公式22）

来计算。类似地，可以估计抗体的关于C-mol的形成速率。为此，抗体的摩尔成分被估计成具有形式上的摩尔质量为的。这里假定：所述摩尔成分对应于如由Villadsen等人[Bioreaction engineering principles (2011),Chapter 3, Elemental and Redox Balances, 第73页, Springer Verlag, ISBN: 978-1-4419-9687-9]说明的那样的蛋白质的平均摩尔成分。整个抗体的摩尔质量已经被估计成。接着，根据如下公式得到该抗体的形成速率：

（公式23）

接着，所述速率的时间变化过程可以被用于选择宏观反应。

5 步骤e)

在步骤e)中，已经生成了宏观反应的EM子集，利用所述EM子集已经尽可能好地重现数据组。为此是需要来自步骤b)的矩阵K。因为超过300000个宏观反应的数目已经导致可能的组合的过大的数目，所以首先执行了与数据有关的预先减少：

为此，已经将在步骤d)中确定的速率用于计算两个外部代谢物的所有组合的收益系数。收益系数的下极限被选择为使得所确定的收益系数的99%都超过该值。上极限已经被选择为使得所确定的收益系数的99%都低于该值。示范性地，在表2中说明了一些所确定的极限以及EM的份额，所述EM的份额的收益系数在所述极限之内。整体上，这样可以将EM的数目减少到大约3000。

表2：外部代谢物，所述外部代谢物的最大和最小收益系数以及EM的份额，所述EM的份额的收益系数在所说明的极限之内。

外部成分	下极限	上极限	在极限之内的EM的份额
				Ala：Asn	-14.6811	-0.1752	99.9134%
Asn：Glc	0.0293	16.1603	64.0488%
				Asp：Ala	-1.1130	-2.1011	74.0839%

紧接着，紧接着与数据有关的减少，还执行了与数据无关的减少。在这种情况下，限定了两个EM的余弦相似性的最大值为0.995。这样，从第一反应开始，从矩阵K中除去了超过所述值的所有宏观反应。从矩阵K（也称作被减少的矩阵K）中保留有大约500个宏观反应，所述大约500个宏观反应继续覆盖通过大约3000个EM展开的解空间的超过95%的容积。

此外，在选择过程之前，实现在根据Niu等人的新陈代谢网中说明的成分（所述成分对应于来自a）的新陈代谢网的外部代谢物）与来自c)的所确定的成分浓度的匹配。除了脯氨酸之外，所有由Baughman等人测量的浓度也在根据Niu等人的新陈代谢网中被考虑。为了能够使用对脯氨酸浓度的测量，或者能使用另一经简化的网络（所述网络包含脯氨酸作为外部代谢物），或者对现有的新陈代谢网的扩展是可能的。

在下文，同样忽略了如下成分，所述成分虽然出现在所计算的宏观反应中，但是所述成分中没有数据能支配。与此相应地，矩阵K的相对应的行是从该矩阵被划去的。划去相对应的行并不意味着不使用所述行的输入和输出。所述成分继续存在于新陈代谢网中，只是缺少了为此的能使所述成分能被调和的测量。在该例子中，忽略了如下成分的输入和输出：精氨酸、谷氨酸、甘氨酸、组氨酸、亮氨酸、赖氨酸、蛋氨酸、铵、氧气、苯丙氨酸、丝氨酸、苏氨酸、色氨酸、酪氨酸和缬氨酸。

现在，在该方法的下面的步骤中，使用被减少的矩阵K（所述矩阵K描绘了背景知识）和来自d)的速率以及来自c)的测量数据（所述测量数据形成过程知识），用来获得来自K的宏观反应的尽可能小的子集L。

已经将按照本发明的“对所选择的宏观反应的反应速率的线性估计”用于质量评定标准。

类似于速率地，这里已经使细胞数目和抗体的测量值标准化到C-mol。这是必要的，以便宏观反应的大小与测量的大小一致。

利用基因算法来选择子集。在这种情况下，在对该基因算法的目标函数的计算中，已经解决了在“对所选择的宏观反应的反应速率的线性估计”中谈及的线性优化问题。所述线性优化问题的最小误差平方的这里所计算出的最终的和同时是所述宏观反应的相应的选择的目标函数的值。

为了选择来自K的子集L的大小，已经以在L中的不同数目的宏观反应重复地执行了所述优化，该数目是模型复杂性与重现精确度之间的妥协。为了确定多少反应对于重现来说足够，或者可以针对变化的数目的宏观反应来重复子集L的选择，或者反应数目的惩罚项将直接被添加给基因算法的目标函数。在这种情况下，已经利用预先限定的数目（10、7、5、4和3）的宏观反应来执行多个优化。在图9中，误差平方的以基因算法找到的最小和相对于宏观反应的数目来绘制。已经显示出：在这种情况下，为了足够好地呈现过程变化，少于七个宏观反应太少。在表3中说明了所选择的宏观反应。

表3：宏观反应的所选择的子集（L）。没有下划线的成分在该模型中不被考虑，因为为此没有进行测量：

0.474 丙氨酸 + 0.474 蛋氨酸

→ 0.158 天冬酰胺 + 0.316 天冬氨酸 + 0.632 甘氨酸

+ 0.158 色氨酸

0.015 丙氨酸 + 0.00789 精氨酸 + 0.0304 天冬酰胺 + 0.0161 葡萄糖

+ 0.0236 谷氨酰胺 + 0.00375 谷氨酸 + 0.00366 组氨酸

+ 0.00953 异亮氨酸 + 0.015 亮氨酸 + 0.112 蛋氨酸

+ 0.00626 苯丙氨酸 + 0.0154 丝氨酸 + 0.0109 缬氨酸

→ 0.963 X (生物量) + 0.00276 天冬氨酸 + 0.24 甘氨酸

+ 0.0208 色氨酸

0.295 天冬酰胺 + 0.147 谷氨酸

→ 0.295 天冬氨酸 + 0.885 甘氨酸 + 0.147 乳酸

0.00753 精氨酸 + 0.113 天冬酰胺 + 0.0603 葡萄糖 + 0.0225 谷氨酰胺

+ 0.0824 组氨酸 + 0.00909 异亮氨酸 + 0.00597 苯丙氨酸

+ 0.0216 色氨酸+ 0.00431 酪氨酸 + 0.0104 缬氨酸

→ 0.918 X (生物量) + 0.061 丙氨酸 + 0.0865 天冬氨酸

+ 0.343 甘氨酸 + 0.0631 蛋氨酸

0.0654 精氨酸 + 0.412 天冬氨酸 + 0.00991 葡萄糖 + 0.0145 谷氨酰胺

+ 0.554 甘氨酸 + 0.00226 组氨酸 + 0.00588 异亮氨酸

+ 0.00926 亮氨酸 + 0.00706 赖氨酸 + 0.0649 苯丙氨酸

+ 0.0095 丝氨酸 + 0.00671 缬氨酸

→ 0.594 X (生物量) + 0.049 丙氨酸 + 0.395 天冬酰胺

+ 0.0503 苏氨酸 + 0.0388 色氨酸

0.0077 精氨酸 + 0.179 天冬氨酸 + 0.0157 葡萄糖 + 0.104 谷氨酰胺

+ 0.216 甘氨酸 + 0.00357 组氨酸 + 0.00929 异亮氨酸

+ 0.0146 亮氨酸 + 0.0112 赖氨酸 + 0.038 酪氨酸 + 0.0106 缬氨酸

→ 0.939 X (生物量) + 0.0624 丙氨酸 + 0.152 天冬酰胺

+ 0.0183 色氨酸

0.0342 精氨酸 + 0.211 天冬氨酸 + 0.00762 葡萄糖 + 0.0195 谷氨酰胺

+ 0.244 甘氨酸 + 0.00452 组氨酸 + 0.0546 异亮氨酸

+ 0.0185 亮氨酸 + 0.0171 赖氨酸 + 0.00406 蛋氨酸

+ 0.0178 酪氨酸 + 0.0203 缬氨酸

→ 0.457 X (生物量) + 0.804 IgG (抗体) + 0.185 天冬酰胺

+ 0.0153 色氨酸。

在所示出的宏观反应中，说明了来自a)的新陈代谢网的所有外部代谢物。然而，所述模型的部分只是有下划线的外部代谢物，因为只有针对这些外部代谢物存在来自c)的测量数据。

6 步骤f)

对于宏观反应的所选择的集合，已经确定了在时间上的反应速率。在该例子中，已经利用按照本发明的方法“对所选择的宏观反应的反应速率的线性估计”，通过估计反应速率r(t)来逼近在图10中示出的测量值。该方法的结果是各个（容积）反应速率的分段地线性的变化过程。通过除以活细胞数目的经插值的变化过程X_v(t)，获得了在表3中示出的宏观反应的细胞特定的反应速率r(t)。这样获得的反应速率r(t)在图10中描绘。

7 步骤g)

针对所有在表3中示出的宏观反应，已经假定按照公式24的通用动力学：

（公式24）。

在这种情况下，所述宏观反应已经通过Monod动力学来实现。也就是说，在每个反应k中针对每个基质i按照公式25：

（公式25）

来引入极限。在这种情况下，是最大反应速度，N_l是所考虑的极限的数目，C_i是成分i的浓度，K_m,k,i是所属的Monod常数而n_i是反应顺序的Hill参数。它们的值在步骤h)和j)中被适配。

从对来自f)的反应速率r(t)的分析中得到其它项。在该例子中，除了基质极限之外也考虑按照公式26的抑制：

（公式26）

针对该极限，也必须适配参数K_l,k,i和n_i的值。所述反应的所使用的动力学项在表4中说明。

表4：来自L的所选择宏观反应的动力学项

。

8 步骤h)

对于每个反应速率来说，利用在表4中说明的动力学和在该动力学下被考虑的浓度C ^int(t)的经插值的值已经可以以代数方式来计算反应速率的变化过程。

对于每个反应i来说，使所述动力学的参数与在步骤f)中确定的反应速率r_i(t)单独地适配。在该例子中，针对在反应i中出现的参数的优化的目标函数是：

（公式27）

所有被计算的的这样适配的变化过程与相对应的r_k(t)一起在图11中示出。前者的变化过程以虚线示出，后者的变化过程以实线示出。能识别出的是，所述变化过程在定性上一致。这意味着：利用所选择的动力学也可以令人满意地重现该过程的动力学。在所述建模步骤中，该信息是非常有帮助的，因为在不令人满意的重现的情况下，可以重复要快速执行的步骤g)（选择其它动力学）和h)（参数值估计），直至达到了所希望的适配度。因此，步骤i)这里是必要的。

9 步骤j)

利用来自c)的测量数据实现对模型参数值p的继续适配。为此，所有参数已经同时被优化。此外，到目前位置没有被考虑的过程包括编程性细胞死亡和细胞溶解在内。这些在微分方程中是必需的，所述微分方程描述了活细胞数目和总细胞数目的增长：

（公式28）

（公式29）

针对编程性细胞死亡的所选择的动力学是：

（公式30）

（公式31）

细胞溶解速率K_l已经被假定为在过程期间恒定。在该步骤中，除了反应速率的参数之外，还确定了由于编程性细胞死亡和细胞溶解而引入的参数（临界乳酸浓度）、（最大衰亡速率）、（用于描述乳酸浓度对衰亡速率的影响的Monod参数）以及K_l（细胞溶解速率）。在该例子中，从数据组的初始值出发，已经通过ODE系统的数值求解来确定所估计的浓度的变化过程。在此，在所测量的浓度C ^m(t)与所估计的浓度之间的差已经以常见的方法利用如下目标函数来最小化：

（公式32）

一般来说，以总共33个参数p要困难地执行所述优化，因为目标函数具有多个局部最优值。而如果在所述参数的从步骤h)已知的初始值方面开始确定性优化算法（诸如列文伯格-马夸尔特（Levenberg-Marquardt）法），那么强烈地提高了成功前景。在图12中描绘了所适配的过程变化。在表5中描绘了所适配的参数。

动力学以及编程性细胞死亡和细胞溶解的参数

K_m,Glc,1	14.6	K_m,Gln,7	0.0187
				K_m,Ala,1	3.41	K_m,Asp,7	0.872
K_m,Glc,2	0.0508	r_1,max	9.47
				K_m,Gln,2	0.00881	r_2,max	9.91
K_m,Asn,2	1.38	r_3,max	57.6
				K_m,Ala,2	2.19	r_4,max	21.7
K_m,Glc,3	7.13	r_5,max	0.345
				K_m,Asn,3	6.84	r_6,max	49.4
K_m,Xt,4	0.0315	r_7,max	3.03
				K_m,Glc,4	1.29	K_IAsn,1	16.1
K_m,Gln,4	2.19	K_I,Lac,3	0.681
				K_m,Asn,4	1.68	K_I,Asp,5	9.74
K_m,Glc,5	100	K_I,Asn,5	1.10
				K_m,Gln,5	28.2	μ_d,max	0.125
K_m,Asp,5	102	K_d,Lac	1.01
				K_m,Glc,6	0.0451	C_Lac,cr	1.22
K_m,Gln,6	0.791	K_l	0.00843
				K_m,Asp,6	1.06	K_m,Glc,7	0.0145

10 步骤k)

已经利用来自表5的所属的参数值输出了所述模型，所述模型有矩阵L、来自表4的动力学以及编程性细胞死亡的动力学组成。

符号列表

_（下划线）	表示向量
		_i（索引i）	表示向量的第i个元素
_k（索引k）	表示向量的第k个元素
		[ ]	表示在括号中的成分的浓度
C	浓度
		ΔC	浓度差
C^Int	经插值的浓度
			所估计的浓度（例如通过对微分方程求解）
C_s	被变动的浓度
		C_cr	临界浓度
C^m	所测量的浓度
		D	稀释速率
q	所确定的细胞特定的分泌和摄取速率
			所确定的细胞特定的分泌和摄取速率，所述分泌和摄取速率已经从任意的单位换算到
r	所确定的反应速率
			所估计的反应速率（例如通过计算反应动力学）
	动力学的极限
		r_max	反应动力学的参数
N	化学计量矩阵
		N_p	外部化学计量矩阵
K	包含宏观反应的矩阵
		E	包含所有基本模式的矩阵
X_t	总细胞数目
		X_v	活细胞数目
	生长速率
			衰亡速率
	生长速率，所述生长速率已经从任意的单位被换算到
		K_d	细胞溶解速率
K_I	抑制极限的参数
		K_M	基质极限的参数
n	抑制极限或基质极限的Hill参数
		L	宏观反应的被用于该模型的子集
p	模型参数
		S	基质
SSR_q	特定的分泌或摄取速率的平方剩余的和
		SSR_C	浓度的平方剩余的和
SSR_r	反应速率的平方剩余的和

附图说明

对附图的描述：

图1示出了测量数据的变动：所示出的是所测量的参量(C_i(t))的实际变化过程，所述实际变化过程在稀释速率(D(t))改变时突然改变。被变动的变化过程(C_i,s(t))仅仅由于由细胞造成的改变而出现。

图2示出了两个特定速率q₁和q₂的通量图。等高线说明了如下频率，所述速率的相应的组合在所测量的数据中以所述频率出现。

图3示出了解空间的三维图示，所述解空间通过EM的正线性组合来展开。用黑色示出了整个集合的解空间，用灰色示出了子集的解空间。

图4示出了两个特定速率q₁和q₂的通量图。示出了集合L的宏观反应的2维投影，作为向量。

图5示出了根据Niu等人的新陈代谢网的示意图。在这种情况下，将细胞的边界示出为箱子。线粒体的细胞内部的边界用虚线来表征。外部成分用索引“xt”来表征。箭头和点状箭头表征反应。

图6示出了根据Baughman等人的利用杂交瘤细胞的发酵的测量数据。在这种情况下，总细胞数目（total Cells）根据活细胞（vital Cells）和死细胞（dead Cells）的和来计算。缩写GLC、GLN、ASP、ASN、LAC、ALA和PRO表示基质葡萄糖和氨基酸谷氨酰胺、天冬氨酸、天冬酰胺、丙氨酸和脯氨酸以及新陈代谢产品乳酸。缩写MAB表示单细胞繁殖的抗体的产品而BM表示生物量。

图7示出了生长和衰亡速率以及细胞特定的分泌和摄取速率。除了之外的所有细胞特定的速率都以为单位来说明。速率以来说明。

图8示出了以利用所选择的反应集合“对所选择的宏观反应的反应速率的线性估计”来逼近的浓度。为此，总细胞数目(X_t)以及抗体浓度(MAB)已经被换算到C-mol。

图9示出了在子集中的宏观反应的数目(n_R)上绘制的误差平方的所确定的最小和（“Minimum error（最小误差）”）。

图10示出了宏观反应的利用按照本发明的方法“对所选择的宏观反应的反应速率的线性估计”所确定的反应速率r(t)。

图11示出了宏观反应的利用按照本发明的方法“对所选择的宏观反应的反应速率的线性估计”所确定的反应速率r(t)（实线）连同以代数方式计算的反应速率（虚线）。

图12示出了所测量的浓度C ^m(t)（点）与所模拟的过程变化（实线）的比较。所述浓度以[mM]为单位说明。例外是活细胞数目和总细胞数目（以[10⁹个细胞/l]为单位的X_v/X_t）以及抗体的浓度（以[10^-4mM]为单位的mAb）。

Claims

1.一种用于创建利用生物体的生物反应的模型的计算机实现的方法，所述方法包括如下步骤：

a. 所述生物体的所选择的新陈代谢路径、所述新陈代谢路径的化学计量以及可逆性特性作为背景知识被输入到所述方法中，而且基本模式根据所述输入来计算；

b. 来自a)的所述基本模式在矩阵K中合并，其中所述基本模式将来自a)的新陈代谢路径合并成宏观反应，而且所述矩阵K包含所有宏观反应的化学计量和可逆性特性；

c. 用于利用所述生物体的生物反应的测量数据被输入；

d. 借助于插值法，基于来自c)的所输入的测量数据来计算所输入的新陈代谢路径的对于所述生物体来说特定的速率、即一个或多个输入参量和输出参量中的分泌和摄取速率；

e. 重要的宏观反应以来自a)的基本模式的子集的形式通过如下方式来选择：

i. 与数据无关地和/或与数据有关地预先减少来自a)的基本模式的数目，

ii. 借助于根据数学质量评定标准的算法和子集在矩阵L中的合并，利用来自c)的测量数据和/或来自d)的一个或多个速率，从来自c) i.的所述预先减少中选择所述子集，

iii. 可选地以图形方式呈现所述子集；

f. 借助于插值法，基于来自c)的所输入的测量数据和/或来自d)的速率来计算所述子集r(t)的宏观反应的反应速率；

i. 根据所述宏观反应的化学计量来设计通用动力学，

ii. 对所述宏观反应的影响参量根据来自f)的反应速率来确定，

h. 可选地，根据来自g)的动力学，使每个宏观反应的模型参数值单独地与每个宏观反应的来自f)的所计算的反应速率进行首次适配；

j. 所述模型参数值与来自c)的测量数据适配；

k. 矩阵L、来自g)的动力学以及来自j)的模型参数值形成所述模型并且被输出和/或被转用到过程控制或过程研发模块中。

2.根据权利要求1所述的计算机实现的方法，其中，在步骤d)中，也计算所述生物体的生长速率、特别优选地也计算所述生物体的衰亡速率。

3.根据权利要求1或2之一所述的计算机实现的方法，其中，在步骤g)中，基于对来自f)的反应速率的分析对所述动力学进行单独的适配。

4.根据权利要求1至3之一所述的计算机实现的方法，其中，在步骤h)中，通过多个适配方法的组合来进行对来自g)的动力学的参数值的适配。

5.根据权利要求1至4之一所述的计算机实现的方法，其中，在步骤e) ii.中，为了选择宏观反应的子集而执行对所选择的宏观反应的反应速率的线性估计。

6.根据权利要求1至5之一所述的计算机实现的方法，其中，在步骤e) ii.中，为了选择宏观反应的子集，与进化算法相结合地执行对所选择的宏观反应的反应速率的线性估计。

7.根据权利要求1至6之一所述的计算机实现的方法，其中，所述测量数据在步骤d)中应用插值法之前予以变动，以便实现对在没有饲料高峰的情况下的恒定的消耗的描述。

8.根据权利要求1至7之一所述的计算机实现的方法，其中，在步骤f)中，执行对所选择的宏观反应的反应速率的线性估计。

9.根据权利要求1至8之一所述的计算机实现的方法，其中，在步骤e) i.中，执行与数据有关的预先减少，而且针对所述与数据有关的预先减少使用利用NNLS对所选择的宏观反应的反应速率的线性估计的方法。

10.根据权利要求1至9之一所述的计算机实现的方法，其中，在步骤e) iii.中，借助于通量图来检查对宏观反应的子集的选择的可靠性。

11.根据权利要求1至10之一所述的计算机实现的方法，其中，在步骤e) ii.中，利用来自c)的测量数据，从来自e) i.的预先减少中进行选择。

12.一种计算机程序，用于执行根据权利要求1至10之一所述的方法步骤。

13.一种软件，用于执行根据权利要求1至11之一所述的方法步骤。