JP6816003B2 - 発酵モデルを作成するためのコンピューター実行式の方法 - Google Patents

Description

本発明は、オーガニズムを用いた、生体反応、特に発酵または全細胞触媒作用のモデルを作成するためのコンピューター実行式の方法に関する。

本明細書の文脈における「オーガニズム」とは、生体反応で用いられる、哺乳動物細胞、酵母菌類、細菌類、藻類といった、動物または植物細胞の培養物を指す。

発酵プロセスの知覚によるモニタリング、および、例えばＢａｙｅｒ Ｔｅｃｈｎｏｌｏｇｙ Ｓｅｒｖｉｃｅｓ ＧｍｂＨ社製のクオリティ・バイ・デザイン分析の自動化プラットフォーム、ＢａｙｃｈｒｏＭＡＴ（登録商標）を用いたプロセスから得たサンプルの分析からは、バイオリアクター内のプロセスの状態に関して様々なリアルタイムの情報が提供される。通常は、細胞数、細胞の生存、炭素源（例えばグルコース）、アミノ酸またはＯ_２、生成物および副産物（例えば、乳酸もしくはＣＯ_２）などの基質の濃度、温度および／もしくはｐＨなどのプロセスパラメーターまたは生成物の特徴が決定される。これらのデータは、例えば従来技術から得た、計算データおよび／または補間で補ってもよい。これらのデータを合わせて、本明細書の文脈での測定データまたはプロセスノウハウが形成される。

本明細書の文脈では、オーガニズムに対する背景知識は、特異的および非特異的反応も含めたオーガニズムの生化学的反応、および、特に細胞内での反応、すなわち、基質、代謝物（代謝ネットワークのノードとも呼ばれる）、生成物およびそれらの間で起こる生化学的反応からなるオーガニズム特有の代謝ネットワーク（ＭＮ）を表現するマクロ反応に関する知識を指す。これらの生化学的反応は、以下の項目によって定義される：
（ａ）化学量論
（ｂ）可逆性（生物学的条件下で）
（ｃ）化学量論的ネットワークへの統合。

現在まで、主に測定データがプロセスの定性的モニタリングのために使用されてきた。ここで、解決のために動的プロセスモデルを必要とする技術的な問題を選択して提示する。

本明細書の文脈でのプロセスノウハウの技術的な使用の１つにより、バイオリアクター内でのプロセス状態の、モデルベースの推定がもたらされる。「拡張カルマンフィルタ」といった方法によって、プロセスの量の連続的な推定が可能になるが、これに関して、非連続的な方法での測定がある［Welch G, Bishop G. 1995. An Introduction to the Kalman Filter. Chapel Hill, NC, USA: University of North Carolina at Chapel Hill.］。測定不能な量の経過も、他の測定から計算できる。これには、その根底となるプロセスを正しく記述するプロセスモデルが必要である。

さらなる使用は、モデルベースの最適プロセス制御である。最適プロセス制御では、モデルベースの予測閉ループ制御回路における生成物の量、生成物の特徴もしくは副産物の形成または他のターゲットとなる量に対してプロセス制御を最適化するために、動的プロセスモデルを使用する。例えば、Frahm et al.がハイブリドーマ細胞培養に関して実証している［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 for parameter precision: State of the art. Chemical Engineering Science 63(19):4846-4872］。モデルの複雑性に関する上記必要条件の他に、前記技術的な使用では、動的プロセスモデルが開発段階の際にすでに存在している必要がある。プロセス開発にかかる時間を最小限に抑えるために、既存のプロセスノウハウからできるだけ早く作成が可能でなくてはいけない。

したがって、背景知識とプロセスノウハウを使用して動的プロセスモデルを作成することが可能な方法を提供する必要性があった。このモデルを、例えば、状態推定、最適プロセス制御またはモデルベースの実験計画との関連で使用できるようにするには、モデルの複雑性が低い必要がある。依存性、すなわちプロセス量またはプロセス状態のプロセス挙動に対する影響は、デザインスペース内で十分正確に定量化されるべきである。プロセス状態に関する利用可能なすべての情報がこの目的のために使用されなくてはいけない。必要に応じて、生成物の特徴のモデルベースの記述をモデルに統合することが可能でなくてはいけない。プロセスノウハウが存在する領域をデザインスペースと呼ぶ。この方法は、上述の生体反応に適用可能であり、そのような動的モデルの開発時間を実質的に短縮するものでなければならない。動的モデル開発への現在までのアプローチは、プロセスモデルの完成まで数カ月から数年かかる。経験上、本アプローチは開発時間を数週間まで短縮する。

本明細書の文脈における典型的な生成物の特徴は、例えば、タンパク質のグリコシル化パターンまたはタンパク質の完全性が挙げられるが、これらに限定されない。現在まで、上記の文脈で使用された動的モデルはこの特性を有していない。本アプローチは、生成物の特徴の単純なモデルベースの統合を可能にする。

発酵のモデルベースのプロセス制御は、例えば、ハイブリドーマ細胞培養に関してFrahm et al.によって実証されている（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）。ここで、根本的なプロセス量の制御はモデルベースである。生成物の特徴の統合はない。この細胞の数学的モデルは、この特定のプロセスのために考案されたものであり、同じオーガニズムまたは他のオーガニズムまたは同じオーガニズムの株を用いたプロセスにはかなりの労力を要さないと移行できない。細胞内の反応形式での背景知識は、モデルでは明示的に考慮されていない。追加の測定量のモデルへの統合、よって、プロセス状態に関する情報の完全な利用は、この場合、かなりの労力を要さないと可能にならない。したがって、このアプローチは、他のプロセスに移行できず、取得したデータを完全に利用することができない個別の解決策となる。上記の方法は、予測されるモデル開発の時間と、同じオーガニズムまたは他のオーガニズムでの他のプロセスへの解決策の移行の労力を必要性のために、上記の技術的問題を解決するものではない。

グリコシル化などの生成物の特徴も組み込んだより広範なモデル化が、Kontoravdi et alの研究で見つけることができる。主要な代謝を記述するこのモデルは、細胞内の反応形式の背景知識を組み入れず、同じオーガニズムまたは他のオーガニズムのプロセスにも移行できない。この場合、追加の測定量をモデルに統合することはできない［Kontoravdi C, Asprey SP, Pistikopoulos EN, Mantalaris A. 2007. Development of a dynamic model of monoclonal antibody production and glycosylation for product quality monitoring. Computers & Chemical Engineering 31(5-6):392-400.］。この方法もまた、プロセス状態に関する情報の完全な利用を可能にせず、長いモデル開発時間を必要とし、他のオーガニズムまたは株に移行できない。したがって、この方法は技術的問題の解決策にはならない。

Jedrzejewski et al.およびJimenez et al.のヌクレオチド糖代謝の統合を伴うグリコシル化のモデルは、内部代謝中間体の平衡方程式の形式の背景知識を組み込む［Jedrzejewski PM, del Val, Ioscani Jimenez, Constantinou A, Dell A, Haslam SM, Polizzi KM, Kontoravdi C. 2014. Towards Controlling the Glycoform: A Model Framework Linking Extracellular Metabolites to Antibody Glycosylation. International journal of molecular sciences 15(3):4492-4522.; Jimenez del Val, Ioscani, Nagy JM, Kontoravdi C. 2011. A dynamic mathematical model for monoclonal antibody N-linked glycosylation and nucleotide sugar donor transport within a maturing Golgi apparatus. Biotechnology progress 27(6):1730-1743］。しかし、プロセス制御にこのモデルを使用する場合、モデル全体の複雑性と細胞内の代謝中間体の不十分な観察度が欠点となる。さらに、主要な代謝のモデルでは、他のプロセスへの移行またはプロセス状態に関する情報の完全な利用が可能にならない。したがって、この方法は技術的問題の解決策にはならない。

バイオプロセスのためのモデルの柔軟な作成が、Leifheit et al.によって提案されている［Leifheit J, Heine T, Kawohl M, King R. 2007. Rechnergestutzte halbautomatische Modellierung biotechnologischer Prozesse [Computer-aided semi-automatic modelling of biotechnology processes]. At - Automatisierungstechnik 55(5)］。このモデルは、背景知識なしで、プロセスノウハウをもとに作成される。この手順は、同じオーガニズムまたは他のオーガニズムでの様々なプロセスに使用することができる。使用者本人が指定するマクロ反応に基づいている。マクロ反応の正確な化学量論がこの方法で決定される。この方法は、少数の状態量または測定量について記述される。追加の状態量または測定量の統合は、方法の複雑性を大幅に増加する。例えばＢａｙｃｈｒｏＭＡＴ（登録商標）プラットフォームによって提供されるような包括的なデータの原則が使用される場合、この方法は実現可能ではなくなる。この方法では、生成物の特徴を統合することはできない。したがって、上記技術的問題の解決策にはならない。

オーガニズムの既知の代謝（化学量論的）ネットワークからエレメンタリーモード（ＥＭ）として得られるマクロ反応の形式の背景知識の使用を、Provostが説明している［Provost A. 2006. Metabolic design of dynamic bioreaction models. Faculte des Sciences Appliquees, Universite catholique de Louvain, Louvain-la-Neuve, Louvain-la-Neuve, p. 81 ff., p. 107 ff. p. 118 ff.］。この方法は、様々なオーガニズムまたは同じオーガニズムの株に使用することができる。プロセスモデルに関するマクロ反応は、プロセスノウハウを使用して選択される。しかし、所定の数のマクロ反応がそれぞれランダムに選択され、プロセスセグメントが定義される。この場合、開示される方法は、エレメンタリーモードの多くの可能な組み合わせのうちの１つを提供する。マクロ反応の数、したがってモデルの複雑性は固定されており、変更できない。この方法は、各プロセスセグメントに対して別個のモデルを作成する。個々のマクロ反応の動態の選択は、選択されたマクロ反応の化学量論を考慮に入れる。しかし、動力学パラメーター（モデルパラメーター）はプロセスデータに対して調整されない。代わりに、プロセスセグメントモデルを別々に使用して、プロセスデータが変更される。反応のランダムな選択は包括的なデータベースをもとに行うことができるが、動態の選択および選択された動態に関して記述されるこのアプローチは、プロセスの経過またはプロセスにおけるオーガニズムの挙動を表現することはできない。さらに、複数のプロセスセグメントモデルの使用により、プロセスモデルの複雑性に不必要な増加がもたらされる。依存性、すなわちプロセス量またはプロセス状態の対するプロセス挙動に対する影響は、この方法では定量化できない。再度述べるが、生成物の特徴の統合はない。よって、この方法は上記技術的問題の解決策にはならない。

したがって、プロセスノウハウおよび測定データに基づくモデルの迅速かつ効率的な提供、ならびに、背景知識を考慮した生成物のターンオーバーおよび重要な生成物の特徴の最適化を可能にし、上記の欠点を有さない方法を提供する必要性があった。

この目的は、バイオリアクター内でのオーガニズムとの生体反応のモデルを作成する方法によって達成され、以下のとおりである。

本明細書は、オーガニズムとの生体反応、特に発酵または全細胞触媒作用のモデルを作成するためのコンピューター実行式の方法を提供し、以下の工程を含んでなる。

ａ．オーガニズムの選択された複数の代謝経路と、これらの代謝経路の化学量論的特性および可逆特性とを背景知識として方法に組み込む工程。すなわち、オーガニズムの１つまたは複数の代謝ネットワークを方法に組み込む。複数のエレメンタリーモード（ＥＭ）はこの入力によって計算される。

ｂ．ＥＭを行列Ｋに結合し、ＥＭは、ａ）からの代謝経路をマクロ反応に結合する工程。したがって、この行列Ｋは、背景知識から得たすべての可能なマクロ反応の化学量論的特性および可逆特性を含む。

ｃ．オーガニズムとの生体反応に関する測定データ（プロセスノウハウとも呼ばれる）を入力する工程。

ｄ．補間法を用いて、オーガニズムに関する入力した代謝経路の特異速度（１つまたは複数の入力量および出力量の分泌および摂取速度）を、ｃ）で入力した測定データに基づいて計算する工程。好ましくは、オーガニズムの増殖率、特に好ましくは死亡率も計算する。

ｅ．ｂ）のエレメンタリーモードのサブセットの形式で、関連するマクロ反応を、
ｉ．ｂ）のＥＭの数のデータ非依存性および／またはデータ依存性の事前縮小と、
ｉｉ．ｃ）の測定データおよび／またはｄ）の１つまたは複数の速度、好ましくはｃ）の測定データを使用して、数学的品質基準に従ったアルゴリズムによりｅ）ｉ．で事前縮小したサブセットを選択し、サブセットを行列Ｌに結合することと、
ｉｉｉ．任意でサブセットを図で示すことと、
により選択する工程。

ｆ．補間法を使用して、サブセットのマクロ反応の反応速度
を、ｃ）の入力測定データおよび／またはｄ）の速度に基づいて計算する工程。

ｇ．ｅ）ｉｉ．のサブセットのマクロ反応の動態を、以下の中間工程、
ｉ．一般的な動態をマクロ反応の化学量論から考案する工程と、
ｉｉ．マクロ反応に影響を及ぼす量をｆ）の反応速度から決定する工程と、
ｉｉｉ．ｇ）ｉ．の一般的な動態を、ｇ）ｉｉ．で決定した影響を及ぼす量を定量化する複数の項によって拡張する工程と、
用いて考案し、それによって複数のモデルパラメーターを定義する工程。

ｈ．任意で、ｇ）の動態に関して、ｆ）で計算した反応速度に対するモデルパラメーターの値の第１の調整を、各マクロ反応について別々に実行し、調整の質をチェックする工程。

ｉ．任意で、所定の調整の質に達するまで、工程ｇ）およびｈ）を繰り返す工程。

ｊ．モデルパラメーターの値を、ｃ）の測定データに調整する工程。

ｋ．行列Ｌ、ｇ）の動態、およびｊ）のモデルパラメーターの値でモデルを形成し、出力として提示し、および／または、プロセス制御モジュールまたはプロセス開発モジュールに移行する工程。

通常、プロセス制御モジュールは、バイオリアクターを制御するために従来使用される分配制御システムとオンラインでつながっている。

通常、プロセス開発モジュールは、オフラインでのプロセス最適化や実験の計画に使用される。

本発明による生体反応モデル化は、内部代謝プロセスを単純化して表現した代表的なマクロ反応の想定に本質的に基づいている。反応の選択には、生化学的な背景知識とプロセスノウハウの両者が必要である。

この方法の第１の工程で、代謝ネットワークの反応、その化学量論的特性および可逆特性は、使用者によりユーザーインターフェースを介して、または理想的には、オーガニズムに関する背景知識が保存されたデータベースモジュールからオーガニズムおよびその蓄積された代謝経路を選択することによって自動的に入力される。代謝ネットワーク（先行技術では化学量論的ネットワークとも呼ばれる）およびその各反応の特性は、オーガニズムの背景知識を表す。好ましくは、代謝ネットワークは、オーガニズムにとって重要な代謝経路の反応、例えば解糖反応を含む。特に好ましくは、選択には外部反応を含む。本明細書の文脈での外部反応は、細胞の外側の構成要素の少なくとも１つ、典型的には少なくとも１つの入力量および／または少なくとも１つの出力量（生成物、副産物など）を含む。特に好ましくは、代謝ネットワークは、細胞増殖を記述する反応を、例えば内部代謝物から外部生物量への単純化した反応の形式で含む。限定はされないが、実施例における図５および表１は、適用可能な代謝ネットワークを説明する。

エレメンタリーモードは、入力したオーガニズムの代謝経路から計算し、１つ以上の化学量論的ネットワークに結合される。各エレメンタリーモードは、代謝経路から得た反応速度、すなわち、代謝ネットワークの内部および外部反応の線形結合であり、内部代謝物については「定常状態」の状態が満たされ、また、反応の可逆性または不可逆性が考慮される。内部代謝物に関して「定常状態」の状態を考慮に入れた反応の線形結合により、内部代謝物は蓄積し得ない。

本明細書の文脈での内部反応は、細胞内だけで起こる。

内部構成要素の外部化によって、すなわち実際には内部である構成要素を入力量または出力量として分類することにより、外部化された内部構成要素に関連する内部反応を外部反応としてモデル化することが可能であり、したがって、この場合の内部代謝物の「定常状態」の状態を回避することが可能である。

本明細書の文脈でのマクロ反応とは、１つまたは複数の入力量から１つまたは複数の出力量へつながるすべての反応を一緒にまとめたものである。よって、各エレメンタリーモードは、マクロ反応を記述する。Leifheit et al.の方法と比較すると、本明細書の文脈でのマクロ反応は、入力される背景知識をもとに決定される。

エレメンタリーモード（ＥＭ）は、対応するアルゴリズムで構成される、行列Ｅ、好ましくは行列構築のモジュールに結合される。エレメンタリーモード行列を構築するために公知のアルゴリズムを用いてもよい。限定はしないが、ＭＥＴＡＴＯＯＬが一例として挙げられる［Pfeiffer T, Montero F, Schuster S. 1999. METATOOL: for studying metabolic networks. Bioinformatics 15(3):251-257.］

ＭＥＴＡＴＯＯＬでは、入力した内部および外部反応を記述する第１の行列Ｅが作成される。

工程ｂ）では、（外部）化学量論行列Ｎ_ｐを使用して、行列（Ｅ）から、可能なマクロ反応からなる行列Ｋを作成する。

行列ＥのＫへの変換は、Provostにより公知である［Provost A. 2006. Metabolic design of dynamic bioreaction models. Faculte des Sciences Appliquees, Universite catholique de Louvain, Louvain-la-Neuve, Louvain-la-Neuve, p. 81］。

行列Ｋの列ベクトルは、マクロ反応を記述する。行ベクトルは、マクロ反応の構成要素（入力および出力量）を記述する。マクロ反応の化学量論は、行列Ｋへ入力される。

代謝ネットワークの文脈で可能な各反応速度は、これらのマクロ反応の正の線形結合として表現され得る。

ＥＭをプロセスモデルの基礎として使用することは先行技術により公知である（例えばProvost A. 2006. Metabolic design of dynamic bioreaction models. Faculte des Sciences Appliquees, Universite catholique de Louvain, Louvain-la-Neuve, Louvain-la-Neuve, p. 87, p. 118 ff. and Gao, J. et al. (2007). Dynamic metabolic modeling for a MAB bioprocess. Biotechnology progress, 23(1), 168-181)。

更なる工程ｃ）では、オーガニズムとの生体反応に関して利用可能な測定データ（プロセスノウハウ）が入力される。通常は、細胞数、細胞の生存、炭素源（例えばグルコース）、アミノ酸またはＯ_２、生成物および副産物（例えば、乳酸もしくはＣＯ_２）などの基質の濃度、温度および／もしくはｐＨなどのプロセスパラメーターまたは生成物の特徴が決定される。この入力は、例えば、測定データの保存のためのデータベースモジュールから選択し、選択したデータをデータベースモジュールに連結したデータ分析モジュールに送ることにより、使用者によって手動または自動で実行され得る。

これらの測定データから、基質および（副）産物の細胞に特異的な分泌および摂取速度（一緒に特異速度
と呼ぶ）ならびに任意でオーガニズムの増殖および死亡率（μ（ｔ），μ_ｄ（ｔ））が工程ｄ）で計算される。この計算には、補間法を用いて、生細胞数、細胞総数および培地濃度の補間が必要となる。測定された量の経時変化はこれから決定することができる。計算された速度
は、オーガニズムの観察可能な動的挙動に関する経時的な情報を提供する。

上記速度を計算するために、１つまたは複数の異なる補間法を組み合わせて使用し得る。

一例として、Leifheit et al.は、測定された量、例えば、細胞総数、生細胞数、または他の培地濃度の経時変化を、スプライン補間した測定データを利用して測定データから決定することを説明する［Leifheit, J., Heine, T., Kawohl, M., & King, R. (2007). Rechnergestutzte halbautomatische Modellierung biotechnologischer Prozesse [Computer-aided semi-automatic modelling of biotechnology processes]. At-Automatisierungstechnik, 55(5), 211-218］。この方法は、参照により本明細書に組み込まれる。

上記の速度
は、これらの経時変化から計算される：
例えば、オーガニズムの増殖率μ（ｔ）ｍ（ｔ）を、細胞総数Ｘ_ｔ（ｔ）および生細胞数）Ｘ_Ｖ（ｔ）のスプライン補間された値、ならびに、そこから計算可能な細胞総数の経時変化
から、式２使用して計算することができる：
式中、Ｄ（ｔ）は希釈率である。

死亡率μ_ｄ（ｔ）は、μ（ｔ）の経過がわかっている際は、Ｘ_ｖ（ｔ）の経過および生細胞数の経時変化の経過
から、式３を用いて計算できる：
別の構成要素ｉの特異速度ｑ_ｉ（ｔ）は、生細胞数）Ｘ_ｖ（ｔ）および構成要素の濃度）Ｃ_ｉ（ｔ）のスプライン補間された値、ならびに、Ｃ_ｉ（ｔ）のスプライン補間値から決定され得る経時変化の経過
から、式４を用いて計算され得る：
本方法の好ましい実施形態では、工程ｃ）の測定データは、第１の補間の前に次のようにして準備される：細胞によって引き起こされたのではないすべての濃度変化を考慮し、測定データから濃度変化の連続経過を得るために、測定値を移動させる（本明細書では「シフト」と呼ぶ）。濃度測定が移動する量ΔＣ_ｉ（ｔ）は、式５に従って計算することができる：
式中、Ｄ（_τ）は希釈率である。次いで、シフトした濃度経過Ｃ_ｉ，ｓ（ｔ）を式６に従って得る：
したがって、シフトした濃度経過Ｃ_ｉ，ｓ（ｔ）を示す微分方程式は、式４および式６から得られ、式７をもたらす：
測定データのこの準備（シフト）は、特に流加プロセスで、供給が開始または停止されるとき（供給ピーク）に、計算された特異速度の突然の変化を防止する。

図１は、本明細書の文脈における測定データの準備／シフトを示す。

本方法の特別な実施形態では、準備したデータは次いで、Leifheit et al.の方法を用いて細胞総数の勾配
を計算するために使用される。これは微分方程式８に従ってスプライン補間により近似される：
特に好ましくは、溶解因子Ｋ_ｌ（式９）を用いて溶解が計算に含まれる。溶解因子は、例えば、プロセスの経過で一定であると仮定され得る。

シフトした細胞総数Ｘ_ｔ，ｓ（ｔ）の低下は、溶解によって説明することができ、その結果、増殖率μ（ｔ）の負の値を回避することができる。

好ましくは、準備したデータは死亡率μ_ｄ（ｔ）の計算にも使用される。

好ましくは、特異速度
の可能な組み合わせは、例えば図２のように、フラックスマップ図に示される。この提示方法は、計算された特異速度
の良好な概要を提供する。ここで、輪郭線は、どの領域が生理的に重要であるかを示す。

特異速度
および任意でさらなる速度μ（ｔ）およびμ_ｄ（ｔ）が異なる桁数および単位を有する場合、これらは通常、単純化して、同じ単位を有する特異速度ベクトル
に結合される。例えば、グラム[g]単位で測定される高分子の特異速度は、単位
を有する。この高分子の組成を、例えばそのＣ−モル含量に基づいて推定する場合、特異速度が単位
を有するように、その特異速度を[g]から[C−mol]へ変化させて表現し得る。

特異速度
は、方法の次の工程ｅ）のための基礎の１つ、すなわち関連するマクロ反応の選択を形成する。

工程ｅ）において、ＥＭのサブセット（Ｌ）がデータに基づいて選択され、これは数学的品質基準に従ってｄ）の特異速度
および／またはｃ）の測定データを良好に表現することができる。プロセスモデルの複雑性を可能な限り小さくするために、サブセット（Ｌ）内のＥＭの数はできるだけ小さくする必要がある。しかし、サブセットＬは、プロセスノウハウを適切に記述するようにしなくてはいけない。

ＥＭの選択は、ａ）の元のＥＭセット（Ｋ）と比較して解空間のサイズを減少させるが、決定された細胞の生理学的に重要な領域をさらに含む。

図３は、元のＥＭのセット（Ｋ）がサブセット（Ｌ）に縮小した解空間を表現したものである。

工程ｅ）では、計算された特異速度
およびｃ）の測定データは、通常、対応するアルゴリズムで構成された、関連するマクロ反応を選択するためのモジュールに送られる。

工程ｅ）ｉ．において、行列Ｋのデータ非依存性および／またはデータ依存性の事前縮小があらゆる所望の順序で行われる：
データ非依存性の事前縮小は、好ましくは、幾何学的縮小によって行われる。これには、ランダムに選択されたＥＭに関するすべての他のモードに対してコサイン類似度をすべて計算することが伴う。最大類似度を有するＥＭをセットから除去する。この手順は、所定数のＥＭに達するまで繰り返される。所望の数は、通常、方法のために予め定義される。解空間の体積は、制御変数として使用してもよい。驚くべきことに、元の体積と比較して、張られる体積の９０〜９８％、好ましくは９２〜９５％を維持しながら、マクロ反応の数の明確な縮小が可能であることがわかった。

データ依存性の事前縮小は、ＥＭの産出係数（Ｙ^ＥＭ）とｄ）の特異速度
から計算された産出係数（Ｙ^ｍ）を比較することによって行うことができる。ｋ番目のＥＭ産出係数
は、式１０によって、外部代謝物ｉおよびｊの対応する化学量論係数を割ることで決定される。ｋ番目のＥＭについては、行列の成分はＫ_ｊ，ｋおよびＫ_ｊ，ｋである。

化学量論係数Ｋ_ｊ，ｋ＝０である場合、産出係数は決定できない。

産出係数
は、式１１に従い、ｄ）で測定または計算された２つの細胞特異速度間の比率
を意味する。

産出係数Ｙ^ｍから、２つの外部構成要素ｉおよびｊの可能な組み合わせのそれぞれの上限値と下限値を決定することが可能である。例えば、２つの外部代謝物ｉおよびｊの最小産出係数
を下限値として使用し、
の最大値を上限値として使用することが可能であるが、他の限界値も可能である。上記上限値超えまたは下限値未満の産出係数
を有するＥＭは、行列Ｋから除去する。ＥＭの産出係数
が決定できない場合は、行列Ｋに残る。好ましくは、データ依存性の事前縮小に関して１５ページに記載されているように、本発明の方法、「ＮＮＬＳを用いる選択したマクロ反応の反応速度の線形推定」を使用することも可能である。事前縮小の文脈での方法が当該箇所で説明される。データ依存性の事前縮小におけるプロセスデータを方法で使用する利点は、この縮小がプロセスに関連し、したがってより効果的かつ集中的であることである。

工程ｅ）ｉｉ．において、マクロ反応のサブセットはアルゴリズムを使用して選択される：選択は、ｄ）の特異速度
および／またはｃ）の測定データがサブセット（Ｌ）でどれくらい良好に表現され得るかを定量化することが可能な品質基準、および、サブセットを選択するためのアルゴリズムを必要とする。

計算された特異速度
を、サブセットＬで表現するために使用される品質基準は、Soons et al.,によると、式１２に従った、特異速度の平方残差の和）（ＳＳＲ_ｑ）であり得る［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, 7-9 July 2010, 156-161］。

（ＳＳＲ_ｑ）の値は可能な限り小さくなくてはいけない。

（ＳＳＲ_ｑ）の最小化には、以下のような非負最小二乗アルゴリズムを利用して、考慮された各時点t_iに関してベクトル
を予め決定することを必要とする：
式には次の境界条件が伴う：
この方法の利点は、式１２〜１４による計算が、多くのＥＭを有する非常に大きなサブセットに対しても実行され得ることでもある。顕著な欠点は、この計算には計算された特異速度
が必要であることである。比速度は補間された測定値から得られるので、真の値に関しては不確実性が大きい。場合によっては、測定の不正確さが、計算された特異速度
に非常に大きな影響を及ぼす可能性がある。結果として、品質基準ＳＳＲ_ｑも、かなりの不確実性を伴ってしか決定され得ない。表現の質に関する情報に加えて、この方法はまた、式１３および１４による最小化の結果として、サブセットLの反応速度
の推定経過ももたらす。

Leighty, R. et al. が、測定値（濃度測定）を、経時的な体積測定の反応速度の線形推定によって直接近似する別の方法を記述する。線形最小二乗解法を用いて線形最適化問題を解くことによって、反応の経過は、支持点間で線形的に進むと仮定して早く推定することができる［Leighty, R. W., & Antoniewicz, M. R. (2011), Dynamic metabolic flux analysis (DMFA): a framework for determining fluxes at metabolic non-steady state. Metabolic engineering, 13(6), 745-755］。この方法は可逆的なマクロ反応（出典で言及される「フリーフラックス」など）のみに関し、さらに希釈効果（したがって細胞によって引き起こされたのではない濃度変化）を考慮することはできない。マクロ反応の大きさと測定値の寸法が一致しない場合、これらの測定値を反応速度の推定に使用することはできない。これは、例えば外部生物量形成の形式での細胞増殖がマクロ反応の一部であり、細胞の乾燥質量のみが測定値から既知である場合である。したがって、この形式では、この方法は、不可逆的なマクロ反応および流加プロセスの使用には適していない。

参照により本明細書に組み込まれるLeighty et al.らの概念を用いて、本発明に従って準備した（シフトした）データを用いて、この方法を流加プロセスに適用することもできる。さらに、線形最適化問題の境界条件として、マクロ反応の反応速度に対して下限値を加えることによって、この方法は、エレメンタリーモードのような不可逆反応にも使用することができる。マクロ反応および測定値の寸法が一致しない場合、適切な相関によって測定値の寸法をマクロ反応の寸法に調整することができる。Leighty et al.の線形推定と本明細書による改良の組み合わせは、以下、「選択したマクロ反応の反応速度の線形推定」と呼ぶ。

したがって、測定データが元のＥＭセットＫのサブセットＬの選択されたマクロ反応で適切に表現され得るかどうかを検証することが可能である。本方法で決定されたシフトした濃度
と、ここで計算されるシフトされた濃度
との間の式１５による平方残差の最終的な和ＳＳＲ_ｑは、測定データがそのサブセットでどれくらい良好に表現されることができるかを示す。

ＳＳＲ_ｑの値が小さいほど、サブセットＬは良好である。この方法は、特異速度の間違っている可能性のある事前の決定がなくても、サブセットの品質の迅速な検証が可能であるため、Soons et al.の方法よりも流加プロセスのモデル化に特に好ましい。推定された反応速度が支持点間で直線的に進行すると仮定し、測定偏差は反応速度の推定にほとんど影響しない。この方法の欠点は、検証中のサブセットＬのサイズが線形最適化問題の解によって制限されることである。サブセット内の反応の最大数は、利用可能な測定数を支持点の数で割った数に等しい。

表現の質に関する情報に加えて、この方法はまた、サブセットの反応速度
の推定経過ももたらす。

本発明の「選択したマクロ反応の反応速度の線形推定」の好ましい実施形態では、線形最適化問題は、線形最小二乗解法の代わりに、Lawson et al.による非負最小二乗解法（ＮＮＬＳ）を使用することによって解決される。［Lawson, C.L. and R.J. Hanson, Solving Least Squares Problems, Prentice-Hall, 1974, Chapter 23, p. 161.］。これにより、この方法でより大きなサブセットの品質も確認することが可能になる。この場合、マクロ反応の最大数は、利用可能な測定数を支持点の数で割った数よりも顕著に大きい可能性がある。この、「選択したマクロ反応の反応速度の線形推定」と非負最小二乗解法の使用との組み合わせは、以下で、「ＮＮＬＳを用いる選択したマクロ反応の反応速度の線形推定」と呼ぶ。

「ＮＮＬＳを用いる選択したマクロ反応の反応速度の線形推定」の本発明の方法は、さらに、工程ｅ）ｉ）におけるＥＭの事前縮小に関連してさらなるデータ依存性の方法として用いられ得る。この目的のために、マクロ反応の非常に大きなセットKを使用することが可能である。この方法の結果は、第１に、ＳＳＲ_ｃの値であり、第２に、反応速度
の経過である。関連する速度
の値が小さいＥＭは、行列Ｋから削除される。この手順は、所定数のＥＭに達するか、またはＳＳＲ_ｃの値が指定の閾値を超えるまで繰り返す。

サブセットを選択するアルゴリズムは、例えばProvost et al. およびSoons et al.により公知である［Provost A. 2006. Metabolic design of dynamic bioreaction models. Faculte des Sciences Appliquees, Universite catholique de Louvain, Louvain-la-Neuve, Louvain-la-Neuve; 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, 7-9 July 2010, 156-161］。

Soons et alは、２段階最適化法におけるＥＭサブセットの形成を開示する。ランダムに選択された様々なＥＭの場合、ＳＳＲ_ｑの値は上記のようにそれぞれの場合において最小化される。１番小さく最小化された値ＳＳＲ_ｑを持つセットが選択される。しかし、ＥＭの数が大きい場合、可能な組み合わせの数が著しく増加するため、ランダムな選択は効果的ではない。例えば、２００００個のＥＭのセットから１０個の反応を選択する場合、２．８・１０^３６２．８・１０^３６個超えの組合せが存在する。そうすると、最適な組み合わせを見つける確率はかなり小さくなる。品質基準ＳＳＲ_ｑの使用により、この方法は測定の不確実性や測定偏差に対しては弱い。

Provostは、エレメンタリーモードのすべての可能な正の線形結合を
の様々な特定の値に対して決定する代替アルゴリズムを開示し、ここで：ＳＳＲ_ｑ（ｔ＝ｔ_ｉ）である。次に、これらの多数の可能な組み合わせからランダムに組み合わせを選択する。どの場合も、この方法では１つのベクトル
のみを使用し、全体の経時変化は用いない。したがって、完全なプロセスのためのＥＭの選択は不可能である。ランダム選択では、ベクトル
は表現され得るが、残りのプロセスがそれによってどの程度表現され得るかは決定されない。この方法の別の欠点は、全ＥＭの解空間内に存在しないベクトル
を使用できないことである。近似解は決定できない。これは特に不確実性の伴う測定および特異速度の場合には大きな欠点である。品質基準ＳＳＲ_ｑの使用により、この方法も同様に測定の不確実性や測定偏差に対しては弱い。

したがって、この方法のさらなる好ましい実施形態では、関連するマクロ反応を選択するため、すなわちＥＭサブセットＬを選択するために、進化的、特に遺伝的アルゴリズムが使用される。そのようなアルゴリズムは、例えば、Baker et al.により公知である［Syed Murtuza Baker, Kai Schallau, Bjorn H. Junker. 2010. Comparison of different algorithms for simultaneous estimation of multiple parameters in kinetic metabolic models. J. Integrative Bioinformatics:-1-1.］。特に好ましくは、ＥＭの様々な組み合わせについて対応する値ＳＳＲ_Ｃを計算するために用いられる「選択したマクロ反応の反応速度の線形推定」方法での目的関数において、遺伝的アルゴリズムを使用すること可能である。または、ランダムな選択も用いられ得る。工程ｉｉ）の完了後、行列Ｌは必要なマクロ反応を含む（工程ｉｉｉ）。

あってもなくてもよい工程ｉｉｉ）では、ＥＭサブセットＬの有効性が図で検証される。この場合、ＥＭサブセットＬを投影するものとして工程ｄ）のフラックスマップを使用することが可能である。図４は、６つのＥＭのサブセットを投影したフラックスマップを示す。ＥＭサブセットＬが有効である場合、測定データはＥＭサブセット内に位置したままとなる。この表現により、選択の有効性を図で迅速に確認することが可能になる。

さらなる工程ｆ）では、ｄ）の特異速度
および／またはｃ）の測定データを使用して、サブセットＬのマクロ反応の反応速度を計算する。
の計算は、Soons et al.に従い、ｅ）に記載される特異速度
に基づいてなされ得る［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, 7-9 July 2010, 156-161］；好ましくは、
の計算は、本発明の「選択したマクロ反応の反応速度の線形推定」を利用してｃ）の測定データに基づいて行われる。

この方法の工程ｇ）では、マクロ反応の動態が考案される。決定された動態は、プロセス状態のそれぞれの反応速度
に対する動的影響を定量化しなければならない：
決定するモデルパラメーターは動態から見出す。

工程ｇ）ｉ．において、一般的な動態はマクロ反応の化学量論から考案する。マクロ反応ｉの基質については、モノータイプの制限が想定される。最大反応速度に様々な制限
を掛ける：
式中、モノー定数Ｋ_{ｍ，ｋ，ｉ}およびヒル係数ｎ_ｉは等式のパラメーターであり、その第１の値を手動で入力する。通常、モノー定数Ｋ_{ｍ，ｋ，ｉ}はそれぞれの最大測定濃度の１／１０に設定され、ヒル係数ｎ_ｉは、１の値に設定される。反応化学量論からの一般的な動態の決定は、Provostまたはby Gao et al.によって説明される［Provost A. 2006. Metabolic design of dynamic bioreaction models. Faculte des Sciences Appliquees, Universite catholique de Louvain, Louvain-la-Neuve, P/. 126 ; [Gao, J., Gorenflo, V. M., Scharer, J. M., & Budman, H. M. (2007). Dynamic metabolic modeling for a MAB bioprocess. Biotechnology progress, 23(1), 168-181］。この方法は、参照により本明細書に組み込まれる。これらの方法は、反応のそれぞれの基質に対してモノータイプの基質制限を用いる。ProvostまたはGaoは開示していないものの、毒性の生成物による阻害も、この方法を用いた反応化学量論から誘導され得る。

工程ｇ）ｉｉ．において、ｆ）で決定された反応速度
に影響を及ぼす量が決定される。これは、プロセス状態（すなわち、マクロ反応の化学量論からは導き出せないｐＨ、反応器温度、分圧などの生体反応条件を含む）を記述する全ての量を考慮することを伴う。影響を及ぼす量は手動で決定され、例えば、部分的最小二乗法のような統計的方法を用いる。この目的のために、プロセス状態（行列に組み合わされる）とｆ）の反応速度
との間の相関が決定される。

工程ｇ）ｉｉｉ．において、ｇ）ｉｉ．で決定した影響が定量化され、ｉ．での動態が対応する項によって拡張される。ｇ）ｉｉ．で分かった反応速度に対するプロセス状態の量の影響には、項
を与えてもよい。項
は、プロセス状態に応じて、０と１の間の値をもたらす任意の所望の関数である。ｇ）ｉ．で確立した反応の一般的な動態に、この項を掛ける。

例えば、構成要素ｉの濃度と反応ｋとの間に負の相関が見られる場合、これは、ｉの濃度（Ｃ_ｉ）による反応速度ｋの影響を意味する。これは、例えば、Haldaneによる阻害速度論で提供することができる：
式中、Ｋ_{ｌ，ｋ，ｉ}は阻害定数を示し、さらなるモデルパラメーターであり、その第１の値は手動で入力され、通常それぞれの最大測定濃度の１０分の１に設定される。

あってもなくてもよい工程ｈ）において、動態のモデルパラメーター
の値は、ｆ）で決定されたマクロ反応の反応速度
に調整される：
これを、以下でモデルパラメーター値の推定と呼ぶ。この工程では、式２〜４による１つ以上の微分方程式の数値解を省略し得る；通常３〜１０個のパラメーターを持つ独立した群では、モデルパラメーターの値を各マクロ反応ｋに対して個別に調整し得る。調整は、ガウス・ニュートン法といった従来の方法で実行される［Bates DM, Watts DG. 1988. Nonlinear regression analysis and its applications. New York: Wiley. xiv, 365.］。

各マクロ反応に対して別々であるこのモデルパラメーター値の推定は、工程ｉ）および工程ｊ）に関して特に効果的であり、なぜなら、迅速に実行され、また、工程ｊ）で行うモデルパラメーター値のｃ）の測定データへの調整のために改善された開始値を提供するからである。

調整の質は、例えば、式２０による平方残差の和SSR_ｒを用いて計算される：
ＳＳＲ_ｒの値が小さいほど、良好な調整である。あるいは、調整の質は、
およびｒ_ｋの図の比較によって検証される。

あってもなくてもよい工程ｉ）では、ｇ）で選択したマクロ反応の動態の調整の質を確認する。これは、工程ｈ）で計算され、モデルパラメーター値の推定の調整の質を定量化する値ＳＳＲ_ｒに基づく。調整の質が満足のいくものでない場合、所定の調整の質に到達するまで、工程ｇ）およびｈ）を繰り返してもよい。

さらなる工程ｊ）において、ｇ）の動態のパラメーター値は、調整のための通常の方法によりｃ）の測定データに調整することができる。この調整には、工程ｈ）の開始値を使用することが好ましい。モデルパラメーター値の調整は、上述した微分方程式（式２〜４）を組み込んで行われ、例えばガウス・ニュートン法を利用する［Bates DM, Watts DG. 1988. Nonlinear regression analysis and its applications. New York: Wiley. xiv, 365.］またはマルチプル・シューティングアルゴリズムを用いる［Peifer M, Timmer J. 2007. Parameter estimation in ordinary differential equations for biochemical processes using the method of multiple shooting. Systems Biology, IET 1(2):78-88］。

好ましくは、生成物の特徴をモデルに統合することができる。特に好ましくは、副産物または中間体の濃度に依存する生成物の特徴に関して導入することができる。ａ）で入力された代謝ネットワークの外部成分である副産物の濃度は、既にモデルに統合されており、計算することができる。しかし、必要に応じて、他の副産物または中間体を１つまたは複数の別個の代謝ネットワークに結合することもできる。これは、期待される分泌または摂取速度が異なる桁である場合、または指定した代謝プロセスがどの程度の詳細まで検査されるかが異なる場合に有利である。統合モデルの代替として、工程ａ）〜ｊ）では、動態を伴うマクロ反応のセットで、別個の代謝ネットワークの外部構成要素のプロセスの経過を記述する生成物の特徴を計算するための別個のモデルを作成することが可能である。オーガニズムの外側に位置していないが、細胞内でのその蓄積が１つ以上の生成物の特徴に影響を及ぼす副産物または中間体は、ＥＭの計算およびマクロ反応の定式化中に工程ａ）およびｂ）で外部化されてもよく、よって、外部構成要素として分類され得る。細胞内または細胞外の濃度に依存する生成物の特徴の統合は、濃度と生成物の特徴間の定量的または定性的な関係をさらに統合することで行われ得る。

本発明はさらに、本発明による方法を実施するためのコンピュータープログラムまたはソフトウェアを提供する。

本発明による方法を使用して提供されるモデルは、プロセス制御またはプロセス制御の計画およびリアクター内におけるプロセスの調査に使用できる。

限定はしないが、本発明による方法の特別な実施形態を一例として記述する。この方法では、一例として、ハイブリドーマ細胞を用いた発酵のモデルも作成し、その有効性を記述するとおりに試験した。

実施例：ハイブリドーマ細胞培養のモデル化
１ 工程ａ）
代謝ネットワークの形式の背景知識を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.)。当該文献に記載される動物細胞の代謝ネットワークは、３７個の内部および外部代謝物をつなげる３５個の反応を含む（図５を参照；表１を参照）。

反応の可逆性は、公表された当該文献では明示的に記述されていない。その代わり、同文献の「代謝フラックス分析」に関するデータを評価し、不可逆性反応を同定するのに使用した。

内部代謝物の化学量論、すなわち、化学量論係数を含む化学量論行列Ｎおよび反応の可逆性に関する情報とともに、ネットワークのすべてのエレメンタリーモード（ＥＭ）をＭＥＴＡＴＯＯＬ ５．１(Pfeiffer et al. METATOOL: for studying metabolic networks, Bioinformatics 199915 (3), pp. 251-257)を用いて評価した。ＥＭの数は、この場合３０００００超である。

２ 工程ｂ）
計算したＥＭを伴う行列Ｅを工程ａ）で得た。行列Ｎに対しても同様に、行列Ｎ_ｐは外部代謝物の化学量論、すなわち、化学量論係数を含む。化学量論的ネットワークの可能なマクロ反応を式２１とで行列Ｋに結合した：

３ 工程ｃ）
プロセスの測定データを、バッチプロセスの経過におけるハイブリドーマ細胞の発酵の様々な測定量を提供するBaughman et al.から利用した（図６参照）［On the dynamic modeling of mammalian cell metabolism and mAb production. In: Computers & Chemical Engineering (2010) 34 (2), pp. 210-222］。測定データを方法に入力した。

４ 工程ｄ）
ｃ）のスプライン補間された測定値
を用いて、増殖および死亡率ならびに特異的な摂取および分泌速度を計算した（参照図７）。所定の溶解因子Ｋ_ｌ＝０．１を方法に入力して溶解を組み込み、プロセス時間の間、一定にした。測定データのシフトは、この場合、さらに追加物のないバッチプロセスのデータであるため、必要ではなかった。よって、すべての濃度変化が追加物ではなく細胞によって引き起こされるため、データは上向きの経過をたどる。

追加情報を、速度
を計算するのに用いる。Baughman et al.のデータセットでも提供される総生物量（
でのＢＭ）および細胞総数を利用して、
の平均Ｃ−ｍｏｌ含有量を計算できた。Ｃ−ｍｏｌに基づいた増殖率は、式２２から計算できた：
Ｃ−ｍｏｌに基づいた抗体の形成速度も同様に推定できる。このため、抗体のモル組成は、
の正式なモル質量でＣＨ_１．５８Ｏ_０．３１Ｎ_０．２７Ｓ_{０．００４}と推定された。ここで、モル組成は、Villadsen et al.が示すたんぱく質の平均モル組成に対応すると想定される［Bioreaction engineering principles (2011), Chapter 3, Elemental and Redox Balances, p. 73, Springer Verlag, ISBN: 978-1-4419-9687-9］。抗体全体のモル質量は
と推定された。抗体の形成速度を次の式から得た：
速度
の経時変化を用いてマクロ反応を選択することができた。

５ 工程ｅ）
工程ｅ）では、マクロ反応のＥＭ−サブセットを作成し、これを用いてデータセットを可能な限り良好に再現した。これには、工程ｂ）の行列Ｋが必要となった。３０００００超の数のマクロ反応では過剰に大きな数の可能な組み合わせにつながるため、データ依存性の事前縮小を最初に実行した。

この目的のため、工程ｄ）で決定した速度
を２つの外部代謝物の全組み合わせに関する産出係数Ｙ^ｍを計算するのに用いた。産出係数Ｙ_ｉ，ｊの下限値は、９９％の決定産出係数
がその値を上回るように選択された。産出係数の上限値は、９９％の決定産出係数
がその値を下回るように選択された。実施例をもとに、表２では、いくつかの決定限界値と、その限界値内に産出係数
を有するＥＭの割合を示す。全体で、ＥＭの数は、約３０００まで縮小できた。

データ依存性の縮小の後、次いで、データ非依存性の縮小を実行した。この場合、２つのＥＭのコサイン類似性に関して０．９９５の最大値を定義した。第１の反応から始まって、この値を超えるすべてのマクロ反応を行列Ｋから除去した。約５００個のマクロ反応が行列Ｋ（縮小した行列Ｋとも呼ぶ）から残り、それでも約３０００個のＥＭの張る解空間の体積を９５％超えでカバーする。

選択プロセスの前に、Niu et al.による代謝ネットワーク（ａ）の代謝ネットワークの外部代謝物に対応）で指示された構成要素を、ｃ）で得た構成要素の測定濃度で調整した。プロリンは除いて、Baughman et al.が測定した濃度もすべてNiu et al.による代謝ネットワークにおいて考慮される。プロリン濃度の測定結果が利用できるように、外部代謝物としてプロリンを含む別の単純化ネットワークを使用する、または、既存の代謝ネットワークを拡張することが可能である。

計算されたマクロ反応で実際に起こったが、そのデータが利用可能ではない構成要素に関しては、以下では無視した。したがって、行列Ｋの対応する行を行列から削除した。対応する行の削除は、これらの入力または出力が細胞によって使用されていないという意味ではない。代謝ネットワークには存在するが、調整するために利用可能な測定結果がないということである。本実施例では、アルギニン、グルタミン酸、グリシン、ヒスチジン、ロイシン、リシン、メチオニン、アンモニウム、酸素、フェニルアラニン、セリン、トレオニン、トリプトファン、チロシンおよびバリンの入力または出力は無視した。

本方法の次の工程では、背景知識を表す縮小した行列Ｋならびにプロセスノウハウを形成するｄ）の速度
およびｃ）の測定データを用いて、Ｋから、マクロ反応の最小サブセットＬを得る。

本発明の「選択したマクロ反応の反応速度の線形推定」を品質基準として使用した。

速度
と同様、細胞数および抗体の測定値をＣ−ｍｏｌに標準化させた。これは、マクロ反応の寸法が測定値の寸法と一致するように必要である。

サブセットを遺伝的アルゴリズムで選択した。この遺伝的アルゴリズムの目的関数の計算で、「選択したマクロ反応の反応速度の線形推定」で扱う線形最適化問題を解いた。ここで計算した線形最適化問題の最小二乗誤差の最終的な和は、同時に、マクロ反応の特定の選択に対する目的関数の値であった。

ＫからサブセットＬのサイズを選択するのに、Ｌ中の異なる数のマクロ反応で最適化を繰り返し実行した。その数はモデルの複雑性と再現性の正確さの間での妥協である。再現に十分な反応の数を決定するために、サブセットＬの選択を、様々な数のマクロ反応に対して繰り返してもよく、または、反応回数のペナルティ項を遺伝的アルゴリズムの目的関数に直接加えてもよい。この場合、複数の最適化を所定の数のマクロ反応（１０、７、５、４および３）で実行した。遺伝的アルゴリズムで見つかった最小誤差を、マクロ反応の数に対して図９にプロットする。この場合、７個未満のマクロ反応は、プロセス経過を十分良好に表現するには少なすぎることが分かった。選択したマクロ反応を表３に示す。

示すマクロ反応では、ａ）の代謝ネットワークのすべての外部代謝物を示す。しかし、ｃ）での測定データが利用可能であった下線部の外部代謝物のみがモデルの部分を形成する。

６ 工程ｆ）
マクロ反応の選択したセットに関して、経時的反応速度を決定した。この実施例では、本発明の方法、「選択したマクロ反応の反応速度の線形推定」を用いて、図１０で示す測定値を反応速度
の推定により近似させた。本方法の結果は、個別の（体積測定の）反応速度の区分的線形経過である。生細胞数Ｘ_ｖ（ｔ）の補間された経過で分割し、表３に示すマクロ反応の細胞異性反応速度
を得た。こうして得られた反応速度
を図１０に示す。

７ 工程ｇ）
表３で示すすべてのマクロ反応に関して、式２４に従った一般的な動態を想定した：
この場合、動態は、モノー動力学により実現し、すなわち、各基質ｉに関する各反応ｋに対して、式２５に従った制限を導入した：
式中、ｒ_ｋ，_ｍａｘは最大反応速度であり、Ｎ_１は考慮される制限の数であり、Ｃ_ｉは構成要素ｉの濃度であり、Ｋ_{ｍ，ｋ，ｉ}は関連するモノー定数であり、ｎ_ｉは反応次数に関するヒルパラメーターである。これらの値は工程ｈ）およびｊ）で調整される。

さらなる項は、ｆ）の反応速度
の分析から分かる。この実施例では、基質の制限に加えて、式２６による阻害も考慮した。

この制限についても、パラメーターＫ_{ｌ，ｋ，ｉ}およびｎ_ｉの値を調整する必要があった。この反応で用いられる動態の項を表４に示す。

８ 工程ｈ）
各反応速度について、反応速度の経過
は、表４に与えられた動態および動態において考慮した濃度の補間値
を用いて代数的に計算することができた。

これらの動態のパラメーターは、反応iごとに、工程ｆ）で決定した反応速度r_i(t)に別々に調節した。反応iで起こるパラメーターを最適化するための目的関数は、本実施例では次の通りであった：
このようにして調整した、計算された
すべての経過を、対応するｒ_ｋ（ｔ）と一緒に図１１に示す。前者の経過は破線で示され、後者の経過は実線で示されている。経過は定性的に合っていることがわかる。つまり、選択した動態によって、プロセスの動力学も満足に再現できることを意味する。再現が不十分であれば、所望の調整度に到達するまで、迅速な工程ｇ）（他の動態の選択）およびｈ）（パラメーター値の推定）を繰り返すことができるので、この情報はこのモデル化工程において非常に有用である。よって工程ｉ）はここでは必要はなかった。

９ 工程ｊ）
モデルパラメーター値
を、ｃ）の測定データでさらに調節した。このため、全パラメーターを同時に最適化した。さらに、従来は検査されなかったアポトーシスおよび溶解のプロセスも含んだ。これらのプロセスは、生細胞数および細胞総数の進行を記述する微分方程式で必要となる：
アポトーシスを記述するために選択した動態は次の通りである：
溶解速度Ｋ_ｌは、プロセスの間一定と想定された。反応速度のパラメーターに加え、アポトーシスおよび溶解によって導入されるパラメーター、Ｃ_{Ｌａｃ，ｃｒ}（臨界ラクテート濃度）、μ_{ｄ，ｍａｘ}（最大死亡率）、Ｋ_{ｄ，Ｌａｃ}（ラクテート濃度の死亡率に対する影響を記述するモノーパラメーター）およびＫ_ｌ（溶解速度）も本工程で決定した。本実施例では、推定濃度
の経過をデータセットの出発値からＯＤＥシステムの数値解法により決定した。測定濃度
および推定濃度
の間の差を次の目的関数を用いた通常の方法で最小化した：
合計３３個のパラメーター
を用いると、原則、この最適化は、目的関数が多くの局所最適を有するため、実行が難しい。決定論的最適化アルゴリズム、例えばＬｅｖｅｎｂｅｒｇ−Ｍａｒｑｕａｒｄｔアルゴリズムを工程ｈ）の公知のパラメーターの開始値に対して開始する場合、成功率はるかに大きくなる。調整したプロセス経過を図１２に示す。調整したパラメーターを表５に示す。

１０ 工程ｋ）
行列Ｌ、表４の動態および表５の関連パラメーター値を伴うアポトーシスの動態からなるモデルを、出力として作成した。

符号の説明

図１は、測定データのシフトを示す：希釈率（Ｄ（ｔ））の変化がある際に突然変化する測定量（Ｃｉ（ｔ））の実際の経過を示す。シフトした経過（Ｃ_ｉ，ｓ（ｔ））は、細胞によって引き起こされた変化のみによって上昇する。 図２は、２つの特異速度ｑ_１およびｑ_２のフラックスマップを示す。等高線は、測定データにおいて速度の特定の組み合わせが生じる頻度を示す。 図３は、ＥＭの正の線形結合の張る解空間を３次元で示す。完全なセットの解空間は黒で示され、サブセットの解空間は灰色で示される。 図４は、２つの特異速度ｑ_１およびｑ_２のフラックスマップを示す。セットＬのマクロ反応の２次元的な投影がベクトルとして示される。 図５は、Niu et al.の代謝ネットワークの概略図を示す。ここでは、細胞の境界をボックスで示す。ミトコンドリアの細胞内境界を破線で示す。外部構成要素には下付きの「ｘｔ」を付けている。矢印と点線の矢印は反応を示す。 図６は、Baughman et al.のハイブリドーマ細胞を用いた発酵の測定データを示す。生細胞と死細胞の合計から、細胞総数（全細胞）を計算する。略称ＧＬＣ、ＧＬＮ、ＡＳＰ、ＡＳＮ、ＬＡＣ、ＡＬＡおよびＰＲＯは、基質グルコースおよびアミノ酸であるグルタミン、アスパラギン酸、アスパラギン、アラニンおよびプロリンならびに代謝による生成物の乳酸を示す。略語ＭＡＢはモノクローナル抗体の生成物を示し、ＢＭは生物量を示す。 図７は、増殖および死亡率および細胞に特異的な摂取および分泌速度を示す。ｑ_ＭＡＢを除くすべての細胞特異速度は で表される。速度ｑ_ＭＡＢは、
で表される。
図８は、選択された反応セットを伴う「選択されたマクロ反応の反応速度の線形推定」で近似させた濃度を示す。細胞総数（Ｘ_ｔ）および抗体濃度（ＭＡＢ）を、このためにＣ−ｍｏｌに変換した。 図９は、サブセット内のマクロ反応の数（ｎ_Ｒ）に対してプロットした最小誤差を示す。 図１０は、本発明の方法「選択されたマクロ反応の反応速度の線形推定」で決定したマクロ反応の反応速度 を示す。
図１１は、本発明の方法「選択されたマクロ反応の反応速度の線形推定」で決定したマクロ反応の反応速度 （実線）と代数的に計算された反応速度
（破線）を一緒に示す。
図１２は、測定された濃度 （点）とシミュレートしたプロセス経過
（実線）との比較を示す。濃度は、［ｍＭ］で表される。生細胞数および細胞総数（［１０^９ｃｅｌｌｓ／ｌ］でのＸ_ｖ／Ｘ_ｔ）および抗体濃度（［１０^−４ｍＭ］でのｍＡｂ）は例外である。

Claims (13)

1. オーガニズムとの生体反応のプロセス制御および評価のためのコンピューター実行式の方法であって、
   オーガニズムとの生体反応のモデルを作成する工程であって、前記モデルは、行列Ｌと、マクロ反応の動態と、複数の計算された値の複数のモデルパラメーターとを含む、工程を含み、前記工程は、
   ａ．前記オーガニズムの複数の代謝経路と、これらの代謝経路の化学量論的特性および可逆特性とを背景知識として入力し、複数のエレメンタリーモードをこの入力から計算する工程であって、前記エレメンタリーモードは、前記代謝経路のマクロ反応を記述する、工程と、
   ｂ．ａ）のエレメンタリーモードと、すべてのマクロ反応の化学量論的特性および可逆特性とを行列Ｋに記述する工程と、
   ｃ．前記オーガニズムとの生体反応に関する測定データを入力する工程と、
   ｄ．補間法を用いて、ａ）の代謝経路の１つまたは複数の入力量および出力量の分泌および摂取速度を、ｃ）の測定データを用いて前記オーガニズムに関して計算する工程と、
   ｅ．ａ）のエレメンタリーモードのサブセットの形式で、関連するマクロ反応を、
   ｉ．ａ）のエレメンタリーモードのデータ非依存性および／またはデータ依存性の幾何学的縮小と、
   ｉｉ．ｃ）の測定データおよび／またはｄ）の１つまたは複数の速度を使用して、数学的品質基準に従ったアルゴリズムによりｅ）ｉ．の幾何学的縮小から前記サブセットを選択し、前記サブセットを行列Ｌに記述することと、
   ｉｉｉ．任意で前記サブセットを図で示すことと、
   により選択する工程と、
   ｆ．補間法を使用して、ｅ）のサブセットのマクロ反応の反応速度
   ［数１]
   を、ｃ）の測定データおよび／またはｄ）の速度を用いて計算する工程と、
   ｇ．ｅ）ｉｉ．のサブセットのマクロ反応の動態を、以下の中間工程、
   ｉ．一般的な動態をｅ）のマクロ反応の化学量論から計画する工程と、
   ｉｉ．ｅ）のマクロ反応に影響を及ぼす因子をｆ）の反応速度から決定する工程と、
   ｉｉｉ．ｇ）ｉ．の一般的な動態を、ｇ）ｉｉ．で決定した影響を及ぼす因子を定量化する前記モデルパラメーターの値によって拡張する工程と、
   を用いて計画、決定する工程と、
   ｈ．任意で、ｆ）で計算した反応速度に対してｇ）ｉｉｉ．のモデルパラメーターの値の第１の調整を、各マクロ反応について別々に実行する工程と、
   ｉ．任意で、所定の調整の質に達するまで、工程ｇ）およびｈ）を繰り返す工程と、
   ｊ．ｇ）ｉｉｉ．、ｈ）、またはｉ）のモデルパラメーターの値を、ｃ）の測定データに調整する工程と、
   ｋ．前記オーガニズムとの生体反応のモデルを構築すると共に、ｅ）ｉｉ．の行列Ｌ、ｇ）の動態、およびｊ）のモデルパラメーターの値を出力し、および／または、プロセス制御モジュールまたはプロセス開発モジュールに送る工程と、
   ｋ）のモデルを用いて、モデルベースの前記生体反応の予測状態評価および／または前記生体反応の閉ループ制御を行う工程と、
   を含む方法。
2. 工程ｄ）で、前記オーガニズムの増殖率を計算し、かつ、最も好ましくは前記オーガニズムの死亡率も計算する、請求項１に記載のコンピューター実行式の方法。
3. 工程ｇ）で、ｆ）の反応速度の分析に基づいて、動態の個別調整を行う、請求項１または２に記載のコンピューター実行式の方法。
4. 工程ｈ）で、ｇ）の動態のパラメーター値の調節を、複数の調整方法を組み合わせて行う、請求項１〜３のいずれか一項に記載のコンピューター実行式の方法。
5. 工程ｅ）ｉｉ．で、前記マクロ反応のサブセットを選択するのに、選択したマクロ反応の反応速度の線形推定を実行する、請求項１〜４のいずれか一項に記載のコンピューター実行式の方法。
6. 工程ｅ）ｉｉ．で、前記マクロ反応のサブセットを選択するのに、選択したマクロ反応の反応速度の線形推定を、進化的アルゴリズムと組み合わせて実行する、請求項１〜５のいずれか一項に記載のコンピューター実行式の方法。
7. 供給ピークを伴わない一定の消費での仕様を達成するために、供給が開始または終了されるときに、計算された特異速度の突然の変化を滑らかにすることにより、前記測定データが工程ｄ）で補間法を使用する前にシフトする、請求項１〜６のいずれか一項に記載のコンピューター実行式の方法。
8. 工程ｆ）で、選択したマクロ反応の反応速度の線形推定を実行する、請求項１〜７のいずれか一項に記載のコンピューター実行式の方法。
9. 工程ｅ）ｉ．で、データ依存性の事前縮小が実行され、ＮＮＬＳを用いる選択したマクロ反応の反応速度の線形推定の方法をこのために使用する、請求項１〜８のいずれか一項に記載のコンピューター実行式の方法。
10. 工程ｅ）ｉｉｉ．で、前記マクロ反応のサブセットの選択の有効性をフラックスマップで試験する、請求項１〜９のいずれか一項に記載のコンピューター実行式の方法。
11. 工程ｅ）ｉｉ．で、ｅ）ｉ．の事前縮小の選択をｃ）の測定データを用いて実行する、請求項１〜１０のいずれか一項に記載のコンピューター実行式の方法。
12. 請求項１〜１１のいずれか一項に記載のコンピューター実行式の方法を実行するためのコンピュータープログラム。
13. 請求項１〜１１のいずれか一項に記載のプロセス工程を実行するためのソフトウェア。