JP5773871B2 - 生物ネットワークのコンピューター実装されるモデル - Google Patents

Description

本発明は、生物ネットワークの反応速度モデルを生成する、コンピューター実装される方法であって、
（ａ）ネットワークトポロジーを選択することであって、該トポロジーのノードは生物学的実体を表し、該トポロジーのエッジは該実体間の相互作用を表す、ネットワークトポロジーを選択すること、
（ｂ）前記相互作用に反応速度則及び反応速度定数を割り当てること、並びに
（ｃ）前記生物学的実体に開始濃度を割り当てること、
を含み、
（ｉ）前記反応速度定数の一部分、及びそれとは独立して前記開始濃度の一部分が実験データであり、
（ｉｉ）前記反応速度定数の残りの部分、及びそれとは独立して前記開始濃度の残りの部分がランダムに選択される、方法に関する。

本明細書において、特許出願及び製造者のマニュアルを含む多数の文書が引用される。これらの文書の開示は、本発明の特許性に関連するとは見なされないものの、参照によりその全体が本明細書に援用される。より詳細には、参照される全ての文書は、各個々の文書が参照によりその全体が本文書に援用されると具体的に個々に指示された場合と同程度、参照により援用される。

ｉｎ ｓｉｌｉｃｏ生物学は、複雑度が高まり続ける生物システムのシミュレーションを可能にしている。これは、疾患に関連した過程への詳細な洞察をもたらし、薬剤開発過程を加速させると共に、治療に対する反応を予測するという可能性を保持している。しかしながら、生物学的過程の詳細な動的モデルを作成する際、実験データが限られていることが主なボトルネックとなる。生物学的経路、ネットワーク、及び過程に関与する分子は多くの場合かなりの程度まで知られているが、これは多くの事例において、これらの分子間の相互作用を定量的に説明する反応速度定数に当てはまらない。未知のパラメーターを推定することは、十分であるとはみなされない。なぜなら、推定は任意であると共にモデルにバイアスをもたらす場合があるためである。

本発明の根底にある技術的課題は、特に実験データがモデルをパラメーター化するのに十分でない事例において、生物システムの時間依存的挙動を予測する手段及び方法を提供することである。

したがって、本発明は、生物ネットワークの反応速度モデルを生成する、コンピューター実装される方法であって、
（ａ）ネットワークトポロジーを選択することであって、該トポロジーのノードは生物学的実体を表し、該トポロジーのエッジは該実体間の相互作用を表す、ネットワークトポロジーを選択すること、
（ｂ）前記相互作用に反応速度則及び反応速度定数を割り当てること、並びに
（ｃ）前記生物学的実体に開始濃度を割り当てること、
を含み、
（ｉ）前記反応速度定数の一部分、及びそれとは独立して前記開始濃度の一部分が実験データであり、
（ｉｉ）前記反応速度定数の残りの部分、及びそれとは独立して前記開始濃度の残りの部分がランダムに選択される、方法に関する。

本明細書において用いられる際の用語「モデル」は、生物システムのｉｎ ｓｉｌｉｃｏ表現を指す。「反応速度モデル」は、生物システムの時間依存的挙動を記述することが可能なモデルである。時間依存的挙動を予測するのに必要な要素は、生物システムの構成要素の転換を含む、生物システムの構成要素間の相互作用を支配する反応速度則及び関連付けられた反応速度定数を含む。これらの構成要素は、本明細書において、「生物学的実体」とも呼ばれる。用語「生物学的実体」は、生物システムにおいて発生し得る任意の分子を含む。好ましい生物学的実体は、下記でさらに詳述する生体分子である。構成要素である生物学的実体は、生物システムのｉｎ ｓｉｌｉｃｏ表現として(an)モデルを表現する。本発明によるモデルは、生物学的実体の開始濃度をさらに含む。反応速度則、反応速度定数、及び開始濃度は合わせて、上記生物ネットワークの時間依存的挙動の予測を可能にする。用語「割り当てる」は、シミュレーションの開始時に、或る特定の特性又は数値を定めるか又は設定することを指す。反応速度則及び反応速度定数は、シミュレーション中に変化しないことが好ましいが、シミュレーション中に想定されるような生物学的実体の濃度は、それぞれの開始濃度と異なる場合があることは自明である。本発明が関係する生物システムは生物ネットワークである。好ましい生物ネットワークは、全ての細胞内相互作用ネットワークを含む。それらの例は、シグナル伝達ネットワーク、転写制御ネットワーク、代謝ネットワーク、感覚及び恒常性ネットワーク、分解ネットワーク、調節ネットワーク、並びにそれらの組み合わせである。

好ましい生物ネットワークは、たとえば受容体リガンド作用、密着結合のような透過性接触、宿主病原体相互作用、及び細胞又は生物間の任意の他の相互作用によって媒介された全ての細胞間相互作用ネットワークも含む。それらの例は、細胞の成長及び分化ネットワーク、血管新生ネットワーク、創傷治癒ネットワーク、炎症及び免疫反応ネットワーク、並びに、結果として機能性組織、器官、生物、及び生物の集合体が生じる、細胞間又は生物間の相互作用の複合ネットワークである。

また、ネットワークは、「グラフ」と呼ばれ、「グラフ」として表されてもよい。ネットワーク又はグラフの例が、本明細書に添付された図１に示されている。より具体的には、当該技術分野において既知であるように、ネットワーク又はグラフはノード及びエッジを含む。ノード及びエッジは合わせてネットワークのトポロジーを形成する。上記ネットワークのノードは、上述した生物学的実体のｉｎ ｓｉｌｉｃｏの対応物であり、上記ネットワークのエッジは上述した実体間の相互作用のｉｎ ｓｉｌｉｃｏの対応物である。本明細書において用いられる際の用語「相互作用」は、任意の種類の相互作用、特に上記相互作用に関与する生物学的実体の量又は濃度に影響を与える場合がある相互作用を指す。より具体的には、用語「相互作用」は、おそらく１つ又は複数のさらなる生物学的実体の影響下における、１つ又は複数の所与の生物学的実体の、１つ又は複数の異なる生物学的実体への転換を含む。他の好ましい相互作用は、たとえば１つ又は複数の他の生物学的実体の作用、存在、又は欠如の結果としての、１つ又は複数の生物学的実体の量又は濃度の減少又は増加を含む。さらに別の好ましい相互作用は、２つ以上の生物学的実体からの複合体の形成である。

換言すれば、本発明による相互作用は、反応を伴うか又は引き起こす。本発明による反応は質量作用反応速度を用いてモデル化することができるが、一般に、任意の他の適切な反応速度則に従うことができる。当該技術分野において既知であるように、質量作用反応速度は、所与の反応に関与する生物学的実体の濃度及び反応速度定数に依拠する。詳細に関しては下記を参照されたい。

従来技術によれば、生物システムの反応速度のモデル化には、関与する全ての生物学的実体又は反応物の全ての反応速度則、反応速度定数、及び開始濃度を知る必要がある。本発明によれば、反応速度モデルであって、反応速度定数の一部のみ及び開始濃度の一部のみが実験データであるか又は実験データから導出され、反応速度定数の残りの部分及びそれとは独立して開始濃度の残りの部分がランダムに、すなわち実験データに依存することなく選択される、反応速度モデルが提供される。この手法は、モンテカルロ手法とも呼ばれる。

反応速度定数を求めるのに適した実験データは、相互作用に関与する生物学的実体の経時変化を含む。しかしながら、上述したように、そのような情報は生物ネットワークの多くの相互作用の場合に利用可能でないか生成するのが困難である。開始濃度に関する実験データは、シミュレートされる生物ネットワークの自然に発生する対応物内、すなわち、たとえば細胞内で測定を実行することによって取得することができる。

既知の反応速度定数及び既知の開始濃度の使用は互いに独立している。したがって、本発明は、全ての開始濃度が実験データであり、上記反応速度定数の一部がランダムに選択される実施形態、並びに全ての反応速度定数が実験データであり、上記反応速度定数の一部がランダムに選択される実施形態を含む。

驚くべきことに、本発明者らは、反応速度定数の一部分又はさらには全てが知られていない事例において、生物ネットワークが反応速度定数の特定の選択に関してロバストであることに気がついた。これは、特に、生物ネットワークによって想定される定常状態及び平衡にあてはまる。特に、反応速度定数を規定する実験データが存在しないときであっても、生物ネットワークの時間依存的挙動によって、偶然による予測をはるかに超える、或る程度の再現可能な予測が生成されることがわかった。これに関して、下記でさらに詳述される、ｉｎ ｓｉｌｉｃｏモデルの統計学的有意性を判断する方法に言及する。同じ観測及び考察を、実験的に知られた開始濃度の部分的欠如又は完全な欠如にも準用する。

本発明は、さらに、生物ネットワークにおいて、生物学的実体の濃度を時間の関数として予測する方法であって、
主要な実施の形態に記載の方法によって生物ネットワークのモデルを生成すること、及び
（ｄ）連立微分方程式を解くことであって、該微分方程式は前記生物学的実体の前記濃度の時間依存性を規定する、連立微分方程式を解くこと、
を含み、それによって前記濃度を時間の関数として取得する、方法に関する。

生物学的実体又は反応物の量又は濃度の時間依存性を、微分方程式によって、より具体的には常微分方程式（ＯＤＥ）によって表すことは、当該技術分野において既知である。

より具体的には、１つの生物学的実体の量又は濃度を制御する相互作用は一般に、以下のようにモデル化することができる。２つの正入力Ａ_１及びＡ_２、並びに１つの阻害性入力Ｉ_１を有する生物学的実体を考える。いずれの正入力も量又は濃度を増加させるのに十分である（Ａ_１∨Ａ_２）、一方、１つの阻害性入力の活動は標的の量又は濃度を減少させるのに十分である（（Ａ_１∨Ａ_２）∧¬Ｉ_１）。いくつかの事例では、ネットワークトポロジー及び実験データに基づいて他の形で相互作用を規定することが好ましい場合がある。この表記において相互作用はブール形式で定式化される。

ＯＤＥモデルの場合、質量作用反応速度が用いられる。分解、輸送、シグナル伝達、及び複合体の会合／解離、並びに翻訳過程及び翻訳後過程のような相互作用を、以下の形式の質量作用反応速度を用いてモデル化した。

上述した反応のグループごとに、１つの偏在性パラメーター（ubiquitous parameter）ｋ_{ｒｅａｃｔｉｏｎｇｒｏｕｐ}を用いた。

ブールモデルの定式化を速度則に移すために、［２８］から導出される手法を用いることができる。特に、以下の基本モジュール、すなわち、実体Ｇに対する正入力Ａの場合、

及び負入力Ｉの場合、

を用いることができる。定数ｃ_Ｘ及びｋ_Ｘは、各実体の調節的役割の個々の特徴を表す。ここで、ｋ_Ｘは阻害剤がないときの活性強度に対応するのに対し、ｃ_Ｘは活性に有意な変化を生成するのに必要な入力の量又は濃度を定める。基本モジュールを組み合わせて、調節的相互作用を含む複合相互作用を定式化することができる。これは、乗算（ブール演算子ＡＮＤに対応する）又は加算（ブール演算子ＯＲに対応する）を用いて行われる。

このため、モデルが遺伝子発現に関連すると想定すると、２つの活性化剤及び１つの阻害剤を用いた上述の例は、ｖ_{ｔｒａｎｓｌａｔｉｏｎ}＝（ｋ_Ａ１・［Ａ１］

となる。

ＯＤＥモデルは、事象を用いて外部入力をオン及びオフにする。濃度を直接変更する代わりに、外部入力を記述するために活性化及び阻害性Ｈｉｌｌ反応速度を用いることができる。これらの反応速度は、いくつかの活性化剤又は阻害剤に依拠するのではなく、シミュレーション時間に依拠する。外部入力の濃度の変化は、以下の微分方程式によって与えられる。

ここで、ｋ_ｄｅｇ・［ｘ］は分解項である。Ｓ_１∈０，１であるため、２つの他の項のうちの一方のみが０に等しくない。このため、Θ及びＳ_１を変更することによって、外部入力が時点Θにおいて活性化される（Ｓ_１＝１の場合）か、又は時点Θにおいて不活性化される（Ｓ_１＝０の場合）。事象を用いてＳ_１及びΘを変更することによって、外部入力の活動を制御することができる。

本発明の根底にある数学的概念及び方法論は、クーン（Kuhn）他の刊行物（２００９）においても詳細に説明されている。クーン他の刊行物（２００９）、この文脈では特に、クーン他（２００９）の「方法」と題するセクションが参照によりすべて援用される。

本発明の方法の好ましい実施の形態では、前記反応速度定数の前記残りの部分は確率分布から選択され、それとは独立して、前記開始濃度の前記残りの部分は確率分布から選択される。それぞれの確率分布は同じであるか又は異なっている場合がある。

用語「確率分布」又は「確率密度関数」は、当該技術分野において既知である。この用語は、特定の事象、この事例では反応速度定数の特定の値又は開始濃度の特定の値を、その発生確率に関連付ける。反応速度定数ごとに、それぞれに利用可能な知識の度合いを反映して選択された確率分布から反応速度定数がランダムにサンプリングされる。それとは独立して、開始濃度に関して選択された確率分布から開始濃度がランダムにサンプリングされる。２つ以上の生物学的実体の開始濃度が選択される事例では、開始濃度を選択するのに同じ確率分布又は異なる確率分布を用いることができ、ここでもまた利用可能な部分的（又は完全な）知識を反映する。さらに説明すると、事前知識が一切ないとき、確率分布は好ましくは全ての未知のパラメーターについて同じである。用語「パラメーター」は反応速度定数及び開始濃度を含む。限られた知識若しくは大まかな知識、又は利用可能な類似の相互作用からの知識しか存在しない事例では、スケーリング係数、及び／又はそのようなタイプの情報を反映するための分布関数の幅の変更によって、その種の知識を考慮に入れることができる。さらに、下記でさらに詳述するように、反応速度則はランダムに、好ましくはここでもまた利用可能な知識に依拠する確率を用いて選択することができる。一般的に言えば、本発明のこの前進法はパラメーター推定の既知の過程とは異なる。パラメーター推定では、モデルパラメーターは、観測に適合する最適パラメーターセットを求める目的で、数学的方法によって推定される。提案される手法では、パラメーターは反復的にランダムに選択され、生成された観測値の有意性が統計学的方法を用いて判定される。

好ましい実施の形態では、前記分布は、対数正規分布である。対数正規分布において、反応速度定数の対数は正規分布であり、すなわちガウス分布に従う。確率分布の妥当性は、分野における用途及び事前知識に依拠する。さらなる適切な確率分布には、一様分布、指数分布、ポアソン分布、二項分布、コーシー分布、ベータ分布、及びガウス分布が含まれる。

本発明による、生物学的実体の濃度を予測する方法のさらなる好ましい実施の形態では、前記反応速度定数の前記残りの部分をランダムに選択し、前記開始濃度の前記残りの部分をランダムに選択し、その後前記連立微分方程式を解くことが反復的に実行される。

この実施の形態は、反応速度定数の異なるセットに対する生物ネットワークの応答の評価を可能にする。反応速度定数の異なるセットは、それらの少なくとも一部が次々とランダムに選択される。驚くべきことに、生物ネットワークの反応速度挙動は、限られた数のパラメーターに依拠し、かつ、ほとんどの反応速度定数の異なるランダムな選択は、生物ネットワークの時間依存的挙動に或る特定の影響を及ぼすものの、上記時間依存的挙動、特に定常状態又は平衡状態を根本的に変えることがないことがわかり、本明細書に添付した実施例において立証されている。

実験的に求められた反応速度定数が完全に欠如しているとき、かつ／又は実験的に求められた開始濃度が完全に欠如しているときであっても、十分なモンテカルロモデル化によって再現可能な予測を達成することができ、これらの予測を既存の知識によって検証することができることは、特に驚くべきである。換言すれば、実験データである反応速度定数の数がゼロであり、生物ネットワークにおける相互作用を規定する全ての反応速度定数がランダムに選択されるという結果になることが意図的に予想される。特定の理論に拘束されることなく、生物ネットワークの特定のトポロジー、すなわち、ノード及びノード間の接続が（ノード間の接続に関連付けられる反応速度定数を無視して）、潜在的な結果の空間を、反復的なモンテカルロモデル化によって系の時間依存的結果の空間の効果的な調査が可能になる範囲で制限すると見なされる。

好ましい実施の形態では、反応速度定数のランダム選択及び常微分方程式系の解決は、複数回行われ、理想的には、利用可能な計算ハードウェアによって制限される実行可能な回数行われる。本発明者らの経験では、現在の計算上の制限に基づいて、１０回〜１０００回の実行が最適であることがわかっているが、１０，０００回以上実行すると、精度が増大するという形でさらなる増分値が生成され続けることもわかっている。

好ましい実施の形態では、全ての利用可能な実験的に導出された反応速度定数及び開始濃度がシミュレーションにおいて用いられ、本発明者らの経験では概して、既知の反応速度定数が一般的にシミュレーションの精度を改善する。しかしながら、別の好ましい実施の形態では、或る特定の既知の反応速度定数は、確率分布から選択された反応速度定数と選択的に置き換えられる。これは、実験的に導出された定数の精度に関して懸念がある事例において、又は実験的に導出された定数が系が定常状態に達するのを妨げる場合に行うことができる。

さらに好ましい実施の形態において、特定の生物システムのシミュレーションのための生物学的実体の反応速度定数及び開始濃度が、同様の系を用いた以前のシミュレーションから導出される。同様の系とは、分析中の系に近い系、たとえば同じであるが特異摂動を伴う生物システムを指す。用語「摂動」は、下記でさらに定義される。定常状態に達するために、以前のシミュレーションが複数のパラメーターセットを用いて実行される。次に、生物学的知識に従って、これらの定常状態のサブセットが選択される。生物学的知識とは、モデル及び既知のパラメーター値によって再現することができる既知のモデル予測を指す。次に、パラメーターセットのこれらのサブセットがシミュレーションに用いられる。

さらに好ましい実施の形態では、反応速度定数はシミュレーションの前に適切な方法によって推定される。

好ましい実施の形態では、上記反応速度定数の少なくとも１０％、少なくとも２０％、少なくとも３０％、少なくとも４０％、又は少なくとも５０％、及びそれとは独立して、上記開始濃度の少なくとも１０％、少なくとも２０％、少なくとも３０％、少なくとも４０％、又は少なくとも５０％が実験データである。明らかに、実験データである、少なくとも６０％、少なくとも７０％、少なくとも８０％、又は少なくとも９０％の反応速度定数を用いることも予想される。

本発明の方法のさらに好ましい実施の形態では、前記反応速度則のさらに一部が知られており、該反応速度則の残りの部分はランダムに選択される。

この実施の形態は、未知の反応速度定数及び／又は未知の開始濃度に加えて、モデルの生物学的実体間の相互作用を支配する反応速度則が部分的に未知である状況にまで及ぶ。反応速度則のランダムな選択は、離散確率分布から実行することができる。確率分布は、反応速度則が離散変数であることの結果として離散している。生物学的実体Ａの生物学的実体Ｂへの転換の反応速度則が未知である事例においてさらに説明すると、反応速度則は、一次反応速度則に関して５０％の確率をもたらし、二次反応速度則に関して５０％の確率をもたらす確率分布から選択することができる。当然ながら、上述したように、大まかであるか又は類似の相互作用から導出される知識から利点を得ることができる。換言すれば、類似の事例において相互作用の９０％が一次反応速度則によって支配され、１０％が二次反応速度則によって支配されることが知られている(it is know)場合、分布は一次反応速度則に関して９０％の確率をもたらし、二次反応速度則に関して１０％の確率をもたらすことができる。

より好ましい実施の形態では、上記反応速度則の少なくとも１０％、少なくとも２０％、少なくとも３０％、少なくとも４０％、又は少なくとも５０％が実験データから導出される。明らかに、実験データから導出される、少なくとも６０％、少なくとも７０％、少なくとも８０％、又は少なくとも９０％の反応速度則を用いることも予想される。

本発明による、生物学的実体の濃度を予測する方法のさらに好ましい実施の形態では、 （ｅ）前記方法は１つ又は複数のさらなる生物ネットワーク上で付随して実行され、
（ｆ）生物学的実体の濃度は、選択された時点において前記生物ネットワーク間で交換される。１つの生物学的実体の量若しくは濃度、又はより多くの若しくは全ての生物学的実体の量若しくは濃度は交換することができる。その量又は濃度が交換される好ましい生物学的実体には、成長因子、サイトカイン、及びホルモンのような細胞間シグナル伝達分子が含まれる。「量の交換」及び「濃度の交換」は、上記さらなる生物ネットワークのうちの１つ、複数、又は全てに対し、上記量又は濃度を利用可能にすることを指す。上記量がさらなる生物ネットワークに利用可能にされると、該さらなるネットワークにおける相互作用を支配する反応速度則に依拠して、それらの量を、該さらなる生物ネットワークにおける入力として用いることができる。

この手法によって、中でも、生物システムの予期される挙動に対応したネットワーク、反応速度則、並びに反応速度定数及び開始濃度の値の範囲の特定の変形を選択することが可能になる。

本発明のこの実施の形態によって、ネットワーク間の相互作用のシミュレーションが可能になる。たとえば、１つのシミュレートされたネットワークは細胞を表すことができ、第２のシミュレートされたネットワークは第２の細胞を表し、２つの細胞は情報及び／又は生物学的実体を交換することが可能である。明らかに、２つのみでなく、５つ、１０個、数百個、又は数千個の生物ネットワークを付随してシミュレートすることができる。そのようなシミュレーションは、多細胞集合体、組織、器官、生物全体又は相互作用する生物群のシミュレーションである。

生物学的実体は、シミュレーションの時間ステップにおいて毎回、又はより大きな間隔、たとえば１つおき、１０個おき、又は１００個おきの時間ステップにおいて交換することができる。

本発明による方法のさらに好ましい実施の形態では、前記生物学的実体の濃度は、野生型と比較して摂動されている。

本明細書において用いられる際、用語「摂動される」又は「摂動」は、野生型からの偏差を指す。対応する実験的設定において、そのような偏差は、所与の生物学的実体の過小発現若しくは過剰発現、又は突然変異に起因して生じる場合がある。他の予想される摂動は、薬剤又は他の物質の投与によって生じる。

好ましい実施の形態では、摂動系における生物学的実体の濃度は実験的に求められる。モデル化される摂動系は、細胞、組織、器官、生物、又は相互作用する生物のグループとすることができ、摂動は、ノックダウン実験、突然変異、病態、又は薬剤の投与によって引き起こすことができる。

本発明によるノックダウン摂動の好ましい実施の形態では、ノックダウンされている１つ又は複数の生物学的実体の濃度（複数可）は、それらの開始濃度の或る特定のパーセンテージに固定される。１０％及び０％が好ましいパーセンテージである。加えて、ノックダウン実体の濃度を増加又は減少させる反応が無効にされるため、生物学的実体はシミュレーション全体を通じて固定の濃度に留まる。摂動される実体の開始濃度は、対数正規分布から選択されるか、又は遺伝子発現、ＲＴ−ＰＣＲ、定量的プロテオミクス技術及びメタボロミクス技術のような実験データから選択される。

本発明による突然変異摂動の好ましい実施の形態では、生物学的実体に対する突然変異の効果は、文献から既知であるようにモデル化されるか、又は突然変異の効果が未知である場合には、突然変異は、生物情報学技術からの推論を用いてモデル化される。たとえば、サイレント突然変異は、野生型の生物学的実体によって効果的にモデル化され、ミスセンス突然変異は多くの場合に生物学的実体の完全なノックダウン（０％）によってモデル化することができ、既知の機能ドメインを損傷する突然変異は、モデル化された生物学的実体と、損傷されたドメインが相互作用することを意図していた生物学的実体との間の適切なエッジを除去することによってモデル化することができ、構成的に活性化する突然変異は、生物学的実体の不活性形態を活性形態に転換する人工的な非可逆反応（エッジ）を付加することによってモデル化することができ、最後に、酵素生物学的実体の酵素効率を変化させることが知られている突然変異は、反応速度定数に効率の変化の既知の因子を乗算することによってモデル化される。これらの全ての事例において、反応速度定数は実験的に求められるか又は対数正規分布から選択される。

本発明による、病態の摂動の好ましい実施の形態では、疾患に関与するとして知られる突然変異の摂動は、上記のセクションにおいて説明したようにモデル化される。加えて、遺伝子発現、タンパク質及びリンタンパク質濃度、代謝産物及びマイクロＲＮＡレベルにおいて体現される活性病態データは、適切な生物学的実体の初期濃度を説明されたレベルに経験的に設定することによってモデルに直接適用される。

本発明による薬剤投与摂動の好ましい実施の形態では、薬剤が相互作用することが知られている生物学的実体に対する該薬剤の適用の効果は、いくつかの異なる形のうちの１つでモデル化することができる。

ａ）薬剤が、酵素、たとえばキナーゼである１つ又は複数の生物学的実体の活性を阻害することによって作用する場合、その活性は、キナーゼごとに実験的に規定されたＩＣ５０を取り、それを、薬剤の既知のＩＣ５０濃度を薬剤のモデル化された細胞濃度で除算してその結果に５０％を乗算することによりキナーゼの反応速度定数に適用することによってモデル化される。モデル化された細胞濃度は、一般に、適用濃度とみなされる。たとえば、５００ｎＭの薬剤適用をモデル化するために、細胞濃度は一般に５００ｎＭであると想定される。しかしながら、モデル化される薬剤に関する因子、たとえば透過性若しくは溶解性、又は薬剤の薬理作用に影響を及ぼす、モデル化される細胞に関する因子、たとえば上方調節されたＰＧＰ若しくはＭＤＲ−１が知られている場合、モデル化される細胞濃度は、適用される濃度の数分の１に設定することができる。これが行われる場合、実験によるデータに基づくことが好ましい。

ｂ）薬剤が、タンパク質間相互作用を阻害することによって作用する場合、その相互作用、すなわちエッジの反応速度定数は、ａ）におけるような薬剤の既知のＩＣ５０によって変更される。

ｃ）薬剤が、ほとんどの従来の抗腫瘍薬のように非機構的である場合、薬剤の適用の効果は、薬剤に対する感度及び耐性を担うことが知られる、ネットワーク内の生物学的実体を作動させることによってモデル化される。たとえば、シスプラチン及びカルボプラチンのような、白金を用いた化学療法薬は、ＤＮＡをキレートすることによって作用し、この結果、細胞ＤＮＡ修復ネットワークが過度に刺激されることになる。加えて(Iin addition)、細胞は、ＰＧＰ及びＭＤＲ−１を過剰発現することによって、又はＤＮＡ損傷感知及び修復並びにアポトーシスネットワークにおける突然変異を獲得することによって、これらの薬剤に対して耐性を有するようになる。したがって、ＤＮＡ損傷感知及び修復経路における選択された生物学的実体は、構成的に活性化されたものとして系内でモデル化することができる。加えて、非機構的薬剤の効果は、ＲＮＡ及びタンパク質発現パターンにおける変化に基づいて間接的にモデル化することができる。このようにしてモデル化される実体(entitites)は、実際の薬剤の実際の細胞への適用の遺伝子発現又はプロテオメトリック実験から取られることが好ましいが、文献から薬剤応答に関して知られていることからモデル化することもできる。

本発明の方法のさらに好ましい実施の形態では、初期条件は、
（ａ）生物学的実体の実験的に求められた濃度、及び／又は
（ｂ）実験的に求められた突然変異データ
を含む。

本明細書において用いられる際、用語「初期条件」は、シミュレーションの開始時の、上記生物学的実体の性質、数、状態、及び濃度を指す。

本発明は、本発明による生物学的実体の濃度を予測する方法の統計学的有意性を求める、コンピューター実装される方法であって、該統計学的有意性を求める方法は
（ａ）本発明による生物学的実体の濃度を予測する方法を実行すること、
（ｂ）ステップ（ａ）において取得された生物学的実体の濃度と、同じ生物学的実体について実験的に求められた濃度との間の一致度を求めること、
（ｃ）前記生物ネットワークのトポロジーをランダム化すること、
（ｄ）ステップ（ｂ）において取得した前記ランダム化された生物ネットワークに対し、本発明による生物学的実体の濃度を予測する方法を実行すること、
（ｅ）ステップ（ｄ）において取得した生物学的実体の濃度と、同じ生物学的実体について実験的に求められた濃度との一致度を求めること、
（ｆ）ステップ（ｂ）において取得した結果をステップ（ｅ）において取得した結果と比較することであって、ステップ（ｂ）における一致度の方が高いことは、本発明による生物学的実体の濃度を予測する方法が、実験的に求められた濃度を偶然よりも良好に予測可能であることを示す、比較すること、
を含む、方法をさらに提供する。

本発明のこの態様は、本発明による、生物学的実体の濃度を時間の関数として予測する方法を検証する方法に関する。該方法は、確率分布から反応速度定数をランダムに選択する効果を生物ネットワークのトポロジーのランダム化と比較する。好ましくは、生物ネットワークを前記ランダム化することは、前記ネットワークのエッジをスワップすることによって有効にされる。

用語「上記ネットワークのエッジをスワップすること」は、該ネットワークにおける接続の変更を指す。たとえば、生物学的実体の濃度を時間の関数として予測するのに用いられる生物ネットワークにおいて、エッジ１がノードＡ及びノードＢを接続し、エッジ２がノードＣ及びノードＤを接続する場合、スワップされたエッジを有するネットワークの例は、エッジ１がノードＡ及びノードＤを接続し(connections)、エッジ２がノードＢ及びノードＣを接続するネットワークである。

上記で言及した「一致度を求めること」は、利用可能な場合には経時変化を用いて、又は定常状態濃度若しくは平衡濃度のような最終濃度を用いて実行することができる。

好ましい実施の形態では、エッジの１０％、２０％、３０％、４０％、５０％、６０％、７０％、８０％、９０％、又は全てのエッジがスワップされる。

本発明による方法のさらに好ましい実施の形態では、前記実体は生体分子であり、好ましくは遺伝子を含む核酸、タンパク質を含む（ポリ）ペプチド、小分子、並びに生体分子の複合体及び代謝産物から選択される。小分子には、サッカリド、アミノ酸、脂質、ヌクレオチド、ヌクレオシド、並びにそれらの代謝産物及び誘導体が含まれる。

本発明による方法のさらに好ましい実施の形態では、上記モデルは、境界条件、好ましくは生理学的状態を表す境界条件を含む。境界条件は、反応速度定数に対して、かつ／又はグラフの接続性に対して影響を有する場合がある。

好ましい境界条件は、生物学的実体又は刺激の存在又は恒常性を含む。好ましい境界条件は、所与の量の１つ又は複数の薬剤の存在を含むことができる。他の好ましい境界条件は、細胞外シグナル伝達勾配、並びに細胞間の伝達及び生理学的シグナルによって課された境界条件である。

本発明は、先行する請求項のいずれか１項に記載の方法を実行するように構成されたコンピュータープログラムをさらに提供する。

本発明によるプログラムを含むコンピューター可読データキャリアがさらに提供される。本発明による方法を実行する手段を含むか、又は本発明によるプログラムがインストールされたデータ処理装置も提供される。

本発明は、生物ネットワークの部分的に未知のパラメーターを求める、コンピューター実装される方法であって、該パラメーターは、ネットワークトポロジー、反応速度則、反応速度定数、及び／又は開始濃度から選択され、該方法は、観測される特性と予測される特性との差異を最小にすることであって、該予測される特性は、本発明による生物学的実体の濃度を予測する方法によって予測されるような濃度を含む、最小にすることを含む、方法をさらに提供する。たとえば、或る特定の生物学的実体の定常状態濃度又は平衡濃度は、容易に実験的に求めることができる。本発明による生物学的実体の濃度を予測するｉｎ ｓｉｌｉｃｏの方法に関して、これらの生物学的実体の予測される定常状態濃度又は平衡濃度は、多様な程度で、未知のパラメーターに割り当てられる値に依拠する。観測される特性と予測される特性との間の差異を最小にすることによって、未知のパラメーターは、取得される反応速度モデルが、実験とシミュレーションとの間の差異が最小にされた特性を最も良好に再現するモデルであるという意味で最適化される。この最小化は、連続最適化、すなわち、本発明による生物学的実体の濃度を予測する方法の実行中の、上記差異の最小化を伴う場合がある。用語「最適化」は、非線形回帰型手法の使用等による、上記で規定したパラメーター、すなわち、たとえば反応速度定数及び／又は開始濃度の最適化に関する。さらに、最適化は、代替的に又は付加的に、ネットワークのトポロジー及び／又は反応速度則の変更のような不連続ステップを伴う場合がある。好ましくは、予測値と実験値との間の一致は、全ての利用可能な測定値に関して最適化される。観測値と予測値との間の差異、又は観測されるパラメーターと予測されるパラメーターとの間の差異を最小にする汎用計算手段及び方法は当該技術分野において既知である。好ましい実施の形態では、上記ネットワークトポロジー及び／又は上記反応速度則は完全に既知である。

本発明は、１つ又は複数の実験を選択する、コンピューター実装される方法であって、
（ａ）本発明による生物学的実体の濃度を予測する方法を実行することによって複数の実験をｉｎ ｓｉｌｉｃｏで実行することであって、該実行することは未知のパラメーターの異なる選択を用いて実験ごとに反復して行われ、該パラメーターは、ネットワークトポロジー、反応速度則、反応速度定数、及び／又は開始濃度から選択される、実行すること、並びに、
（ｂ）前記複数の実験から、前記異なる選択に依拠して、前記生物学的実体の濃度を予測する方法が予測濃度の最大の分散をもたらす１つ又は複数の実験を選択すること、
を含む、方法をさらに提供する。

用語「実験」は、特記されない限り、実行される実世界の実験に関する。この実験は、生物システムにおいて又は生物システムに対してそれぞれ、１つ若しくは複数の遺伝子のノックダウン、１つ若しくは複数の遺伝子に特異的な１つ若しくは複数のｓｉＲＮＡ（低分子干渉ＲＮＡ）の投与、及び／又は酵素のような１つ若しくは複数のポリペプチドの活性の１つ若しくは複数の修飾因子の投与を実行することを含む実験とすることができる。好ましくは、上記生物システムはｉｎ ｖｉｔｒｏ系である。好ましい修飾因子は、薬剤又は薬剤の開発に適したリード化合物である。

換言すれば、これらの実験を選択する、コンピューター実装される方法であって、モデルの異なる代替形式間を最も効率的に区別し、全ての可能なモデルを用いて実験をモデル化することによって未知のパラメーターを求めるのに最も効果的となる方法が提供される。実験的に求めることができる、濃度のような特性の予測値間の差異が最大である実験は、異なる形式のモデルのうちのいずれが正確である可能性がより高いかに関するほとんどの情報を提供し、したがって、可能であるがおそらく与える情報がより少ない実験のより大きなセットの中から選択されることができる。

１つ又は複数の実験を選択する方法の好ましい実施の形態において、上記ネットワークトポロジー及び／又は上記反応速度則は完全に既知である。好ましくは、上記選択は、上記で詳述したようにランダムな選択である。さらに好ましい実施の形態では、選択された実験（複数可）が実行される。

デビッドソン（Davidson）他［１４、４］から再現される内胚葉系中胚葉ネットワークトポロジーを示す図である。トポロジーは、バイオタペストリ（Biotapestry）［１５］を用いて生成される。 ａｌｘ１−ｍＲＮＡ及びｏｔｘ−ｍＲＮＡのシミュレートされた経時変化を示す図である。ｍＲＮＡ量は任意の単位で与えられ、時間は受精後時間で与えられる。８００回のシミュレーション実行の平均が各時点においてプロットされる。 ｐｍａｒ１発現の摂動によって影響を受ける遺伝子の発現を表す図である。ｍＲＮＡ量は任意の単位で与えられる。検出される摂動の効果は、［２５］に記載されている実験データに従う。 無摂動のｐｍａｒ１−ＭＡＳＯ条件（上）及び無摂動のａｌｘ１−ＭＯＥ条件間の、ｈｅｓｃ−ｍＲＮＡ（左）及びａｌｘ１−ｍＲＮＡ（右）の量を比較する散布図である。ｘ軸は摂動条件のｍＲＮＡ量を表し、ｙ軸は無摂動条件下のｍＲＮＡ量を表す。量は、受精後２５時間（ｈｐｆ）の時点でａ．ｕ．単位で測定される。８００個の異なるパラメーターセットの値がプロットされる。対角線からの偏差は、摂動条件の発現がより高い（対角線の下）か、又は摂動状態の発現がより低い（対角線の上）ことを示す。 アポトーシスの誘導の前（左）及び後（右）の切断カスパーゼ９の量を示す図である。緑は８０％カスパーゼ３ノックダウンに対応し、青は対照に対応し、赤は２つの間の差異に対応する。バー（ｙ軸）の長さは、複数のモンテカルロ−シミュレーションをそれぞれの結果（ｘ軸）と共に示す。初期タンパク質濃度はランダムに選択され、４００個全てのシミュレーションに適用される(apply)。 上枠（Ａ）は、アポトーシスのシミュレートされた誘導の前（左）及び後（右）のカスパーゼ９（未切断）の量を示す図である。下枠（Ｂ）左は、スタウロスポリン(staurosporin)によるアポトーシスの誘導を用いない場合（０ｈ＋）及び用いた場合（６ｈ＋）の、カスパーゼ３ｓｉＲＮＡを９６時間形質移入されたＷＩ３８細胞におけるカスパーゼ９の発現を、アポトーシスの誘導を用いない場合（０ｈ−）及び用いた場合（６ｈ−）の形質移入されていない細胞と比較したウエスタンブロットである。下枠（Ｂ）右は、ＡＩＤＡを用いたウエスタンブロットの定量化を示す図である。値はアクチンを用いて標準化される。 注釈及びネットワークモデル過程を示す図である。（Ａ）は、癌治療の最小ネットワークを含むシグナル伝達経路の一部の図である。（Ｂ）は、ＰｙＢｉｏＳ系内のコンピューター可読反応系におけるこれらの機能的相互作用の説明の図である。 提案されるモンテカルロ手法のフローチャートである。定常状態シミュレーションが正常状態及び摂動状態（薬剤による阻害）に関して実行される。阻害される成分を除く全てのモデル成分に関して、正常状態及び摂動状態が、反応速度パラメーターの同じセット及び初期値を用いて初期化される。その後、モデルは該モデルの定常状態に入るようにシミュレートされる。この手順は、所与のランダム分布からサンプリングされる異なるパラメーターベクトルを用いてｋ回反復される。グラフィック検査に関して、それぞれの対照及び処理シミュレーション実行の定常状態の結果は、モデルの全ての成分についてヒストグラムを用いてプロットされる。さらに、対照及び処理間の結果の分布は、対照／処理値のセットごとにコルモゴロフ−スミルノフ検定によって計算されるＰ値によって定量化される。 全体統計を示す図である。（Ａ）は、阻害実験ごとに有意に変化したモデル成分の数によって、処理の個々の感度を示す図である。（Ｂ）は、ヒストグラムによって、阻害実験の集団（panel）に対するモデル成分の感度を示す図である。このヒストグラムは、この研究において分析される２４個の阻害実験のうちの所与の数（Ｘ軸）によって有意に変更されるモデル成分の数（Ｙ軸）を示す。（Ｃ）２４個の異なる阻害実験の対数比の主成分分析（ＰＣＡ）を示す図である。ＰＣＡはＪ−Ｅｘｐｒｅｓｓ Ｐｒｏ（ベルゲン、Molmine社）を用いて行われた。（７６７個のうちの）２つの第１の主成分によって全分散の６７％が明らかになる。 ＡＫＴシグナル伝達を標的にする特定の阻害実験を示す図である。ＡＫＴによって調節される複合ネットワークを標的にする４つの薬剤／小分子阻害剤（ペリホシン、ウォルトマンニン、メトホルミン、及びラパマイシン）を示す図である。変化は増殖、成長、及びアポトーシスに影響を与える。阻害は、点線（blunted line）によって示される。ＩＲＳはインスリン受容体基質であり、ＡＫＴはタンパク質キナーゼＢであり、Ｐｌ３Ｋはホスファチジルイノシトール３−キナーゼであり、ｍＴＯＲはラパマイシンの哺乳類標的である。 薬剤標的阻害の直接的及び間接的効果を示す図である。選択された対照状態（上側バー）対処理状態（下側バー）の定常状態頻度のヒストグラムを示す。選択結果は、図１５に示される特異的阻害実験の直接的な効果（左枠）及び間接的な効果（右枠）を示す。（Ａ）は、阻害されたＡＫＴ２の、「ホスホ−ＧＳＫ３β」の濃度に対する直接的な効果を示す図である。（Ｂ）は、ＡＫＴ２阻害の間接的な効果が、「Ｓ６Ｋリン酸化」の阻害を示すことの図である。（Ｃ）は、ＰＤＫ１が、「ＰＩＰ３：リン酸化ＰＫＢ」の阻害における直接的な効果を示すことの図である。加えて、（Ｄ）は、「ホスホ−ＧＳＫ３β」の阻害が、治療ドメインにおいて間接的な効果とみなされることを示す図である。（Ｅ）は、ＩＲＳが「ホスホ−ＩＲＳ：Ｐｌ３Ｋ」の形態でＰｌ３Ｋの阻害を通じて直接的な効果を与えることを示す図である。（Ｆ）は、活性化ＩＲＳの阻害の結果、「ＴＳＣ２−Ｐ」リン酸化の阻害が生じることを示す図である。（Ｇ）は、治療ドメインＡＫＴＰＩ３Ｋが、「ホスホ−ＧＳＫ３β」に対する直接的な効果を示すことの図である。（Ｈ）は、治療ドメインＡＫＴＰＩ３Ｋが、「Ｓ６Ｋ１−Ｐ」に対する間接的な効果を示すことの図である。合計２５０回のシミュレーション実行が、阻害実験ごとに行われた。実際には、或る特定の量のシミュレーション実行を超えた実行時間制約に起因して複数のシミュレーション実行が失敗することにより、結果としてのシミュレーション実行回数はこれより小さい（平均２２０回）。 ＡＫＴ阻害を伴う薬剤作用の共発現クラスターを示す図である。異なる阻害実験（列）及びモデル成分（行）に対する、処理状態及び対照状態の平均定常状態値の対数比（底２）を示す図である。赤色は阻害に対するモデル成分の上方調節を示し、緑色は下方調節を示す。対類似度測定としてピアソン相関が用いられ、更新規則として平均連結法が用いられた。Ｊ−Ｅｘｐｒｅｓｓ Ｐｒｏ（ベルゲン、Molmine社）を用いてクラスター化が実行された。 ＡＫＴ（ＲＡＣβセリン）に対して高い感度を示すモデル成分のクラスターを示す図である。この研究に用いられた２９個の異なる単一阻害実験にわたって同様の感度パターンを示すモデル成分のクラスターを示す図である。薬剤標的ごとの感度値（列）が全てのモデル成分（行）について計算された。−２０〜２０の段階の感度値がプロットされる。色の段階で示されるように、高い感度成分は赤又は緑である。Ｊ−Ｅｘｐｒｅｓｓ Ｐｒｏ（ベルゲン、Molmine社）を用いてクラスター化が実行された。以下の実施例は、本発明を説明するが、限定的なものとして解釈されるべきでない。 摂動シミュレーションの定性的結果を示す図である。緑色の細胞は、列の摂動に起因して行内の遺伝子の発現が増加したことを示し、赤は発現が低減したことを示す。細胞は、或る特定の領域に含まれた効果を強調するために網掛けにされる。影響を受ける領域は個々の細胞において言及される（Ｅ：内胚葉、Ｍ：中胚葉、Ｐ：ＰＭＣ、Ｔ：全ての領域の合計）。一時的な発現情報は表から除外された。それによって、遺伝子が何らかの所与の時点において有効にされた場合、その効果は遺伝子一般に割り当てられる。 シミュレーションと実験データとの比較の定性的概観を示す図である。表を短くするために、異なる時点の比較結果が結合されている。符号化は以下の通りである。緑は摂動及び有効にされた遺伝子の全ての実験データが整合していることを示し、緑及び黒のパターンは、所与の摂動及び遺伝子に関して不整合よりも整合の方が多く存在することを示し、灰色は、同程度の数の整合及び不整合が存在することを示し、赤及び黒のパターンは整合よりも不整合の方が多く存在することを示し、赤は不整合のみを示す。最終行は、各列の平均を示す。

実施例１：内胚葉系中胚葉ネットワーク
背景
受精卵からの成体の発生は、基本過程と同様に複雑である。発生は、周囲の細胞からのシグナルによって駆動される個々の細胞の特異化、並びに細胞運動性及び卵割事象を含む。主な調節的入力は多数の細胞によって生成されるが、これらの巨視的効果を生成する微視的事象は、細胞レベルで精密に調節されなくてはならない。このため、発生メカニズムは、シグナル伝達事象及びタンパク質相互作用、並びに遺伝子調節及び細胞間相互作用を含む。これらの発生メカニズム、及び異なる種間のこれらのメカニズムの差異を理解することによって、進化メカニズムを見通すことができる［１］。さらに、発生中のいかなる摂動も、何らかの形で生物において顕在化する可能性が高い。

発生の科学的分析は、ウニを用いて１８９０年代に始まった［２］。ウニは歴史的理由からのみでなく、その興味深い進化位置からもモデル生物である。進化プログラムの基本過程は、哺乳類の発生にも同様のものがあると予期される［３］。

今日、ウニの発生の理解における最も野心的な取り組みは、アメリカムラサキウニ（Ｓ．ｐｕｒ．）における内胚葉系中胚葉ネットワークである［４］。このネットワークは、初期のウニ胚における発生を制御及び調節する複雑な遺伝子調節ネットワーク（ＧＲＮ）、主に内胚葉、中胚葉、及び第一次間充織（ＰＭＣ）領域の特異化の構造を捉えようと試みている。遺伝子間の相互作用は、摂動実験、主に、モルホリノ基置換アンチセンスオリゴヌクレオチド（ＭＡＳＯ）注入が、特定の遺伝子のｍＲＮＡ翻訳［５］、ｍＲＮＡ過剰発現、及びエングレイルドリプレッサードメインとの融合を効果的にノックダウンすることから推測される。このネットワーク（図１）は静的であり、相互作用のタイプ、タイミング、及び強さに関する規則を規定する必要がある。本発明の方法を用いて、内胚葉系中胚葉ネットワークが常微分方程式（ＯＤＥ）を用いてモデル化された。実験データは主に摂動実験［４］に基づき、詳細な研究は限られた数の遺伝子に関してのみ実行された［６、７、８、９、１０、１１、１２、１３］ため、実験データは結果としてのモデルを完全にパラメーター化するには不十分である。

データがわずかであるにも関わらず時間依存的挙動を予測するために、ランダムにサンプリングされたパラメーターを用いてシミュレーションを実行した。パラメーター変動に対する各ネットワーク成分の感度を検出することに加えて、ネットワークが基づく実験データに関する摂動から、パラメーター変動に対してロバストな効果が検出される。これは同じパラメーターセットを用いて、正常条件及び摂動条件の下でモデルをシミュレートすることによって達成される。結果としてのシミュレーション対の比較によって、用いられたパラメーター値とは独立した効果が示される。

比較結果は、モデルトポロジーがどの程度実験データを再現することができるかを示す。ネットワーク内の全ての相互作用がパラメーター変化に対して不変であることを一見して予期することはできない。

モデルのランダム化されたバージョンを用いて検証が行われる。元のモデルを用いて実行したのと同じ分析が該モデルのランダム化されたバージョンに適用される。元のモデルが実験データを再現する程度を、ランダム化されたバージョンが該データを再現する程度と比較することによって、元のモデルの妥当性が推論される。

方法
ネットワーク妥当性の推論
ネットワークトポロジーの妥当性を推論するために、ネットワーク構造を含む入力ファイルがＰｙＢｉｏＳ［２９、３０］においてＯＥＤモデルに転換された。結果としてのモデルは、Ｐｙｔｈｏｎ［３１］で構築され、容易にシミュレートすることができる。

問題となっている異なる胚領域のより現実的なシミュレーションを作成するために、全て互いから独立している、同じ元モデルの３つのコピーを含むモデルが作成された。異なる外部入力を系に適用するため、これらの入力によって異なる領域を判別することができる。ＰｙＢｉｏＳ出力を用いて、モデルのシミュレーションのために、転写調節のためのランダムパラメーターのセットがサンプリングされる。モデルは、７０個の時間ステップにわたってシミュレートされた。ここで、１つのタイムステップは受精後１時間（ｈｐｆ）に対応する。外部で設定された入力を用いて胚領域を判別するため、これらの入力はタイマーとしての役割を果たし、実験データに従って発現の時間フレームを確立する。

ノックダウン実験をモデル化するために、問題となっているｍＲＮＡ生成の速度則を設定することによってＰｙＢｉｏＳシンタックスにおけるモデルが、ｖ_{ｔｒａｎｓｃｒｉｐｔｉｏｎ}からｖ_{ｔｒａｎｓｃｒｉｐｔｉｏｎ}×０．０５に変更された。過剰発現実験の場合、元のモデルにおいて到達した最大発現率が任意の単位で約２／ｈであるため、ｖ_{ｔｒａｎｓｃｒｉｐｔｉｏｎ}はｖ_{ｔｒａｎｓｃｒｉｐｔｉｏｎ}＋２に変更された。ランダムにサンプリングされたパラメーター値を用いると、この値は任意となり、ＭＯＥ実験の分析における誤りの原因となり得る。

上述した方法で変更された任意の数のモデル及び元のモデルは、同じパラメーターセットを用いてシミュレートすることができる。この研究において用いられるパラメーターセットは、σ＝１．５、μ＝０．５の対数正規分布からサンプリングされた。これは、本発明による好ましい分布である。この分布は、一つには、ＯＤＥ系の数値的安定性を妨げ得るあまりに極端なパラメーター値を回避するために選択された。この研究では、シミュレーションのために８００個の異なるパラメーターセットを用いた。

全ての可能なパラメーターセットＳの下で、有効にされた可能性がある遺伝子Ｘのそれぞれにつき１つの特異摂動ｐに関するシミュレーション結果のセットごとに、次のステップが計算される。

Ｔ＝（１４，１９，２５，３３，４５，６６）となるような、［１９］において与えられる時間間隔を反映する時点のリストＴが規定される。

条件ｐ及び無摂動条件Ｃ下の全てのパラメーターセットＳに関する結果が対で整列される。

ｔ∈Ｔごとに、最も近いシミュレーション時点が求められ、双方の条件に関するシミュレーション結果、及びパラメーターセットｓ∈Ｓごとの比ｒ_{ｘ，ｓ，ｔ}＝［ｘ］_{Ｃ，ｓ，ｔ}／［ｘ］ｐ，ｓ，ｔが記憶される。これらの値を用いて、摂動及び遺伝子ごとに、全てのパラメーターセットＳにわたる平均ｒ_{ｘ，Ｓ，ｔ}及び分散ｖａｒ_{ｘ，Ｓ，ｔ}を計算する。これが３つ全ての領域及び３つの領域の和について計算された。

１つの遺伝子ｘに関する値が、視覚定位に関して図１５に散布図で示されている。摂動の効果は、或る特定の閾値を上回っているか又は下回っている発現比の数として測定される。具体的には、或る特定の閾値ｚ_ｕを上回っている比又は或る特定の閾値ｚ_ｌを下回っている比を数える。

１つの所与の時点ｔについて、Ｔ個より多くのＳ_Ｒに関してｒ_{ｘ，ｓ，ｔ}＞ｚ_ｕ（ただしｚ_ｕ≧１）である場合、遺伝子ｘは時点ｔにおいて摂動ｐによって有意に下方調節されるといわれている。他方で、Ｔ個より多くのＳ_Ｒ及び１つの特定のｔに関してｒ_{ｘ，ｓ，ｔ}＜ｚ_ｌ（ただしｚ_ｌ＜１）である場合、ｘは時点ｔにおいて摂動ｐの下で有意に上方調節されるといわれている。この分析において、ｚ_ｕ＝ｚ_ｌ＝１、及び閾値Ｔ＝９０％を設定する。有意でない変化から判別するために、Ｓ_ＲのうちのＴ個について０．８≦ｒ_{ｘ，ｓ，ｔ}≧１．２５である場合の効果は有意でないものとみなした。計算された全ての摂動に関する全てのシミュレーション結果を用いて、それらを利用可能な実験データと比較すると、すなわち、実験結果にマッチするシミュレーション結果における上方調節／下方調節／無差別の割合として、モデルがどの程度実験データを再現するか、すなわちモデルの妥当性を計算することが可能である。異なる領域の独立したシミュレーションから全体的な妥当性を推論するために、領域のうちの少なくとも１つにおける実験データとの整合を、モデルと実験データとの間の有意な整合として数えることを選択した。

シミュレーションのセットが実験データの所与のセットを再現する程度は、ｚ_ｕ、ｚ_ｌ、Ｔの選択及び有意でない効果に関する閾値に依拠する。系の知識が限られていることに起因して、比較的緩い閾値を設定した(set chose)。モデルのサイズがより小さいか、又は詳細で包括的なモデルが利用可能であるサブセットが用いられる場合、これらの閾値を厳しくすることを選択してもよい。

ロバスト性の推論
内胚葉系中胚葉ネットワークモデルの異なる成分のロバスト性の推論は、無摂動モデルのシミュレーション結果のみを用いる。全ての利用可能なパラメーターセットからのシミュレーション結果を比較することによって、遺伝子ごとのシミュレーション結果がパラメーター値に依拠して異なる程度が抽出される。これらの変動の程度は、パラメーター値に対する遺伝子の発現のロバスト性である。

遺伝子の発現の１つの真の経時変化が存在すると想定され、これはサンプリングされたパラメーターセットによって必ずしも発見されない。それにもかかわらず、異なるパラメーターセットに関するシミュレーション結果は、経時変化の個々の時点を分析するとき、真の経時変化の周りに正規分布することがさらに想定される。

この想定によって、異なる時点において、遺伝子ごとにｍＲＮＡ濃度の平均μ及び分散σ^２を計算することが可能になる。平均に対する分散の程度を測定するために、絶対平均値又は絶対分散値とは独立して、時点と遺伝子との間でσ^２／μで比較することができる相対変動ｖａｒ_ｒｅｌ＝を計算することができる。これは基本的には、対応する平均に対する分散を与える。

遺伝子ごとに十分な数の時点を選択するとき、ｖａｒ_ｒｅｌのリストを用いて、遺伝子の発現の全体的なロバスト性を得ることができる。ここでは、遺伝子ごとに最大ｖａｒ_ｒｅｌを選択することによって行われる。結果としての値が実質的に異なり、異なるレベルのロバスト性を示す場合がある。ｖａｒ_ｒｅｌ及びそれらの平均は、ロバストな遺伝子と脆弱な遺伝子との区切りを求める絶対値としては解釈されない。なぜなら、これらの区切りは未知である実際のパラメーター値に大きく依拠するためである。

ネットワークのランダム化
ＰｙＢｉｏＳ形式のＯＤＥモデル。これは、閾値モデルの場合、２つの遺伝子をランダムに選択し、領域ごとに２つのランダムに選択された入力をこれらの遺伝子と交換することによって行われる。このランダム化手順を３０００回反復する２つのランダムネットワーク、及びランダム化ステップを３００００回反復する１つのモデルを作成した。３つの領域間の境界が維持されることに留意されたい。このランダム化アルゴリズムを用いて、ノード及びエッジの数、平均ノード次数及び次数分布のようなネットワークの全体特徴は、個々の配線が変更される一方で保持されることに留意されたい。ランダム化されたネットワークが構築された後、上述したようにその妥当性を判断することができる。この分析において、３０００個のランダム化ステップを用いて２つのランダム化モデルを構築し、それぞれを１００個の異なるパラメーターセットを用いてシミュレートした。より強力なランダム化の効果を推論するために、３００００個のランダム化ステップを用いる１つのモデルも構築し、該モデルを２０個のパラメーターセットを用いてシミュレートした。

結果及び検討
内胚葉系中胚葉ネットワークの動的モデル
内胚葉系中胚葉ネットワーク［４］は、ウニＳ．ｐｕｒにおける内胚葉、中胚葉、及びＰＭＣの分化を駆動する遺伝子間の予測される調節的相互作用を示している。このネットワークの精緻化されたバージョンが基本データ［１９］と共に利用可能である［１４］。このデータに加えて、いくつかの遺伝子の調節的相互作用が詳細に研究された［６、７、８、９、１０、１１、１２、１３］。

摂動及び利用可能な経時変化データ（言及したソースにおけるもの）は一般に、胚全体について測定されるが、定性的データが利用可能であり、或る特定の遺伝子が或る特定の領域のみで発現されるか、又はさらには複雑な時空間発現パターンに従うことが示される［７、６、２０］。この差次的発現は、胚の細胞間の直接的相互作用及び間接的相互作用によって駆動及び実行される［２１、２２］。

各領域は、遺伝子の発現、ＴＦ及びシグナル伝達分子の存在量に関して異なるものの、各細胞が同じ遺伝子情報を含む、すなわち、発生の初期段階においてヒストン修飾が起こらないと想定する。さらに、各領域が、一様な数の細胞から構成される、すなわち、同じ領域内の細胞がＴＦの同じ組み合わせを含み、同じ遺伝子を発現すると想定する。

これらの想定によって、１つの細胞のみをモデル化することにより各領域（内胚葉、中胚葉、及びＰＭＣ）をモデル化することが可能になる。このため、同じメカニズムの３つの複製を含むモデルを構築する。このため、モデル化された細胞間の差次的発現は、細胞間シグナル伝達における差異及び異なる開始条件からのみ生じることができる。内胚葉、中胚葉、及びＰＭＣは胚全体を構成するものではなく、細胞間シグナル伝達の正確なメカニズム及び胚内の拡散率は未知であるため、これらのメカニズムは除外することにした。基本メカニズムを除外したにも関わらずこれらの異なるシグナル伝達事象の効果を含めるために、［１９］に含まれるデータを用いて、モデル自体から出現しない各領域に対する一時入力を模倣する。これらの適切な初期条件、及び細胞外過程を模倣する人工的入力を所与とすると、モデルは、内胚葉細胞、中胚葉細胞、及びＰＭＣ細胞の挙動を再現することが可能であるはずである。

この系に対し、単一の遺伝子のノックダウン（ＫＤ）又は過剰発現の形態で摂動を加えると、各単一の領域に対する効果を求めることができる。効果を結合して実験データにおいて測定されるような全体効果にするための、所与の時点における胚の正確な組成を知らないため、領域のうちの少なくとも１つにおける摂動の効果の検出によって、それに従う胚全体の挙動を示すことを提案する。そのために、このモデルを用いて、既存のトポロジーを検査することができる一方、新たな実験データと、代替的な仮説の検査とを正規に統合することも容易になる。

そのような用途において、常微分方程式（ＯＤＥ）のモデルが広く用いられている。それらのモデルは定常状態の分析及び経時変化の詳細なシミュレーションを可能にする［２３］。それらのモデルは、より単純なモデル化の枠組みよりも詳細な結果を生成するため、ＯＤＥのモデルは、モデル化された系に関する、より詳細な情報も必要とする。

本明細書において上記で説明した反応速度則を用いて、内胚葉系中胚葉ネットワークのトポロジーを用いるＯＥＤモデルを構築した。図１を参照されたい。このモデルは、細胞タイプごとに合わせて５４個の遺伝子及び１４０個の変異種から構成される。ｍＲＮＡ及びタンパク質の転写、翻訳、分解、並びに複合体の会合及び解離を含めて、モデルは、細胞タイプごとに、２８７個のパラメータによって制御される合計２７８個の反応を含む（翻訳、ｍＲＮＡ分解、及びタンパク質分解は１つの偏在性パラメーターによって制御される）。ここでもまた、これらの数はパラメーター推定がこのネットワークに関して現在実現不可能である理由を示している。

図２は、異なる遺伝子及び領域に関してシミュレートされた経時変化を示している。これらの経時変化は実験による経時変化を再現しない場合があるが、これはこの分析の目標では決してなく、実際にパラメーター推定を必要とした。それにもかかわらず、これらの経時変化は、本発明者らのモデルが差次的発現を生成することが可能であることを実証する。ａｌｘ１−ｍＲＮＡ量の経時変化は、ＰＭＣ条件下でのみ発現されることを明確に示している一方、ｏｔｘは３つ全ての領域において発現されるが、各領域において異なる程度で発現されることを明白に示している。

方法の評価
本明細書において説明される方法は、１２個の遺伝子から構成される内胚葉系中胚葉ネットワークのサブモデルに対して検査された［２４］。該サブモデルはわずかに変更を加えられ、該サブモデルに関して、摂動データではなく既知の経時変化データを再現するためのパラメーターが推定された。推定されたパラメーターが真のパラメーターであり、かつネットワークの動的挙動が正確であると想定すると、これを、ランダムにサンプリングされたパラメーターを適用してネットワークのトポロジー特徴を抽出するためのベンチマークとして用いることができる。

このサブネットワークに関して、推定されたパラメーターを用いて取得されたシミュレートされた摂動効果が、ランダムにサンプリングされたパラメーターを用いたシミュレーションから取得されたシミュレートされた摂動効果と比較された。推定されたパラメーターを用いてモデルをシミュレートすることによって達成される効果は、真の実験データの代わりに基準として用いられた。

このサブネットワークの分析によって、ランダムにサンプリングされたパラメーターを有するモデルを用いた摂動効果の検出の受入可能な精度が示される（双方の設定において摂動効果の７３．８％が等しく検出された）。偽陰性（推定されたパラメーターの使用において効果が検出され、サンプリングされたパラメーターセットの使用において効果が検出されない）が８．１％を構成し、偽陽性が１５．７％を構成し、２．１％が２つの設定において効果が逆の符号を有して検出されたことを示す。

これは、サンプリングされたパラメーターセットを用いてＯＥＤモデルのトポロジー特徴を推論することによって、少なくとも小さなネットワークに関して満足のいく結果がもたらされることを示している。多数のシミュレーションに必要とされる計算時間は、選択されるネットワークトポロジー及び反応速度を含む複数の要因に依拠するため、推定されたパラメーターの使用とサンプリングされたパラメーターとの間の正確な比較は、この単一の比較に基づいて達成することができない。それにもかかわらず、パラメーター推定は、パラメーター相互作用及び利用可能なデータに依拠してサンプリングされたパラメーターを用いて効果的にシミュレートすることができる系に関して実現不可能であり得る。

ランダムなパラメーター変動に対する遺伝子発現のロバスト性
方法のセクションにおいて詳述したように、遺伝子の発現のロバスト性は、相対変動（ｖａｒ_ｒｅｌ＝σ^２／μ）を計算することによって、対応する変動を用いて様々な時点において異なるパラメーターセットを用いて発現の平均を比較することによって測定される。無摂動条件下で８００個のパラメーターセットを用いてモデル内の遺伝子ごとにｖａｒ_ｒｅｌを計算した。遺伝子ごとに、全ての時点の最も高いｖａｒ_ｒｅｌ（最も低いロバスト性に関する）が遺伝子のロバスト性を示すものとして用いられる。一般に、ロバスト性は、より初期の測定点において、より高い。これは、問題となっている遺伝子に到達し該遺伝子を生み出すのに時間がかかる、分析される遺伝子の上流の遺伝子の発現における変動に起因する。

表１は、ネットワークにおける各遺伝子のｖａｒ_ｒｅｌに対する概観を与えている。ｖａｒ_ｒｅｌは用いられる８００個のパラメーターセットに関して０．２１〜２５．６９の範囲にわたり、ネットワークにおける遺伝子が、実際は、ランダムなパラメーター変化に対するそのロバスト性において実質的に異なることを示している。

表１：ロバスト性によってソートされた、異なる遺伝子のロバスト性の概観。各遺伝子のロバスト性に関連付けられるスコアは、方法セクションにおいて説明されたように相対変動（ｖａｒ_ｒｅｌ）である。８００個の異なるパラメーターセットからの結果が用いられる。最小スコア（最大のロバスト性を示す）に関してソートする以外はクラスター化又はグループ化はこの表に適用されない。第３の列は遺伝子の入次数（入来する相互作用の数）を与え、第４の列はノード出次数（遺伝子が有する調節的相互作用の数）を示す。

ノード入次数とロバスト性との間の相関（発現のロバスト性が調節因子が増えると共に増加することを示す）を検出することもできないし、ノード出次数とロバスト性との間の相関（重要な調節的遺伝子が他の遺伝子よりもロバストであることを示す）を検出することもできなかった。ネットワークにおける遺伝子の位置（初期内中胚葉、初期ＰＭＣ、後期内胚葉及び中胚葉、後期ＰＭＣ）とその(it's)ロバスト性との間の相関性を見出すこともできなかった。

これは、現在のアーキテクチャがランダム摂動に対するロバスト性においていかなる特定の遺伝子のセットを優遇することがないか、又はネットワークが、同様のロバスト性を用いて表現される遺伝子のセットを検出するには小さすぎることを示す。

異なる摂動のシミュレーション
内胚葉系中胚葉ネットワークの妥当性を推論することを可能にするために、異なる摂動実験をシミュレートする必要があった。ノックダウン実験及び過剰発現（ＯＥ）実験のモデルを効率的に作成するために、転写の速度則、たとえばｖ_{ｔｒａｎｓｃｒｉｐｔｉｏｎ}を、ノックダウンの場合にはｖ_{ｔｒａｎｓｃｒｉｐｔｉｏｎ}・０．０５に、過剰発現の場合はｖ_{ｔｒａｎｓｃｒｉｐｔｉｏｎ}＋２に変更した。

エングレイルドリプレッサードメインを用いた遺伝子の融合の結果発生する効果をシミュレートしなかった。なぜなら、これは全ての影響を及ぼされる速度則を注意深く調整することが必要とするためである。

１０００個のパラメーターセットを用いて、元のモデル及び摂動モデルの双方をシミュレートした。非常に硬いＯＤＥを生成したサンプリングされたパラメーターセットから生じる数値の問題に起因して、８０１個のシミュレーション実行のみが終了した。

たとえば、ａｌｘ１、ｅｔｓ１、ｔｂｒ、及びｈｅｓｃのようないくつかの遺伝子に対するｐｍａｒ１摂動の効果を提示する。異なる摂動条件下のｍＲＮＡ量の平均を経時変化図３としてプロットする。効果を視覚化する別の方法が図４に与えられている。ここでは、或る特定の時点における１つのｍＲＮＡの無摂動及び摂動条件下の量が、すべてのパラメーターセットに関してプロットされる。視覚的効果は、ｈｅｓｃ発現に対するＰｍａｒ１の阻害的役割を示すと共に、ａｌｘ１、ｅｔｓ１、及びｔｂｒに対するＨｅｓＣの阻害的役割を示す、［２５］に説明される所見に従う。本発明者らのデータは、この阻害が取り除かれる場合に、ａｌｘ１、ｅｔｓ１、及びｔｂｒがモデルに従ってより高い率で発現されることをさらに示している。

全てのシミュレートされた摂動実験を視覚化したものが図１４に与えられている。結果は、摂動された遺伝子の標的を示す摂動によって分析することができるか、遺伝子の調節因子を示す有効にされた遺伝子によって分析することができる。列方向の分析によって、実際に異なる胚領域に関連付けられる同時調節された遺伝子のグループが明確に同定される。たとえば、ａｌｘ１、ｄｒｉ、ｅｔｓ１、ｈｎｆ６、及びｐｍａｒ１のノックダウンはＰＭＣ条件下のみにおいて検出可能な効果を有する。調節的効果の一致パターンを有する行は、同様の領域／調節因子に関連付けられる遺伝子のグループを示す。図１の左下部分からの針状のマトリックス遺伝子のグループは、異なる遺伝子の摂動に対するそれらの感度を考慮すると、右下部分の二次間充織細胞遺伝子のグループと大幅に異なる。ＧａｔａＥ及びＨｏｘは概して反対の効果を有し、これはＨｏｘがｇａｔａｅ発現に対して有する阻害的役割によって説明することができる。ｈｏｘ−ＫＤの甚大な効果を、ネットワークトポロジーにおけるその比較的限定された役割と比較すると、多数の検出された効果がｇａｔａｅ発現の阻害に起因することが明らかである。ネットワークトポロジーから明らかであるが、これはシミュレーション結果のみから判別することがむしろ困難である。シミュレーション結果のみを考慮すると、Ｈｏｘ及びＧａｔａＥは共に全ての有効にされた遺伝子を並行に調節していると誤解される場合がある。

内胚葉系中胚葉ネットワークトポロジーの妥当性
内胚葉系中胚葉ネットワーク（図１）の妥当性は、シミュレーション結果において正確に再現される実験結果の比に関して評価される。計算分析の結果は定性的結果しか生じないため、実験データは定量的データから定性的データに変換されることを必要とした。既存の実験データのうちのいくつかは、データが不明瞭であるために除外された。たとえば、ｇａｔａｅ−ＭＡＳＯに応答したｆｏｘｂの発現は、−２．８／−２．４／−３．４、−２．４／ＮＳ、＋３．７である（ただし、スラッシュは異なる実験からの測定値を示し、カンマは反復測定を示し、ＮＳは有意でないと考えられる値を示す）［１９］。このため、妥当性は不明瞭でない実験結果のみを用いて計算される。元の内胚葉系中胚葉ネットワークのモデル、並びに［１９］に示されるＭＡＳＯ及び過剰発現実験に関連するモデルを用いて、元のモデル及び摂動モデルにおける遺伝子の発現が比較される。本発明者らのモデルは、内胚葉特異的濃度、中胚葉特異的濃度、及びＰＭＣ特異的濃度を含み、胚の正確な組成は未知であるため、実験データとシミュレーション結果のうちの任意のものとの間の合致するものを考察して、再現された摂動を示す。

実験データに従う信頼可能な計算に必要な複数のパラメーターセットを試験した。モデルは多数のパラメーターを含むため、全ての可能なパラメーターセットを含むパラメーター空間は巨大である。驚くべきことに、２０個ほどの少ないパラメーターセットであっても、８００個のパラメーターセットを用いて得られた結果と同様な結果を生じるのに十分であった。これらの発見は、表３に要約され、本発明者らの分析において検出されるネットワークのトポロジーの、すなわちパラメーター不変の特徴に対し示唆を与える。

実験データと比較した定量的結果にわたる概観が表２に与えられ、定性的概観が図１５に与えられる。いくつかの摂動の効果が非常に良好に再現されている（ｇａｔａｅ−ＫＤ、ａｌｘ１−ＫＤ、及びｂｒａ−ＫＤ）。

表２：異なる摂動実験に関するシミュレーション結果と実験データとの間の比較の概要。実験データとシミュレーション結果との間の合致するものの数が、列１において与えられ、全ての可能な合致するものの割合としての合致するものの数が第２の列に与えられている。表の底部に全ての摂動の平均が与えられている。

図１５は、これらの観測を裏付けるが、多くの場合に実験データに従って摂動に反応する遺伝子（ｐｋｓ、ｎｒｌ、ｆｖｍｏ、ａｌｘ１、及びｂｒａ）、及び、ｓｍ５０、ｓｍ２７、及びｆｉｃｏｌｉｎのような、その予期されない挙動が、異なるパラメーターセットにより変化する多数のアップストリーム相互作用によって生じる可能性がある、実験データに従って反応することがほとんどない遺伝子をさらに示す。ｆｏｘｂ、ｆｏｘａ、及びｅｖｅのような、重要な調節的遺伝子である他の遺伝子も、非常に多くの場合に実験データの再現に失敗する。遺伝子の双方のセットが実験データを不十分に再生するが、第２のセット内の遺伝子は、上述したものとは対照的に重要な調節的役割を有するため、これらの遺伝子をより差し迫って精緻さを欠いていると見なす。

表３：用いられるパラメーターセット数及び実験データに対する全体的な一致度
要約すると、結果はｅｔｓ１及びｓｏｘｂ１を伴う調節的相互作用の詳細な分析の必要性を指示し、一方で他の遺伝子の調節も同様にチェックする必要がある（ｆｏｘｂ、ｆｏｘａ、及びｅｖｅ）。

全体として、モデルは実験データの４２％を再生する。これはモデルが実験データとの一致度を増大するために精緻化を必要とすることを示す。この精緻化は、より多くの実験データに大きく依拠しなくてはならない。なぜなら、分析された遺伝子に対するモデル化されたＭＡＳＯ摂動の５９２０個の可能な効果のうちの小さな部分しか実験データに関連付けられないためである。

ランダムネットワークとの比較
方法セクションにおいて説明したように、元のＯＤＥモデルにおけるエッジをランダムにスワップすることによって、比較可能な特徴を有する２つのランダムなネットワークが構築された。

ランダム化されたネットワークは、サンプリングされたパラメーターセットを用いて、対照条件及び摂動条件の下でシミュレートされた。これらのパラメーターセットは元のモデルと同様であるが、前のセクション及び表３において説明された所見に起因して、ランダム化されたモデルの場合、元のモデルの場合の８００個ではなく１００個のパラメーターセットのみを用いた点が異なる。ランダム化されたモデルは、それらの一時入力においてのみ異なる３つの同一のサブモデルも含んでいた。ここで分析された２つのランダム化されたモデルは、実験データの２０．１５％及び２３．５％のみを再現することが可能であった。３０００個のエッジの代わりに３００００個のエッジを同様の結果と切り換えた１つのランダム化されたモデルも検査した。

これは元のモデルの内胚葉部分の範囲内であるが、いずれのランダム化されたモデルも、完全なモデルのＰＭＣ部分ほど有意にこの一致度を超えない（表４を参照のこと）。また、８個のランダム化されたネットワークのうちのいかなる３つのから成る組み合わせも、実験データとの全体一致度が２５％を超えることはなかった。したがって、ランダム化されたネットワークの計算される一致度は、モデルの一般的なトポロジー特徴（サイズ及び接続性等）、実験データ、及び一時入力のみに依拠すると想定する。これは、特定の一時入力を用いてここで説明されたモデルのような比較可能な一般的な特徴のモデルと、ここで用いられた実験データとの間の一致度が、偶然のみで２５％未満になると予期することができることを示す。

表４：異なる胚領域に関するシミュレーション結果を用いて、実験データとの達成された一致度が示されている。この表に関して８００個全てのパラメーターセットが用いられている。

このため、偶然により予期される実験データとの一致度は、本発明者らのモデルを用いて観測された一致度の約半分であり、ランダムネットワークと比較して内胚葉系中胚葉ネットワークの妥当性が有意に改善されることを示す。

本発明の方法を用いてウニ内胚葉系中胚葉ネットワークにおいて得られる結果が、刊行物クーン他（２００９）においてさらに詳細に説明されている。クーン他の刊行物（２００９）、及び本文脈では特にクーン他（２００９）の「結果(Reults)及び検討」と題するセクションが、参照によりその全体が援用される。

実施例２
アポトーシス
アポトーシスネットワークのｉｎ ｓｉｌｉｃｏ表現は、９７個の微分方程式及び１１３個（全て未知）の反応速度パラメーターを含む。予測される濃度は、ノックダウン実験から得られた実験データと比較された。カスパーゼＣ３はｓｉＲＮＡを用いてＷｉ３８細胞においてノックダウンされた。シミュレーションにおいてカスパーゼ３の実験的なノックダウンを反映するために、シミュレーションにおけるカスパーゼ３の初期濃度は、対照状況においてカスパーゼ３の濃度の２０％に設定された。いずれの事例においても、すなわち対照及びノックダウンにおいて、アポトーシス経路の時間依存的挙動の４００個のシミュレーションが実行された。結果が以下の表５に示される。

表５：モデルが平衡状態に到達した後の、アポトーシス正常細胞及びカスパーゼ３ノックダウン（２０％）細胞におけるアポトーシス経路のタンパク質及びタンパク質複合体の濃度。値は４００個のシミュレーションの平均に対応する。

カスパーゼ９及び切断カスパーゼ９の実験結果との比較。実験の設定において、カスパーゼ３のノックダウンによって、切断カスパーゼ９の量の有意な減少、及び（切断されていない）カスパーゼ９のわずかな増加が生じた。これらの結果は、本発明による方法を用いて行われた予測と良好に一致する。図５及び図６に示すように、カスパーゼ９から切断されたもの(from)の濃度が減少する一方、カスパーゼ９はそれに付随して増加する。

実施例３
大型癌ネットワークに対するノックダウン効果の計算予測の生成
癌治療の効果を予測するための最小ネットワーク
癌治療のための最小相互作用ネットワークの開発及び注釈における第１のステップとして、癌に関連する遺伝子及びそれらの機能相互作用に関する大規模な情報を利用してきた。癌はおそらく、複数の遺伝子及び経路を伴う最も複雑な疾病のうちの１つであり（ビルド（Bild）他（２００６）、ハナハン（Hanahan）及びワインバーグ（Weinberg）（２０００）、ワインバーグ（２００７））、たとえばアポトーシスの回避及び反成長シグナルに対する不感受性につながる、細胞生理学における激しい機能変化の症状と見なされる。これらの機能変化は、癌発病及び進行に伴う主要な分子及び経路と関連する。ほとんどの癌研究は、癌遺伝子及び腫瘍抑制遺伝子の突然変異に起因するこれらの経路の異常な活性の結果に集中してきた（キンズラー（Kinzler）及びフォーゲルシュタイン（Vogelstein）（１９９６））。細胞増殖及びアポトーシスの調節に重要なのは、細胞シグナル伝達ネットワーク、すなわち、広大なクロストークを示す複合ネットワークによる、成長シグナル及び死シグナルの認識及び統合である。経路間の正のフィードバックループが、不活性から永久に活性化された状態への遷移を引き起こす場合があり、これは連続した細胞増殖へとつながり、このため癌の発症に寄与する（キム（Kim）他（２００７））。

表６：癌における標的にされた治療薬剤。細胞表面受容体又は開始された経路の下流成分を標的にする異なる抗癌剤の選択。阻害実験は、阻害された単一の又は複数のモデル成分に関する。これらの成分は表９において説明される。

現在の抗癌剤の大部分は、少なくとも名目的に、単一の標的に向けられるが、多くの場合、一緒に多くのオフターゲット効果を発生する。これらの薬剤のほとんどは、想定される標的及び関連する経路に対して直接的な阻害性効果を有する（表６）。たとえば、重要な標的経路はＰＩ３Ｋ／ＡＫＴシグナル伝達経路である。なぜなら、該経路は細胞の成長及び存続の多くの局面にとって重大であるためである。該経路はゲノム異常によって影響を及ぼされ、これは多くの癌において経路の活性化へとつながる（ヘネシー（Hennessy）他（２００５））。したがって、いくつかの癌実体、たとえば前立腺癌に関する前臨床及び臨床試験においてＡＫＴ阻害剤ペリホシンが用いられた（ヴァン アンメルセン（Van Ummersen）他（２００４））。この薬剤は、おそらくＡＫＴのＰＨドメインと相互作用することによって膜への局在を防ぐ。ウォルトマンニン及びＬＹ２９４００２は、ｉｎ ｖｉｔｒｏ及びin vivoでの抗腫瘍活性を示す（ヘネシー他（２００５））。メトホルミン塩酸塩は、筋肉細胞及び脂肪細胞上のインスリン受容体の数及び／又は親和性を増大させ、受容体及び後受容体結合部位におけるインスリンへの感度を増大させる（ＮＣＩ薬剤辞書http://www.cancer.gov）。ラパマイシン（WyethAyerst社ラパミューン）はＡＫＴの下流で機能するｍＴＯＲ（ラパマイシンの哺乳類標的）の特異的阻害剤である（ヘイ（Hay）及びソネンバーグ（Sonenberg）（２００４））。ｍＴＯＲ阻害剤は、臨床試験において、乳癌及び他の固形腫瘍を有する患者に関して試験されている（チャン（Chan）他（２００５）、イダルゴ（Hidalgo）及びロウィンスキー（Rowinsky）（２０００）、ナガタ（Nagata）他（２００４））。細胞において、エベロリムスはイムノフィリンＦＫ結合タンパク質−１２（ＦＫＢＰ−１２）に結合し、ｍＴＯＲ、すなわち主要な調節キナーゼに結合してその活性化を阻害する免疫抑制複合体を生成する。ｍＴＯＲ阻害は、下方調節されたＰＴＥＮによって生じるトラスツズマブ耐性を克服する試みにおいて研究されている。テムシロリムス（Wyeth社Ｔｏｒｉｓｅｌ）は、キナーゼｍＴＯＲの阻害剤であり、Ｓ６Ｋ１のリン酸化を阻止し（フェーヴル（Faivre）他（２００６））、進行性腎細胞癌の治療に用いられる。ソラフェニブ（ネクサバール）は、ＲＡＦ−キナーゼ、ＶＥＧＦＲ−２、ＰＤＧＦＲ、ＦＬＴ−２、及びｃ−ＫＩＴに対する経口マルチキナーゼ阻害剤であり（ストランバーグ（Strumberg）（２００５））、血管形成及び腫瘍増殖を標的にする。これはインターフェロンα又はインターロイキン２に基づく治療に耐性を示す進行性腎細胞癌又は腎臓癌を有する患者の治療に適している。ＭＥＫは癌細胞の成長及び存続に伴うＭＡＰＫ経路の重要な構成要素である。ＰＤ−３２５９０１は、この経路を遮断し癌細胞を破壊するように設計された新たな薬剤である。ＰＤ−０３３２９９１は、サイクリン依存性キナーゼ、特にＣＤＫ４を選択的に阻害し、網膜芽細胞腫（Ｒｂ）タンパク質のリン酸化を阻害することができる。Ｒｂリン酸化の阻害は、Ｒｂ陽性の腫瘍細胞が細胞周期のＳ期に入るのを防ぎ、その結果ＤＮＡ複製が抑制され、腫瘍細胞増殖が減少する。ＡＥＧ３５１５６はＸ連鎖アポトーシス阻害タンパク質（ＸＩＡＰ）、すなわち多くの細胞において過剰発現されるアポトーシスの阻害剤の細胞発現を選択的に阻止する。この化学物質は、細胞毒性薬剤と相乗的に作用して腫瘍細胞におけるＸＩＡＰの総レベルを低減し、アポトーシスに対する腫瘍細胞耐性を克服する。別の化合物、ＦＪ９は、Ｆｒｉｚｚｅｄ−７ Ｗｎｔ受容体とＤｉｓｈｅｖｅｌｌｅｄのＰＤＺドメインとの間の相互作用を妨害し、古典的（canonical）Ｗｎｔシグナル伝達を下方調節し、腫瘍細胞成長を抑制する（フジイ（Fujii）他（２００７））。エンザスタウリン塩酸塩は、ＡＴＰ結合部位に結合して、プロテインキナーゼＣβを選択的に阻害する。

細胞の成長及び生存に不可欠な重要なシグナル伝達経路は、突然変異、増幅、転位を含むゲノム異常に起因してヒト癌において頻繁に活性化される。ＲＡＳ−ＲＡＦ−ＭＥＫ−ＥＲＫ、ＰＩ３Ｋ−ＡＫＴ−ｍＴＯＲ、又はＪＡＫ−ＳＴＡＴシグナル伝達経路のような成長及び生存経路を対象とする、ますます多くの合理的に設計された小分子阻害剤が、今や癌の治療の臨床試験に入ってきている（ヘネシー他（２００５）、マッキュブリー（McCubrey）他（２００８）、ヴァン アンメルセン他（２００４））。それにもかかわらず、多くの阻害剤は、事前に知られていない「オフターゲット」によって生じる予期されない毒性に起因して、又は薬剤標的自体が脱制御に感受性のある複数の機能相互作用に関与するため、臨床試験において失敗する。加えて、標的にされた薬剤の臨床的失敗は、予期されないフィードバックループの存在、代替シグナル伝達経路の補完的上方調節、又は薬剤耐性突然変異によっても生じ、それらの全ては標的阻害の効果を回避し、腫瘍細胞を生存及び増殖させる（ブルヘルト他（２００５））。

したがって、これらの実施例が示すように、予測モデルは、複数の標的の複雑度及び経路間のクロストークに対処するために、関連する機能相互作用の多く（又は理想的には全て）を含むべきである。そのようなモデルは、患者の予期されない副作用又は不感受性を明らかにすることによって、新規な標的薬剤の開発の重要なサポートを提供することができる。しかしながら、これまで、癌過程の計算的モデリングは主に、ＲＡＦ（キム他（２００７））、ＡＫＴ（アラウホ（Araujo）他（２００７））、又はＷＮＴシグナル伝達（キム他（２００７））のような個々の下位経路に焦点を当ててきた。しかしながら、いくつかの分離した経路をより大きなフレームワークに統合することは、経路間のクロストークも捕捉し、薬剤作用の予測に不可欠である場合がある。この観点において、本発明者らは、次の段落で説明される癌に関連する過程の「最小」相互作用ネットワークに対する薬剤作用の効果をシミュレートすることによって計算による予測力を検査した。薬剤、それらの分子標的又は標的のセット、及びそれらが機能する細胞相互作用ネットワークに関する情報を凝集し、次のステップは阻害性効果をコンピューターにおいて翻訳することである。これは、薬剤をそれらの意図される薬剤標的タンパク質に関連させ、薬剤標的に関連付けられる単一の又は複数のモデル成分の阻害によって該薬剤標的の阻害の効果をシミュレートすること(to simulate)によって行われる（表６）。

モデル集団及び注釈ワークフローの自動化
がん治療のための最小予測相互作用ネットワークを構成する主成分の選択は、主に生物学的機能に基づく。薬剤効果の計算モデル化は、経路スキーマ（図７Ａ）を、モデル成分の濃度及びこれらの成分が関与する反応の反応速度に関する情報を保有することができるコンピューターモデルに変換することを必要とする。この過程は、適切なコンピューターオブジェクトの設計、反応の実装、これらの反応に対する反応速度則の割り当て、及びモデル分析を含む（図７Ｂ）。

モデル集団、シミュレーション、及び分析をサポートする計算ツールがシステムズ生物学の基礎を構築する（ビリング（Wierling）他（２００７））。いくつかのモデル化ツール、たとえばＣｅｌｌＤｅｓｉｇｎｅｒ（フナハシ（Funahashi）他（２００３））、Ｅ−Ｃｅｌｌ（トミタ（Tomita）他（１９９９））、Ｇｅｐａｓｉ（メンデス（Mendes）（１９９３）、メンデス（１９９７））、及びその後継機Ｃｏｐａｓｉ（ホープス（Hoops）他（２００６））、ＰｒｏＭｏＴ／Ｄｉｖａ（ギンケル（Ginkel）他（２００３））、Ｖｉｒｔｕａｌ Ｃｅｌｌ（レーブ（Loew）及びシャーフ（Schaff）（２００１）、スレプチェンコ（Slepchenko）他（２００３））、及びＳｙｓｔｅｍｓ Ｂｉｏｌｏｇｙ Ｗｏｒｋｂｅｎｃｈにおいて統合されるツール（ハッカ（Hucka）他（２００２））が最近提案されている。これらのシステムのほとんどは、わずかな反応を用いて小さなモデルを手動で分析するように設計されている。しかしながら、関連する予測モデルの開発は、多数の成分、反応、及び反応速度パラメーターに関する情報を必要とする。このため、モデル化過程の自動化は、いくつかのステップにおいて、たとえば反応ネットワークの注釈、モデルの流れの可視化、微分方程式の数値積分、及びモデル分析において、大きなネットワークを扱う際、サポートを必要とする。

この研究において、モデル化及びシミュレーションシステムＰｙＢｉｏＳ（ビリング他（２００７））を用いてそのような手法を実証する。ＰｙＢｉｏＳ（http://pybios.molgen.mpg.de）は、パスウェイデータベースへのインターフェースを提供することによってモデルの自動生成をサポートし、これによって関連する反応系への高速で自動化されたアクセスが可能になる。癌に関連する反応ネットワークに対する既存の知識の多くは、ＰｙＢｉｏＳシステムに直接取り込むことを可能である、ＢｉｏＣｙｃ（カープ（Karp）他（２００５））、ＫＥＧＧ（カネヒサ（Kanehisa）他（２００６））、及びＲｅａｃｔｏｍｅ（ジョシートープ（Joshi-Tope）他（２００５）、ワストリック（Vastrik）他（２００７））等のパスウェイデータベースに凝集される。

本発明者らの分析は、３つの部分、すなわち癌に関連する経路を含むモデルの確立、たとえば薬剤の効果に起因するモデル成分の阻害の導入、及びモデル分析から構成される。薬剤標的阻害の効果をシミュレートする予測ネットワークの原型は、文献からの情報、及び癌の特質に焦点をあてた２０個の異なるシグナル伝達経路に関する公共のパスウェイデータベースからの情報を組み合わせることによって編集された。現在のところ、ネットワークは１９１３個の個々の反応を有する７６７個の成分を含む（概要は表７を参照されたい）。

表７：注釈付きヒト癌ネットワークに含まれるモデル成分のリスト
^＊Ｃａｎｃｅｒ Ｇｅｎｅ Ｃｅｎｓｕｓ：ｈｔｔｐ：／／ｗｗｗ．ｓａｎｇｅｒ．ａｃ．ｕｋ／ｇｅｎｅｔｉｃｓ／ＣＧＰ／Ｃｅｎｓｕｓ／
^＊＊ラス（Russ）及びランペル（Lampel）（２００５）に基づく創薬可能遺伝子

経路注釈は大方依然として手動で実行されるが、モデル化システムへの静的モデルのアップロードを記憶し容易にする自動化された注釈のためのツールが利用可能である。本発明者らの原型ネットワークの経路注釈は、シグナル伝達経路の情報を収集及び記憶する過程を自動化するＲｅａｃｔｏｍｅ Ｃｕｒａｔｏｒ Ｔｏｏｌを用いて実行された（ｈｔｔｐ：／／ｗｗｗ．ｒｅａｃｔｏｍｅ．ｏｒｇ）。全体ネットワークは、２０個の異なる経路から構成される。それらの経路は、増殖因子（ＥＧＦ、ＮＧＦ、ＩＧＦ−１、ＴＧＦ−β）、細胞増殖（Ｗｎｔ、Ｒｂ、ノッチ受容体、ヘッジホッグ）、サイトカイン（インターロイキン２、ＳＴＡＴ−ＪＡＫ）、炎症（トール様受容体）、アポトーシス（ＴＮＦα、ＦＡＳ、ＴＲＡＩＬ）、及び代謝調節（Ｇタンパク質共役型受容体）のような刺激によって活性化されるシグナル伝達カスケードを構成する。反応はＰｙＢｉｏＳモデル化システムに取り込まれ、対応する数学的ＯＤＥモデルが反応系から自動的に作成された。

モンテカルロ手法：不確実性に直面したモデル化
モデル化は、モデルサイズと予測精度とのトレードオフである。高い精度を有するモデルは、或る程度少ない数のモデル成分に基づいて非常に詳細な（of large detail）計算予測を生成する。しかしながら、これらの予測は、多くの場合にin vivo条件下で関連パラメーターを測定することが困難であることによって支障をきたし、これは関与する異なる経路間のクロストークを無視することによって悪化する。しかしながら、パラメーター適合戦略は、任意のそのような手法の一般的な難点を被るものの、不正確なモデルであっても、適合を生成するために十分なパラメーターを変形することができれば、一般的に極めて良好に適合することができる。特に、医学関連のモデルは、モデル精度に関する結果をもたらす多数のモデル成分を伴う可能性が高い。この観点において提案される戦略はこの目的のために設計される。モデルに関与する反応は、合成反応、複合体及び生成物形成、並びに不可逆質量作用反応速度

によって説明される分解反応のような少ない数の異なる反応タイプから構成される。ただし、ｋは反応速度定数であり、Ｓ_ｉはｉ番目の基質の濃度である。

可逆反応はそれぞれ不可逆質量作用反応速度を用いる別個の前方向反応及び後方向反応によって説明される。複雑な形成及び解離の場合、１００［ａ．ｕ．］の解離定数を有する可逆質量作用反応速度が適用された。合成及び減衰反応が適宜導入された。適用されたモデルの速度則が図８に要約される。

表８：異なるタイプの反応速度
本発明者らのモデル化手法において、同じネットワーク異なる状態、たとえば薬剤標的の異なるセットを不活性化したか又は不活性化していない生物ネットワーク間の差異の予測に焦点を当て、薬剤処理の効果の単純化されたモデルを表す。パラメーターのうちの多くにおける不確実性を補償するために、パラメーターベクトルの成分が適切な確率分布から選択され、利用可能な知識が反映される。ここで説明されるシミュレーション実行のセットでは特に、未知の反応速度パラメーターが対数正規分布から平均μ＝２．５及び標準偏差σ＝０．５でサンプリングされた。各パラメーターベクトル及び入力値の各ベクトルを用いて、通常対照状態及び阻害実験に対応する「処理された」状態を、該２つの状態の初期の差異が１つのパラメーター又は(parameteror)モデルパラメーターの小さなセットのみの変化に起因するようにモデル化した。異なる阻害実験のシミュレーションに関して、表９に列挙されているモデル成分が、標的タンパク質の完全な除去に対応する０．０［ａ．ｕ．］の固定濃度を用いて初期化された。全ての成分の定常状態レベルが到達されるまで、異なるパラメーターベクトルを用いて対照及び処理シミュレーションが２５０回反復された（図８）。

「処理された」状態と「処理されていない」状態との間のモデル挙動の差異が、個々のモデル成分ごとの定常状態濃度を比較することによって分析された。成功したシミュレーション実行にわたるモデル成分の２つの状態の最終定常状態値における差異を、コルモゴロフ−スミルノフ検定を用いて統計学的有意性に関して分析し（コノバー（Conover）（１９７１））、特別な治療によって影響を受けるモデル成分を特定する。

薬剤標的阻害の系統的分析
提案されるフレームワークの検査として、表９に列挙されている異なる阻害実験に従って、無摂動状態におけるモデル成分の挙動を、処理された状態と比較した。薬剤標的のコンピューターによる阻害によって、いくつかのレベルで結果を予測する機会が与えられる。一方で、これは、定常状態濃度に従って統計学的に有意な変化のセットを計算することによって処理の感度の全体評価を可能にする。他方で、これは直接効果（治療における所望の効果）及び間接効果（治療の潜在的な望ましくない副作用）のような、経路成分における特別な変化の同定を可能にする。

表９：抗癌剤の効果をシミュレートする阻害実験（表６を比較されたい）及びそれぞれの実験において阻害された関連付けられるモデル成分

大域分析によって薬剤作用の差異をモニタリングする
図９は、全体統計を要約している。異なる阻害実験に関する複数の有意な変化によって表される摂動感度は或る程度変動する（図９Ａ）。いくつかの阻害ドメイン、たとえば、ＡＫＴ２酵素（モデル成分ＰＩＰ３：リン酸化ＰＫＢ）の活性体の阻害が、異なる阻害実験において単一の標的（ＩＲＳ）又は複数の標的（ｍＴＯＲ、ＩＲＳ、ＡＫＴ、及びＰＩ３Ｋ）の阻害によって６０個より多くの異なるモデル成分に影響を及ぼすのに対し、他の阻害ドメインは非常に特異的であり、たとえばＳＴＡＴ阻害は７６７個のモデル成分のうちの１０個未満に影響を及ぼす。他方で、標的感度は非常に高く、ほとんどのモデル成分は阻害性効果に対してロバストである（図９Ｂ）。モデル成分の最も大きな部分（７６７個のうちの５２０個）は、阻害実験のいずれによっても影響を及ぼされない。最も有意な変化は、単一の阻害実験（７３個）又は５つ未満の異なる阻害実験（２つ以上５個以下について１３８個）においてのみ観測される。モデル成分の小さな部分（３５個）が、阻害実験のうちの多数（８個以上）における効果を示している。多数の薬剤標的阻害によって影響を及ぼされる成分は、シグナル伝達経路において中心的役割を果たすもの、たとえばＧＳＫ３β、及びそのリン酸化形式、若しくはＷｎｔシグナル伝達経路からのＡＰＣ及びアキシンに結合するＧＳＫβ、又はｍＴＯＲシグナル伝達の中心成分である。一般に、選択された薬剤及び阻害ドメイン（表９）は主成分分析（図９Ｃ）によって示されるように３つの異なるグループに入る。特に興味深いのはＩＲＳ、ＡＫＴ２、ＰＩ３Ｋ、ＰＤＫ１のグループ及びそれらの組み合わせである。なぜなら、これらの阻害実験はほとんどのモデル成分に非常に影響を及ぼすためである。

ＡＫＴシグナル伝達を標的にすることによって、直接的な下流効果及び間接的な下流効果が明らかにする
大域分析を補完するものとして、モデル化手法は主要処理に関して詳細な予測を生成する。いくつかの阻害実験は、リン酸化ＧＳＫ３β（ホスホ−ＧＳＫ３β）、及びＴＳＣ１：阻害ＴＳＣ２−１−Ｐ複合体のようなＡＫＴシグナル伝達の直接的な不活性化又は間接的な不活性化を標的にする。これらの成分は、たとえば増殖及びアポトーシスとして癌関連細胞の読み出しに影響を与える。これらの介入点及び読出し点の概略図が図１０に示されている。それぞれのモデル成分の阻害によって、薬剤副作用の経路クロストークの重大性を示す、直接的な効果及び間接的な効果が明らかになる。

図１１は、阻害実験の直接的な効果（左枠）又は間接的な効果（右枠）を示す選択された結果を示している。たとえば、リン酸化ＧＳＫβ（ホスホ−ＧＳＫ３β）の定常状態濃度の低減は、ＡＫＴ２阻害単体（ＡＫＴ２）の直接的効果として、又はＰＩ３Ｋ阻害と組み合わせて（ＡＫＴＰＩ３Ｋ）見て取ることができる（図１１Ａ、図１１Ｇ）。リン酸化ＧＳＫβ（ホスホ−ＧＳＫ３β）の濃度における直接的な低減は、活性ＡＫＴ２及びＰＩ３Ｋ（ＡＫＴ２、ＰＩ３Ｋ）の阻害に起因する。リン酸化ＧＳＫβ（ホスホ−ＧＳＫ３β）の非リン酸化形式は、Ｗｎｔシグナル伝達において重要な役割を果たす。ＰＤＫ阻害（図１１Ｃ）は、ＡＫＴ２のリン酸化において直接的な効果を有し、結果としてＰＩＰ３：リン酸化−ＰＫＢ複合体の下方調節をもたらす。ＰＩＰ３がセリン／トレオニンキナーゼＰＤＫ１及びＡＫＴ２を原形質膜に入れることが既知である。ここでＡＫＴ２はＰＤＫ１媒介リン酸化によって活性化される。ＩＲＳ阻害実験において、ＰＩ３Ｋの阻害は、インスリン受容体基質（ＩＲＳ）のリン酸化（活性）形態の阻害に起因して直接的な効果と見なされ、ＰＩ３Ｋ（図１１Ｅ）の主要活性剤は、それに続く、ホスホ−ＩＲＳ：ＰＩ３Ｋ複合体の定常状態濃度の低減をもたらす。

補間的に、モデル化戦略は多くの間接的な効果を同定する。Ｓ６Ｋ１は、小さな（４０−Ｓ）リボソームサブユニットの成分であり、このサブユニットがリボソーム形成及びしたがってタンパク質合成に関与することを可能にする。Ｓ６Ｋ１のリン酸化体が阻害実験ＡＫＴ及びＡＫＴＰ３Ｋにおいて下方調節されることがわかった（図１１Ｂ、図１１Ｈ）。Ｓ６Ｋリン酸化の阻害は、ＡＫＴシグナル伝達経路におけるＰＩ３Ｋカスケードの下流成分及びｍＴＯＲシグナル伝達に対する下流効果に起因する。ＧＳＫ３βのリン酸化（図１１Ｄ）は、ＰＤＫ（阻害実験ＰＤＫ１）によって間接的に制御される。なぜなら、ＡＫＴの阻害に起因して、ＰＤＫの阻害がホスホ−ＧＳＫ３βの下方調節につながるためである。さらなる間接的な効果として、活性化されたＩＲＳの阻害によって（阻害実験ＩＲＳ）、ＴＳＣ２のリン酸化の低減につながる（図１１Ｆ）。

阻害実験の予測効果の共発現クラスターは文献による知識との一致を示す
多くのモデル成分が、阻害実験の集団全体にわたって同様の発現パターンを示しており、薬剤作用のこれらの共発現クラスターのうちのいくつかを、文献からの以前のデータによって説明することができる。図１２は、ＡＫＴを伴う阻害実験によって影響を及ぼされるモデル成分の特殊な例を示している。

ＡＫＴ活性化は、プレクストリン相同ドメイン（ＰＨＤ）のホスファチジルイノシトール−３，４，５−三リン酸（ＰＩＰ３）への結合、及びその後に続く調節アミノ酸セリン４７３及びトレオニン３０８のリン酸化によって開始される膜への局在によって駆動される（ヴィヴァンコ（Vivanco）及びソイヤーズ（Sawyers）（２００２））。ＡＫＴと原形質膜との病的会合はＡＫＴを癌につなげる共通関連性（common thread）である。ＡＫＴの活性体の阻害に基づく薬剤効果（阻害実験ＡＫＴ２、ＤｖｌＡＫＴ２、ＡＫＴＰＩ３Ｋ、ＩＲＳＡＫＴＰＩ３Ｋ、ｍＴＯＲＩＲＳＡＫＴＰＩ３Ｋ、表９を参照されたい）に関して、下流成分の下方調節を同定することができる（図１２）。さらに、阻害実験ＰＩ３Ｋ及びＰＤＫ１は、ＡＫＴ阻害に基づく阻害実験と比較して同様の阻害された成分を示している。この観測は、ＡＫＴがＰＩ３Ｋ及びＰＤＫ１の効果の一次下流仲介物質として同定される文献データ（ヘネシー他（２００５））と一致する。

ＩＲＳ阻害実験における成分のうちのいくつかは、ＡＫＴ２と同様の挙動を示す。しかしながら、ＩＲＳ阻害実験は、３つの特異成分（ＲＡＣ−βセリン／トレオニンプロテインキナーゼ［ＡＫＴ２］、そのＰＫＢ調節因子［ＰＫＢ：ＰＫＢ調節因子、ＡＫＴ２はＰＫＢの別名である］との複合体、及びＰＩ３Ｋ）の上方調節によって特徴付けられるのに対し、ホスホ−ＩＲＳ：ＰＩ３Ｋ複合体は下方調節される。成分ｅＥＦ２Ｋ−Ｐ、ｅｌＦ４Ｇ−Ｐ、リン酸化４０Ｓ小リボソームタンパク質、ｅｌＦ４Ｂ−Ｐ、ＴＳＣ１：阻害されたＴＳＣ２−１−Ｐ、Ｓ６Ｋ１−Ｐ、活性化されたｍＴＯＲＣ１、阻害されたＴＳＣ２−１−Ｐ、及びリン酸化ＡＫＴ複合体の阻害は、ＡＫＴ阻害に起因する。ＡＫＴ１によるＴＳＣ１：ＴＳＣ２複合体のリン酸化の結果、ＴＳＣ１：阻害されたＴＳＣ２−１−Ｐ複合体、及びプロテオソーム経路を通じたその後の分解がもたらされる。ＰＤＥ３Ｂリン酸化の下方調節はＡＫＴ阻害に起因し、ＰＩ３Ｋ阻害実験においても注目される。これによって、ウォルトマンニンを通じたＰＩ３Ｋ阻害が、ＰＤＥ３Ｂのリン酸化及び活性化、並びにインスリンの抗脂肪分解効果を阻害するという報告（ラーン（Rahn）他（１９９４））が裏付けられる。さらに、ＰＩ３Ｋによる活性化に起因して、ｍＴＯＲがＳ６Ｋをリン酸化及び活性化し、ｍＲＮＡの翻訳を増加させるため、Ｓ６Ｋのリン酸化におけるＰＩ３Ｋの影響が知られている（ローデス（Rhodes）及びホワイト（White）（２００２）、サルティエル（Saltiel）及びカーン（Kahn）（２００１））。

感度分析を用いたネットワーク依存性のモニタリング
感度分析を用いて、異なるモデル成分間の関係及び相互作用を理解する。これらの相互作用を定量化するためのいくつかの方法が提案された。たとえば、代謝制御分析（ＭＣＡ）は基本系における無限小の摂動に対するモデル成分の感度を研究する（クリップ（Klipp）他（２００５））。チョー（Cho）他（２００３）において、マルチパラメトリック感度分析（ＭＰＳＡ）が紹介された。この方法によって、数学モデルがモンテカルロシミュレーションに基づく手法に関して非常に感受性があるパラメーターの同定が可能になる。チアン（Jiang）他（２００７）において、時間依存定量的制御値を計算することによって系の局所挙動が分析された。時間におけるパラメーター感度軌跡を用いて、反応速度パラメーターにおける無限小の変化に対する代謝産物の感度が求められた。

特定のモデル成分の１００％の阻害に加えて、より小さい入力変動におけるモデル成分の挙動を調査するために、感度分析を実行した。一度に２９個のモデル成分を用いて小さな阻害（１０％）が実行され、以前のセクションで説明した異なる阻害実験において分析された。これらのモデル成分のそれぞれについて、その定常状態値の９０％に固定された選択された成分を除いた対照状態の定常状態値を用いて阻害実験のＯＤＥモデルを初期化した。続いて、１０％阻害に関する定常状態の変化が計算された。１０％低減の場合の各モデル成分の定量化された感度が、

によって与えられる。ただし、ｉ＝１，．．．，ｎは全ての成分にわたるインデックスであり、ｊ＝１，．．．，ｍは阻害された標的のインデックスである。したがって、

は制御実行におけるインデックスｉを有する成分の定常状態濃度を表し、

はインデックスｊを有する薬剤標的の１０％阻害モデルにおける成分ｉの定常状態濃度を表す。

結果としての感度値のクラスター化によって、２９個の阻害実験のセットにわたって同様の感度パターンを示すモデル成分のグループが明らかになる。

図１３は、ＡＫＴ２（別名ＲＡＣβセリン）に対し高い感度を示すモデル成分の選択されたクラスターを示している。非活性ＡＫＴ２をわずかに低減することによって活性ＡＫＴ２（ＰＩＰ３：リン酸化ＰＫＢ複合体）、並びにその後の標的、すなわちＴＳＣ１：阻害されたＴＳＣ２及びホスホ−ＧＳＫ３βにおける有意な変化がもたらされる。感度分析によって、ＧＳＫ３βのリン酸化体（ホスホ−ＧＳＫ３β）も、様々な他の薬剤標的（たとえばＳＲＣ、ＧＳＫ３、ＰＩＰ３複合体、ＰＤＫ１、ＰＩ３Ｋ）における小さな変化に感受性があることが明らかになっている。

実施例１で引用される引用文献

実施例３で引用される引用文献

実施例４
薬剤作用の反応速度データを用いた計算予測の生成
本発明によるモンテカルロ戦略は、系／モデルのそれぞれの酵素（たとえばキナーゼ、ホスファターゼ）に対し、薬剤の反応速度結合定数のような薬剤作用に関するサポートデータを用いることによって精緻化することができる。たとえば、カラマン（Karaman）他（２００８）において提供されるような結合定数における反応速度データを以下のように考察することができる。

［ＥＩ］が酵素／阻害剤複合体の濃度であり、［Ｅ］及び［Ｉ］がそれぞれ遊離型酵素及び阻害剤の濃度(the concentrations of the concentrations of)である所与の反応に関して、定常状態における阻害された酵素と遊離型酵素との間の比は、解離定数によって以下のように与えられる。

ただし、ｋ_１及びｋ_−１は、それぞれ解離反応の反応速度定数及び会合反応の反応速度定数(the kinetic constants of the dissociation respective association reaction)であり、Ｋ_Ｄは解離定数である。

ＧＥが遊離型酵素及び阻害された酵素の総キナーゼ濃度である場合、遊離型酵素の濃度ｙを

によって計算することができる。

阻害剤の濃度が各キナーゼの濃度と比較して高い場合、この式は、薬剤によって阻害され得る任意のキナーゼに適用することができる。

阻害剤の全体濃度が細胞内のキナーゼの濃度と同じ範囲内にある場合、個々のキナーゼに結合される阻害剤の量も考慮に入れなくてはならない。

Claims (8)

1. 演算部と記録部とを備えるコンピュータにおいて、癌に対する薬剤の影響を分析するためにコンピュータが行う方法であって、
   （ａ）前記演算部が、癌を表す生物ネットワークの反応速度モデルを受け取り、前記記録部に記録する工程(該モデルは複数個のモデル成分を含む)；
   （ｂ）前記演算部が、前記反応速度モデルに基づいて、該生物ネットワークにおける薬剤の影響をコンピュータ・シミュレートして、ｉｎ ｓｉｌｉｃｏで摂動実験を実施する工程；
   （ｃ）前記演算部が、癌における薬剤の効果を同定するための該摂動実験の結果を、前記記録部に記録する工程、
   を含み、
   摂動は突然変異又は薬剤投与によって引き起こされるモデル成分の過剰発現又は過小発現を意味し、シミュレーションの結果は癌に対する最適な治療を決定するために使用される、方法。
2. 上記モデル成分が、生化学反応、化学反応、生物学的パスウェイ、タンパク質、突然変異遺伝子、創薬につながる遺伝子及び複合体を含む群から選択される、請求項１に記載の方法。
3. 前記癌を表す生物ネットワークが、細胞などの自然に発生する対応物の分析に基づいたものである、請求項２に記載の方法。
4. 前記反応速度モデルは、前記コンピュータにおいて、
   前記演算部が、ノード及びエッジを含むネットワークトポロジーの選択を受け取ること(該トポロジーのノードは遺伝子、ペプチド、タンパク質、小分子、複合体及び代謝産物を含む生物学的実体の群から選択されるモデル成分を表し、該トポロジーのエッジは該モデル成分間の相互作用を表す)、
   前記演算部が、該相互作用についての反応速度則及び反応速度定数の割り当てを受け取ること、並びに
   前記演算部が、該生物学的実体についての濃度の割り当てを受け取ること、
   によって生成され、該反応速度モデルは自然に発生する対応物を反映している、請求項１に記載の方法。
5. 前記ｉｎ ｓｉｌｉｃｏでの摂動実験が、生物ネットワークにおける選択された薬剤の効果を表している、請求項１に記載の方法。
6. 前記ｉｎ ｓｉｌｉｃｏでの摂動実験が、モデル成分における選択された薬剤の効果を表している、請求項５に記載の方法。
7. 前記ｉｎ ｓｉｌｉｃｏでの摂動実験が、ノックダウン実験である、請求項１に記載の方法。
8. 癌に対する薬剤の影響を分析するための装置であって、
   （ａ）癌を表す生物ネットワークの反応速度モデルを受け取り、記録する手段(該モデルは複数個のモデル成分を含む)；
   （ｂ）前記反応速度モデルに基づいて、該生物ネットワークにおける薬剤の影響をコンピュータ・シミュレートして、ｉｎ ｓｉｌｉｃｏで摂動実験を実施する手段：
   （ｃ）癌における薬剤の効果を同定するための該摂動実験の結果を記録する手段、
   を備え、
   摂動は突然変異又は薬剤投与によって引き起こされるモデル成分の過剰発現又は過小発現を意味し、シミュレーションの結果は癌に対する最適な治療を決定するために使用される、装置。

Images