(57)【要約】【課題】推定すべきパラメータの数を削減することで、より大規模な代謝ネットワークについてもべき乗則モデリングを可能とする代謝ネットワーク推定方法及び代謝ネットワーク推定装置を提供する。【解決手段】解析対象の代謝ネットワークに含まれる代謝物質について得られた時系列濃度データから当該代謝ネットワークを記述するヤコビ行列を推定するステップと、推定されたヤコビ行列と、当該代謝ネットワークの定常状態における各代謝物質の流束データと、各代謝物質の定常状態における濃度データとから、べき乗則モデリングで表される代謝ネットワーク式のパラメータを決定するステップとを含む。【選択図】図２
2008015889.tif 000002
【特許請求の範囲】
【請求項１】
解析対象の代謝ネットワークに含まれる代謝物質について得られた時系列濃度データから当該代謝ネットワークを記述するヤコビ行列を推定するステップと、
推定されたヤコビ行列と、当該代謝ネットワークの定常状態における各代謝物質の流束データと、各代謝物質の定常状態における濃度データとから、べき乗則モデリングで表される代謝ネットワーク式のパラメータを決定するステップと
を含む代謝ネットワーク推定方法。
【請求項２】
上記代謝ネットワーク式は、下記式で表される代謝物質xiの変化量（xi'）を表現した式であり、決定するパラメータが下記式におけるαi、βi、gij及びhijであることを特徴とする請求項１記載の代謝ネットワーク推定方法。
【数１】
2008015889.tif 000029
（ここで、αi及びβiは速度定数であり、gij及びhijは反応次数であり、xiは化合物の濃度を示している。nは代謝ネットワークを構成する代謝物質の数を示す正の整数である。i及びjは代謝ネットワークを構成する代謝物質に割り振られた正の整数である（i,j＝1,…,n））
【請求項３】
代謝物質xiの定常状態における濃度xi,0とし、代謝物質xiへと到る流束の合計をFluxiとすると、決定すべきαi及びβiは以下のように計算され、
【数２】
2008015889.tif 000030
推定されたヤコビ行列をJijとし、代謝物質xjの定常状態における濃度xj,0とすると、決定すべきパラメータgij及びhijが下記式で表現されることを特徴とする請求項２記載の代謝ネットワーク推定方法。
【数３】
2008015889.tif 000031
【請求項４】
対象の代謝ネットワークを構成する各反応について、基質となる物質に関するデータと生成する物質に関するデータとを関連付けて記述したモデル構造データを入力するステップを更に含む請求項１記載の代謝ネットワーク推定方法。
【請求項５】
解析対象の代謝ネットワークに含まれる代謝物質について得られた時系列濃度データから当該代謝ネットワークを記述するヤコビ行列を推定するヤコビ行列演算手段と、
前記ヤコビ行列演算手段により推定されたヤコビ行列と、当該代謝ネットワークの定常状態における各代謝物質の流束データと、各代謝物質の定常状態における濃度データとから、べき乗則表現で表される代謝ネットワーク式のパラメータを決定するパラメータ決定手段と
を含む代謝ネットワーク推定装置。
【請求項６】
上記代謝ネットワーク式は、下記式で表される代謝物質xiの変化量（xi'）を表現した式であり、決定するパラメータが下記式におけるαi、βi、gij及びhijであることを特徴とする請求項５記載の代謝ネットワーク推定装置。
【数４】
2008015889.tif 000032
（ここで、αi及びβiは速度定数であり、gij及びhijは反応次数であり、xiは化合物の濃度を示している。nは代謝ネットワークを構成する代謝物質の数を示す正の整数である。i及びjは代謝ネットワークを構成する代謝物質に割り振られた正の整数である（i,j＝1,…,n））
【請求項７】
代謝物質xiの定常状態における濃度xi,0とし、代謝物質xiへと到る流束の合計をFluxiとすると、決定すべきαi及びβiは以下のように計算され、
【数５】
2008015889.tif 000033
推定されたヤコビ行列をJijとし、代謝物質xjの定常状態における濃度xj,0とすると、決定すべきパラメータgij及びhijが下記式で表現されることを特徴とする請求項６記載の代謝ネットワーク推定装置。
【数６】
2008015889.tif 000034
【請求項８】
対象の代謝ネットワークを構成する各反応について、基質となる物質に関するデータと生成する物質に関するデータとを関連付けて記述したモデル構造データを入力する入力手段を更に含む請求項５記載の代謝ネットワーク推定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、代謝ネットワークを表現する代謝ネットワーク式に含まれるパラメータを決定する代謝ネットワーク推定方法及び代謝ネットワーク推定装置に関する。
【背景技術】
【０００２】
代謝システムを数理モデル化し、解析的手法を用いてその性質や仕組みを理解しようという試みが近年行われている。このアプローチの中で最も基本的なものは、代謝ネットワークを構成する各酵素反応における反応速度を微分方程式で記述する手法である。この手法は正確なシミュレーション、詳細な数理解析が可能だが、モデリングに必要な全てのパラメータを推定しようとすると膨大な生化学実験が必要となる。したがって、この手法は、全ての代謝ネットワークに適用することは困難である。
【０００３】
そこで、膨大な生化学実験を必要とせず、代謝ネットワークの数理モデルを構築する手法として、非特許文献１及び２に開示された、一般質量作用則（GMA）及びS-Systemといったべき乗則モデリングを用いる手法を挙げることができる。特に、S-Systemでは、GMAよりもパラメータ数が少ないべき乗則モデル式を使用する。このため、S-Systemはモデリング・フレームワークとしてより妥当である。
【０００４】
非特許文献３には、代謝系を示すヤコビ行列からS-Systemモデルを構築することができる、定常状態条件を基礎としたアプローチが開示されている。ヤコビ行列は、解が一意に決定されないパラメータの検索範囲を最適化段階で拘束するため、非特許文献３に開示された方法によれば、効果的にパラメータ推定を行うことができる。しかしながら、非特許文献３に開示された方法においても、推定に必要とされるパラメータの数が多いといった問題は未解決のままである。
【０００５】
【非特許文献１】Voit EO, Ferreira AiEN: Computational analysis of biochemical systems: a practical guide for biochemists and molecular biologists. Cambridge; New York: Cambridge University Press; 2000.
【非特許文献２】Torres NV, Voit EO: Pathway analysis and optimization in metabolic engineering. New York: Cambridge University Press; 2002.
【非特許文献３】Diaz-Sierra R, Fairen V: Simplified method for the computation of parameters of power-law rate equations from time-series. Math Biosci 2001, 171(1):1-19.
【発明の開示】
【発明が解決しようとする課題】
【０００６】
そこで、本方法は、上述した実情に鑑み、推定すべきパラメータの数を削減することで、より大規模な代謝ネットワークについてもべき乗則モデリングを可能とする代謝ネットワーク推定方法及び代謝ネットワーク推定装置を提供することを目的とする。
【課題を解決するための手段】
【０００７】
上述した目的を達成した本発明に係る代謝ネットワーク推定方法は、解析対象の代謝ネットワークに含まれる代謝物質について得られた時系列濃度データから当該代謝ネットワークを記述するヤコビ行列を推定するステップと、推定されたヤコビ行列と、当該代謝ネットワークの定常状態における各代謝物質の流束データと、各代謝物質の定常状態における濃度データとから、べき乗則モデリングで表される代謝ネットワーク式のパラメータを決定するステップとを含む。
【０００８】
本発明に係る代謝ネットワーク推定装置は、解析対象の代謝ネットワークに含まれる代謝物質について得られた時系列濃度データから当該代謝ネットワークを記述するヤコビ行列を推定するヤコビ行列演算手段と、前記ヤコビ行列演算手段により推定されたヤコビ行列と、当該代謝ネットワークの定常状態における各代謝物質の流束データと、各代謝物質の定常状態における濃度データとから、べき乗則表現で表される代謝ネットワーク式のパラメータを決定するパラメータ決定手段とを含む。
【発明の効果】
【０００９】
本発明に係る代謝ネットワーク推定方法及び代謝ネットワーク推定装置によれば、入力データとして定常状態における各代謝物質の流束データを使用することによって、べき乗則モデル式に含まれるパラメータのうち、推定すべきパラメータの数を大幅に削減することができる。これにより、代謝ネットワーク推定を行うに際してべき乗則モデル式を高速、且つ高い信頼性で出力することが可能となる。
【発明を実施するための最良の形態】
【００１０】
以下、本発明に係る代謝ネットワーク推定方法及び代謝ネットワーク推定装置を、図面を参照して詳細に説明する。
【００１１】
先ず、本発明を適用した代謝ネットワーク推定装置のハードウェア構成について説明する。代謝ネットワーク推定装置のハードウェア構成は、図１に示すように、中央処理装置（CPU）１０１と、読取り専用メモリー（ROM）１０２と、書き換え可能メモリー（RAM）１０３と、ハードディスクドライブ（HDD）１０４と、ハードディスク（HD）１０５と、ディスプレイ１０６と、インターフェース（I/F）１０７と、キーボード等の入力装置１０８とを備えている。また、各構成部はバス１０９によってそれぞれ接続されている。また、ハードウェア構成としては、図示しないが、例えば、各種インターフェースで接続された外部記憶装置、スキャナ及びプリンタを備えていても良い。また、入力装置１０８としては、キーボードに加えてマウスやタッチパネル式の入力パッド等の装置であっても良い。
【００１２】
また、代謝ネットワーク推定装置は、I/F１０７を介して通信回線網とデータ通信可能であって、通信回線網を経由して各種のデータを入力するものであってもよい。ここで通信回線網は有線及び無線に限定されず、またローカルエリアネットワークでもよいし、インターネットであっても良い。
【００１３】
ここで、CPU１０１は、代謝ネットワーク推定装置の全体の制御を司る。ROM１０２は、後述する処理プログラムを記憶している。RAM１０３は、CPU１０１のワークエリアとして使用される。HDD１０４は、CPU１０１の制御下にHD１０５に対する各種データの書き込み／読み出しを制御する。
【００１４】
次に、代謝ネットワーク推定装置の機能的構成について説明する。代謝ネットワーク推定装の機能的構成は、図２に示すように、解析対象の代謝ネットワークに含まれる代謝物質について得られた時系列濃度データ、当該代謝ネットワークの定常状態における各代謝物質の流束データ及び各代謝物質の定常状態における濃度データを入力する入力処理部２０１と、入力処理部２０１で入力した時系列濃度データから当該代謝ネットワークを記述するヤコビ行列を推定するヤコビ行列演算処理部２０２と、入力処理部２０１で入力した各代謝物質の流束データ並びに各代謝物質の定常状態における濃度データ及びヤコビ行列演算部２０２で取得したヤコビ行列からべき乗則モデル式のパラメータを決定するパラメータ決定処理部２０３と、パラメータ決定処理部２０３で決定したパラメータを含むべき乗則モデル式を出力する出力処理部２０４とから構成されている。なお、入力処理部２０１は、解析対象の代謝ネットワークを構成する各反応について、基質となる物質に関するデータと生成する物質に関するデータとを関連付けて記述したモデル構造データを入力することもできる。また、代謝ネットワーク推定装置の機能的構成としては、上記時系列濃度データ、上記流束データ、上記濃度データ及びモデル構造データを格納したデータベース２０５を有していても良い。
【００１５】
次に、代謝ネットワーク推定装置における入力データについて説明する。代謝ネットワーク推定装置における入力データは、代謝物質について得られた時系列濃度データ、各代謝物質の流束データ、及び各代謝物質の定常状態における濃度データである。また、代謝ネットワーク推定装置において、上述したモデル構造データを入力データに含めて定義することもできる。
【００１６】
代謝物質について得られた時系列濃度データは、例えば、キャピラリー電気泳動-質量分析装置を使用して得ることができる。キャピラリー電気泳動-質量分析装置によれば、解析対象の代謝ネットワークに含まれる多数の代謝物質に関する時系列濃度データを網羅的に取得することができる。代謝ネットワーク推定装置は、キャピラリー電気泳動-質量分析装置から出力された代謝物質の濃度データを直接的に入力してもよい。このとき、キャピラリー電気泳動-質量分析装置は、所定の間隔でサンプリングした試料を順次解析にかけることで、代謝物質の濃度データを時系列的に取得することができる。本手法は定常状態における各代謝物質の濃度データを計測する場合にも用いることができる。また、定常状態流束プロファイルの測定方法は、例えば、Zhao J, Shimizu K: Metabolic flux analysis of Escherichia coli K12 grown on 13C-labeled acetate and glucose using GC-MS and powerful flux calculation method. J Biotechnol 2003, 101(2):101-117及びFischer E, Sauer U: Metabolic flux profiling of Escherichia coli mutants in central carbon metabolism using GC-MS. Eur J Biochem 2003, 270(5):880-891を算用することができる。
次に、代謝ネットワーク推定装置における処理ステップについて説明する。
【００１７】
時系列データからのヤコビ行列の算出
代謝ネットワーク推定装置では、第１のステップとしてヤコビ行列演算処理部２０２において代謝物質の時系列濃度データからヤコビ行列を取得する。例えば、Sorribas et alの方法（Sorribas A, Lozano JB, Fairen V: Deriving chemical and biochemical model networks from experimental measurements. Recent Res Devel in Physical Chem 1998, 2:553-573.）に従ってヤコビ行列を取得することができる。
生化学系においてヤコビ行列は以下のように定義される。
【００１８】
【数１】
2008015889.tif 000003
【００１９】
ここで、Jはヤコビ行列であり、δXは小摂動を示し、且つその要素として濃度Xiを含む。
ヤコビ行列の要素は、線形最小二乗適合法を用いて、摂動された時系列から取得される（Diaz-Sierra R, Lozano JB, Fairen V: Deduction of chemical mechanisms from the linear response around steady state. J Phys Chem A 1999, 103(3):337-343参照）。本代謝ネットワーク推定装置は、定常状態条件下で小摂動に対する線形応答が一時的に生じるといった事実に基づいている（Torres NV, Voit EO: Pathway analysis and optimization in metabolic engineering. New York: Cambridge University Press; 2002参照）。よって本手法は、テイラー級数展開の第一次項によって、小摂動による線形応答が十分に表現可能であるということを数学的基礎として、パラメータ変形を行う。
【００２０】
代謝ネットワークのべき乗則モデル式に含まれるパラメータの決定
代謝ネットワーク推定装置では、第２のステップとしてパラメータ決定処理部２０３において、べき乗則モデリングで表される代謝ネットワーク式のパラメータを決定する。具体的に、べき乗則モデル式は、所謂S-Systemのモデル式であって以下のように表現される。すなわち、S-Systemにおいて、対象となる物質の変化量xi'は生成率と分解率の２つの項から構成される。
【００２１】
【数２】
2008015889.tif 000004
（ここで、αi及びβiは速度定数であり、gij及びhijは反応次数であり、xiは化合物の濃度を示している。nは代謝ネットワークを構成する代謝物質の数を示す正の整数である。i及びjは代謝ネットワークを構成する代謝物質に割り振られた正の整数である（i,j＝1,…,n））
定常状態条件では、以下のように簡潔に表現できる。
【００２２】
【数３】
2008015889.tif 000005
【００２３】
ここで、Fluxは、代謝物質xへ定常状態流束の合計であり、V+及びV-はそれぞれ生成項及び消費項である。代謝物質xiの定常状態における濃度xi,0とし、代謝物質xiへと到る流束の合計をFluxiとすると、定常状態条件において、式(2)中の生成項及び消費項は以下のように表現することができる。
【００２４】
【数４】
2008015889.tif 000006
その結果、αi及びβiは以下のように計算される。
【００２５】
【数５】
2008015889.tif 000007
Xiを以下のように定義すると、
【００２６】
【数６】
2008015889.tif 000008
式(2)のヤコビ行列は以下のように表すことができる。
【００２７】
【数７】
2008015889.tif 000009
【００２８】
定常状態条件において、式(8)はαiに式(5)を代入し、βiに式(6)を代入することで以下のように表すことができる。
【００２９】
【数８】
2008015889.tif 000010
従って以下の式を得ることができる。
【００３０】
【数９】
2008015889.tif 000011
【００３１】
すなわち、ヤコビ行列Jij、定常状態における代謝物質濃度xj,0及び定常状態における流束の合計Fluxが与えられると、式（2）は式(10)のように表すことができる。すなわち、代謝ネットワーク推定装置では、定常状態における代謝物質の流束を使用することによって反応次数gij及びhijを決定することができる。ここで、gij及びhijの一方は0である場合が多いので他方の反応次数は容易に算出することができる。また、gij及びhijのいずれも0でない場合には、式(10)を拘束条件としてパラメータ推定を行うことになるが、パラメータ検索範囲を限定するため、比較的に容易に反応次数を推定することができる。
【実施例】
【００３２】
以下、実施例を用いて本発明をより詳細に説明するが、本発明の技術的範囲は以下の実施例に限定されるものではない。
【００３３】
〔実施例１〕
ケーススタディ１＜生化学的な直線経路＞
本発明を、図３に示す酵母ガラクトース代謝モデル（de Atauri P, Orrell D, Ramsey S, Bolouri H: Is the regulation of galactose 1-phosphate tuned against gene expression noise? Biochem J 2005, 387(Pt 1):77-84.）に適用した。本モデルは、５つの代謝物質及び４つの酵素反応からなる。これら４つの酵素反応は、主としてミカエリスメンテン式によって記述される。反応速度式、システムパラメータ及び初期濃度を以下に示した。
【００３４】
４つの酵素反応に関する反応速度式は
【数１０】
2008015889.tif 000012
である。ここで、本モデルにおける各速度表現は以下のように示される。
【００３５】
【数１１】
2008015889.tif 000013
【００３６】
なお、使用されるパラメータは上記de Atauri Pらの論文を参照して入手することができる。計算された定常状態条件は、X1が0.5 mM (固定)であり、X2が0.146 mMであり、X3が0.00703 mMであり、X4が0.817 mMであり、X5が0.243 mMである。また、経路全体の流束は0.0081 mM・s-1である。
【００３７】
本実施例において、S-Systemモデルの全パラメータを、推定する必要なしに決定することができた。生化学モデルは、以下のS-System式に変換することができる。
【００３８】
【数１２】
2008015889.tif 000014
【００３９】
X1は独立変数であるため、上記の反応速度式からは無視することができる。経路の構成から以下の制約が設けられた。すなわち、h22＝g32であり、h23＝g33であり、h33＝h43＝g53であり、h34＝h44＝g54であり、h54＝g44であり、h55＝g45であり、β2＝α3であり、β3＝β4＝α5であり、α4＝β5である。定常状態濃度は、X1,0＝0.50 mM (固定値)であり、X2,0＝0.146 mMであり、X3,0＝0.007 mMであり、X4,0＝0.817 mMであり、X5,0＝0.243 mMである。定常状態流束は、V1＝V2＝V3＝V4＝0.081 mM/minである。ヤコビ行列は以下のようになる。
【００４０】
【数１３】
2008015889.tif 000015
【００４１】
本発明に係る方法によれば、定常状態の代謝物質濃度、定常状態の流束及びヤコビ行列から線形経路のS-Systemモデルを生成することができる。例えば、g21及びg32は以下のように決定される。
【００４２】
【数１４】
2008015889.tif 000016
ここでh22はg32と等価なので、
【００４３】
【数１５】
2008015889.tif 000017
である。
【００４４】
このモデル化工程において、この種の単純な計算によって全てのパラメータを決定することができる。結果として得られるS-Systemモデルは、
【００４５】
【数１６】
2008015889.tif 000018
となる。
【００４６】
〔実施例２〕
ケーススタディ２＜生化学的な分岐経路＞
図４に示すような生化学的な分岐経路を第２の実施例として検討した。この経路は、生化学的経路における典型的な特徴、すなわち分岐、および生成物阻害によるフィードバック制御を含んでいる。図４に示した分岐経路モデルにおける反応速度式は、
【００４７】
【数１７】
2008015889.tif 000019
である。ここで、速度表現は
【００４８】
【数１８】
2008015889.tif 000020
である。
【００４９】
X3はV3及びV4の阻害剤である。パラメータの大部分は、本発明に係る方法によって取得することができる。しかしながら、１つのパラメータは計算によって決定することができなかったため、パラメータ推定方法を使用した。本モデルは、以下のS-System式に変換することができた。
【００５０】
【数１９】
2008015889.tif 000021
【００５１】
経路の構成から以下の制約が設けられた。すなわち、h11＝g21である。定常状態濃度は、X1,0＝0.067であり、X2,0＝0.049であり、X3,0＝0.081であり、X4,0＝0.041である。また、定常状態流束は、V1＝V2＝0.1であり、V3＝V5＝0.043であり、V4＝V6＝0.057である。ヤコビ行列は以下のようになる。
【００５２】
【数２０】
2008015889.tif 000022
【００５３】
半分以上のパラメータは、定常状態濃度、定常状態流束及びヤコビ行列から単純な計算によって取得することができた。α3、g33、β3及びh33については、このような単純な計算によって決定することができなかった。例えばh11は以下の式によって決定することができたパラメータである。
【００５４】
【数２１】
2008015889.tif 000023
g33は計算によって決定することができなかったパラメータであるため、幾つかの制約によって供される推定を必要とした。α3、h33及びβ3は従属したパラメータとして扱うことができるため、g33を決定することで、これら４つのパラメータを決定することができた。
【００５５】
【数２２】
2008015889.tif 000024
【００５６】
さらに、X3は反応V3の基質であるため、h33は正の値をとると定義することができ、g33の検索範囲を限定することができた。g33の検索範囲は0〜-0.95である。検索範囲が限定された、解が一意に決定できないパラメータが１つであるとき、線形適合アルゴリズムは、オリジナルの反応速度モデルに時系列データを適合する際に使用することができる。本モデル化においては、線形最適化法によってg33を決定した。その結果、以下のS-System式が得られた。
【００５７】
【数２３】
2008015889.tif 000025
上述した計算に使用したパラメータのリストを表１に示した。
【００５８】
【表１】
2008015889.tif 000026
【００５９】
摂動後における一過的動的挙動の比較
本発明に係る方法の汎用性を評価するため、様々な摂動に応答したS-Systemモデルの一過的動的挙動を元々のミカエリス・メンテンモデルによる挙動と比較した。この比較は、平均相対誤差（MRE）を計算することで行った。
【００６０】
【数２４】
2008015889.tif 000027
ここで、nはモデル内の代謝物質の数である。ｍは時系列におけるサンプリング点の数である。x’i(t)は作製されたS-Systemモデルから計算された時系列を示す。xi(t)は元々のミカエリス・メンテンモデルから計算された時系列を示す。本実験では10個のサンプリング点が存在する。サンプリング点の間隔は0.5である。ケーススタディ１におけるX2と、ケーススタディ２におけるX1を0〜200％の範囲で摂動する対象とした。100％の摂動に応答した代謝物質の時系列を図５に示した。両実験において、摂動への応答として、S-Systemモデル及び参照モデルにおける同様な動的挙動を観察することができた。
【００６１】
摂動範囲への応答としてのMREの変化を図６に示した。摂動なしでは初期MREが1%以内であった。良好のケーススタディにおいてMREの緩やかな情報が観察されたものの、200%の摂動であっても、MRE の増加は4%以内であった。
【００６２】
〔考察〕
代謝物質濃度に関する時系列データからのべき乗則モデリングは、目的とする代謝ネットワークのサイズに依存するパラメータ推定が必要となる場合が多い。特に、大規模代謝モデルを構築する場合、問題特異的な単純化を余儀なくされる。本発明では、定常状態流束と、代謝物質の時系列データから導き出される代謝ネットワークのヤコビ行列とを使用することによって、解が一意に決定されないパラメータの数を削減することができた。したがって、本発明は、大規模べき乗則モデリングに適したものと言える。実施例１及び２に示したように、運動方程式によって記述される、実在する２つの生化学パスウェイモデル（直線経路と分岐経路）を使用して本発明の効果を検証することができた。
【００６３】
実施例１に示したように、１２個のパラメータを有する線形代謝パスウェイのモデル化するS-Systemにおいて、本発明に係る方法では全てのパラメータを実際に決定することができた。これに対して、Diaz-Sierra and Fairenの論文に記載された方法では、エラー訂正パラメータを含む少なくとも４つのパラメータについては、一意に解を決定することができなかった。実施例２に示したように、１６個のパラメータを有する分岐経路の場合、本発明に係る方法によれば１つを除く他の全てのパラメータを決定することができた。これに対して、Diaz-Sierra and Fairenの方法ではエラー訂正パラメータを含む少なくとも１２つのパラメータについては、一意に解を決定することができなかった。これら実施例によれば、本発明のように、代謝物質濃度に関する時系列データからのS-Systemモデルの構築方法は、解を一意に決定できないパラメータの数を、従前の方法よりもより効果的に削減できることがわかった。さらに、摂動に対する応答実験により、200％の摂動範囲であってもMREが3％程度であることから、本発明に係る方法によって構築されたモデルが、参照モデルと同様な動的挙動を示すことが示された。本発明に係る方法によって生成されたS-Systemモデルは、広範な定常状態点について正確なシミュレーションを行うことができる。この制限は、代謝経路のモデリング及び分析を妨げるものではない。すなわち、本発明は、多くのべき乗則代謝モデリング（例えば、Voit EO, Ferreira AiEN: Computational analysis of biochemical systems: a practical guide for biochemists and molecular biologists. Cambridge; New York: Cambridge University Press; 2000及びTorres NV, Voit EO: Pathway analysis and optimization in metabolic engineering. New York: Cambridge University Press; 2002に開示されたモデリング方法を参照）に適用可能である。
【００６４】
実験ノイズに対する信頼性は、モデリング法の実際の適用場面においては重要な要求事項である。ヤコビ行列は、時系列におけるノイズに対してより感受性が高いので、ヤコビ行列に含まれる数値誤差に応答したロバスト是正が、その使用に関するモデルの構築において重要になる。Diaz-Sierra and Fairenの方法において、誤差訂正パラメータは、実験ノイズの効果を削減するためS-Systemに組み込まれていない。しかしながら、このアプローチでは、解が一意に決定されないパラメータの数が増加してしまう。時系列データにおける実験誤差を削減する別法としては、非線形ニューラルネットワークモデルを使用した適切勾配推定が効果的である（Voit EO, Almeida J: Decoupling dynamical systems for pathway identification from metabolic profiles. Bioinformatics 2004, 20(11):1670-1681）。また、この適切勾配推定は、エラー制御された時系列データを提供することができる。エラー制御された時系列データによれば、高精度のヤコビ行列を取得することができる。一般的な摂動による特異的効果に関する適切推定を取得する方法については、Schmidt H, Cho KH, Jacobsen EW: Identification of small scale biochemical networks based on general type system perturbations. Febs J 2005, 272(9):2141-2151において議論されている。そして、この方法によれば、摂動条件で得られた時系列代謝物質濃度といったデータをより精緻に分析できる。
【００６５】
実験ノイズに対する本発明に係る方法の信頼性を精査するため、数値エラーをヤコビ行列及び定常状態流束に手動で入れ込んだ場合における計算誤差を測定した。シミュレーションされたS-Systemの時間曲線と元々のミカエリス・メンテンモデルとのMREを表２に示した。
【００６６】
実施例２で用いた代謝ネットワークについて、S-SystemモデルのMREと元々のモデルのMREとを測定した。各モデルから得られる時系列データは、S-Systemモデルと元々のミカエリス・メンテンモデルとの間における動的応答の違いを検査するため、状態変数の摂動を含んでいる。ヤコビ行列及び定常状態流束は、100％の数値誤差を含むように再計算された。ヤコビ行列に対する数値誤差はヤコビ行列の全ての要素に等価的に挿入した。摂動対象をX1とした。時系列データは、X1を50％増加させるように摂動した場合のS-Systemモデルから取得した。計算のため10点をサンプル点とした。時系列データにおけるサンプル点の間隔は0.5秒とした。MREは、S-Systemモデルにおける時系列データ及び元々のミカエリス・メンテンモデルにおける時系列データとから計算した。
【００６７】
実施例２に示した生化学的な分岐経路において、ヤコビ行列は100%まで増加し、３つの定常状態流束J1-2、J3-5及びJ4-6は、それぞれ、V1及びV2を介した流束、V3及びV5を介した流束及びV4及びV6を介した流束である
【００６８】
【表２】
2008015889.tif 000028
【００６９】
ヤコビ行列は数値誤差に対して非常に感受性であることは明かである。数値誤差がヤコビ行列の反応次数に直接影響するためである。さらに、定常状態流束プロファイルデータには実験誤差が含まれている場合がある。しかしながら、本発明に係る方法によって作製されたモデルは、ヤコビ行列及び定常状態流束における誤差に対して比較的堅牢である。これは、本発明に係る方法が実際上、利用しやすいことを示唆している。
【図面の簡単な説明】
【００７０】
【図１】本発明を適用した代謝ネットワーク推定装置のハードウェア構成の一例を示すブロック図である。
【図２】本発明を適用した代謝ネットワーク推定装置の機能的構成の一例を示すブロック図である。
【図３】実施例１で使用した生化学系直線経路を示す模式図である。
【図４】実施例２で使用した生化学系分岐経路を示す模式図である。
【図５】(a)及び(b)は実施例１の代謝ネットワークにおいて摂動対象をX2としたときの代謝物質濃度の変化を、それぞれ参照モデル及び本発明のS-Systemモデルについて解析した結果を示す特性図であり、(c)及び(d)は実施例２の代謝ネットワークにおいて摂動対象をX1としたときの代謝物質濃度の変化を、それぞれ参照モデル及び本発明のS-Systemモデルについて解析した結果を示す特性図である。
【図６】(a)は実施例１の代謝ネットワークにおいて摂動対象をX2とし、摂動の範囲を0〜200％としたときの平均濃度誤差(MRE)の変化を示す特性図であり、(a)は実施例２の代謝ネットワークにおいて摂動対象をX1とし、摂動の範囲を0〜200％としたときの平均濃度誤差(MRE)の変化を示す特性図である。
【符号の説明】
【００７１】
101…CPU、102…ROM、103…RAM、104…HDD、105…HD、106…ディスプレイ、107…I/F、108…入力装置、109…バス、201…入力処理部、202…ヤコビ行列演算部、203…パラメータ決定処理部、204…出力処理部、205…データベース
【図１】
2008015889.tif 000035
【図２】
2008015889.tif 000036
【図３】
2008015889.tif 000037
【図４】
2008015889.tif 000038
【図５】
2008015889.tif 000039
【図６】
2008015889.tif 000040