説明

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

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

【発明の詳細な説明】
【技術分野】
【0001】
本発明は、生物ネットワークの反応速度モデルを生成する、コンピューター実装される方法であって、
(a)ネットワークトポロジーを選択することであって、該トポロジーのノードは生物学的実体を表し、該トポロジーのエッジは該実体間の相互作用を表す、ネットワークトポロジーを選択すること、
(b)前記相互作用に反応速度則及び反応速度定数を割り当てること、並びに
(c)前記生物学的実体に開始濃度を割り当てること、
を含み、
(i)前記反応速度定数の一部分、及びそれとは独立して前記開始濃度の一部分が実験データであり、
(ii)前記反応速度定数の残りの部分、及びそれとは独立して前記開始濃度の残りの部分がランダムに選択される、方法に関する。
【0002】
本明細書において、特許出願及び製造者のマニュアルを含む多数の文書が引用される。これらの文書の開示は、本発明の特許性に関連するとは見なされないものの、参照によりその全体が本明細書に援用される。より詳細には、参照される全ての文書は、各個々の文書が参照によりその全体が本文書に援用されると具体的に個々に指示された場合と同程度、参照により援用される。
【背景技術】
【0003】
in silico生物学は、複雑度が高まり続ける生物システムのシミュレーションを可能にしている。これは、疾患に関連した過程への詳細な洞察をもたらし、薬剤開発過程を加速させると共に、治療に対する反応を予測するという可能性を保持している。しかしながら、生物学的過程の詳細な動的モデルを作成する際、実験データが限られていることが主なボトルネックとなる。生物学的経路、ネットワーク、及び過程に関与する分子は多くの場合かなりの程度まで知られているが、これは多くの事例において、これらの分子間の相互作用を定量的に説明する反応速度定数に当てはまらない。未知のパラメーターを推定することは、十分であるとはみなされない。なぜなら、推定は任意であると共にモデルにバイアスをもたらす場合があるためである。
【発明の概要】
【発明が解決しようとする課題】
【0004】
本発明の根底にある技術的課題は、特に実験データがモデルをパラメーター化するのに十分でない事例において、生物システムの時間依存的挙動を予測する手段及び方法を提供することである。
【課題を解決するための手段】
【0005】
したがって、本発明は、生物ネットワークの反応速度モデルを生成する、コンピューター実装される方法であって、
(a)ネットワークトポロジーを選択することであって、該トポロジーのノードは生物学的実体を表し、該トポロジーのエッジは該実体間の相互作用を表す、ネットワークトポロジーを選択すること、
(b)前記相互作用に反応速度則及び反応速度定数を割り当てること、並びに
(c)前記生物学的実体に開始濃度を割り当てること、
を含み、
(i)前記反応速度定数の一部分、及びそれとは独立して前記開始濃度の一部分が実験データであり、
(ii)前記反応速度定数の残りの部分、及びそれとは独立して前記開始濃度の残りの部分がランダムに選択される、方法に関する。
【0006】
本明細書において用いられる際の用語「モデル」は、生物システムのin silico表現を指す。「反応速度モデル」は、生物システムの時間依存的挙動を記述することが可能なモデルである。時間依存的挙動を予測するのに必要な要素は、生物システムの構成要素の転換を含む、生物システムの構成要素間の相互作用を支配する反応速度則及び関連付けられた反応速度定数を含む。これらの構成要素は、本明細書において、「生物学的実体」とも呼ばれる。用語「生物学的実体」は、生物システムにおいて発生し得る任意の分子を含む。好ましい生物学的実体は、下記でさらに詳述する生体分子である。構成要素である生物学的実体は、生物システムのin silico表現として(an)モデルを表現する。本発明によるモデルは、生物学的実体の開始濃度をさらに含む。反応速度則、反応速度定数、及び開始濃度は合わせて、上記生物ネットワークの時間依存的挙動の予測を可能にする。用語「割り当てる」は、シミュレーションの開始時に、或る特定の特性又は数値を定めるか又は設定することを指す。反応速度則及び反応速度定数は、シミュレーション中に変化しないことが好ましいが、シミュレーション中に想定されるような生物学的実体の濃度は、それぞれの開始濃度と異なる場合があることは自明である。本発明が関係する生物システムは生物ネットワークである。好ましい生物ネットワークは、全ての細胞内相互作用ネットワークを含む。それらの例は、シグナル伝達ネットワーク、転写制御ネットワーク、代謝ネットワーク、感覚及び恒常性ネットワーク、分解ネットワーク、調節ネットワーク、並びにそれらの組み合わせである。
【0007】
好ましい生物ネットワークは、たとえば受容体リガンド作用、密着結合のような透過性接触、宿主病原体相互作用、及び細胞又は生物間の任意の他の相互作用によって媒介された全ての細胞間相互作用ネットワークも含む。それらの例は、細胞の成長及び分化ネットワーク、血管新生ネットワーク、創傷治癒ネットワーク、炎症及び免疫反応ネットワーク、並びに、結果として機能性組織、器官、生物、及び生物の集合体が生じる、細胞間又は生物間の相互作用の複合ネットワークである。
【0008】
また、ネットワークは、「グラフ」と呼ばれ、「グラフ」として表されてもよい。ネットワーク又はグラフの例が、本明細書に添付された図1に示されている。より具体的には、当該技術分野において既知であるように、ネットワーク又はグラフはノード及びエッジを含む。ノード及びエッジは合わせてネットワークのトポロジーを形成する。上記ネットワークのノードは、上述した生物学的実体のin silicoの対応物であり、上記ネットワークのエッジは上述した実体間の相互作用のin silicoの対応物である。本明細書において用いられる際の用語「相互作用」は、任意の種類の相互作用、特に上記相互作用に関与する生物学的実体の量又は濃度に影響を与える場合がある相互作用を指す。より具体的には、用語「相互作用」は、おそらく1つ又は複数のさらなる生物学的実体の影響下における、1つ又は複数の所与の生物学的実体の、1つ又は複数の異なる生物学的実体への転換を含む。他の好ましい相互作用は、たとえば1つ又は複数の他の生物学的実体の作用、存在、又は欠如の結果としての、1つ又は複数の生物学的実体の量又は濃度の減少又は増加を含む。さらに別の好ましい相互作用は、2つ以上の生物学的実体からの複合体の形成である。
【0009】
換言すれば、本発明による相互作用は、反応を伴うか又は引き起こす。本発明による反応は質量作用反応速度を用いてモデル化することができるが、一般に、任意の他の適切な反応速度則に従うことができる。当該技術分野において既知であるように、質量作用反応速度は、所与の反応に関与する生物学的実体の濃度及び反応速度定数に依拠する。詳細に関しては下記を参照されたい。
【0010】
従来技術によれば、生物システムの反応速度のモデル化には、関与する全ての生物学的実体又は反応物の全ての反応速度則、反応速度定数、及び開始濃度を知る必要がある。本発明によれば、反応速度モデルであって、反応速度定数の一部のみ及び開始濃度の一部のみが実験データであるか又は実験データから導出され、反応速度定数の残りの部分及びそれとは独立して開始濃度の残りの部分がランダムに、すなわち実験データに依存することなく選択される、反応速度モデルが提供される。この手法は、モンテカルロ手法とも呼ばれる。
【0011】
反応速度定数を求めるのに適した実験データは、相互作用に関与する生物学的実体の経時変化を含む。しかしながら、上述したように、そのような情報は生物ネットワークの多くの相互作用の場合に利用可能でないか生成するのが困難である。開始濃度に関する実験データは、シミュレートされる生物ネットワークの自然に発生する対応物内、すなわち、たとえば細胞内で測定を実行することによって取得することができる。
【0012】
既知の反応速度定数及び既知の開始濃度の使用は互いに独立している。したがって、本発明は、全ての開始濃度が実験データであり、上記反応速度定数の一部がランダムに選択される実施形態、並びに全ての反応速度定数が実験データであり、上記反応速度定数の一部がランダムに選択される実施形態を含む。
【0013】
驚くべきことに、本発明者らは、反応速度定数の一部分又はさらには全てが知られていない事例において、生物ネットワークが反応速度定数の特定の選択に関してロバストであることに気がついた。これは、特に、生物ネットワークによって想定される定常状態及び平衡にあてはまる。特に、反応速度定数を規定する実験データが存在しないときであっても、生物ネットワークの時間依存的挙動によって、偶然による予測をはるかに超える、或る程度の再現可能な予測が生成されることがわかった。これに関して、下記でさらに詳述される、in silicoモデルの統計学的有意性を判断する方法に言及する。同じ観測及び考察を、実験的に知られた開始濃度の部分的欠如又は完全な欠如にも準用する。
【0014】
本発明は、さらに、生物ネットワークにおいて、生物学的実体の濃度を時間の関数として予測する方法であって、
主要な実施の形態に記載の方法によって生物ネットワークのモデルを生成すること、及び
(d)連立微分方程式を解くことであって、該微分方程式は前記生物学的実体の前記濃度の時間依存性を規定する、連立微分方程式を解くこと、
を含み、それによって前記濃度を時間の関数として取得する、方法に関する。
【0015】
生物学的実体又は反応物の量又は濃度の時間依存性を、微分方程式によって、より具体的には常微分方程式(ODE)によって表すことは、当該技術分野において既知である。
【0016】
より具体的には、1つの生物学的実体の量又は濃度を制御する相互作用は一般に、以下のようにモデル化することができる。2つの正入力A及びA、並びに1つの阻害性入力Iを有する生物学的実体を考える。いずれの正入力も量又は濃度を増加させるのに十分である(A∨A)、一方、1つの阻害性入力の活動は標的の量又は濃度を減少させるのに十分である((A∨A)∧¬I)。いくつかの事例では、ネットワークトポロジー及び実験データに基づいて他の形で相互作用を規定することが好ましい場合がある。この表記において相互作用はブール形式で定式化される。
【0017】
ODEモデルの場合、質量作用反応速度が用いられる。分解、輸送、シグナル伝達、及び複合体の会合/解離、並びに翻訳過程及び翻訳後過程のような相互作用を、以下の形式の質量作用反応速度を用いてモデル化した。
【0018】
【数1】

【0019】
上述した反応のグループごとに、1つの偏在性パラメーター(ubiquitous parameter)kreactiongroupを用いた。
【0020】
ブールモデルの定式化を速度則に移すために、[28]から導出される手法を用いることができる。特に、以下の基本モジュール、すなわち、実体Gに対する正入力Aの場合、
【0021】
【数2】

【0022】
及び負入力Iの場合、
【0023】
【数3】

【0024】
を用いることができる。定数c及びkは、各実体の調節的役割の個々の特徴を表す。ここで、kは阻害剤がないときの活性強度に対応するのに対し、cは活性に有意な変化を生成するのに必要な入力の量又は濃度を定める。基本モジュールを組み合わせて、調節的相互作用を含む複合相互作用を定式化することができる。これは、乗算(ブール演算子ANDに対応する)又は加算(ブール演算子ORに対応する)を用いて行われる。
【0025】
このため、モデルが遺伝子発現に関連すると想定すると、2つの活性化剤及び1つの阻害剤を用いた上述の例は、vtranslation=(kA1・[A1]
【0026】
【数4】

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

【0030】
ここで、kdeg・[x]は分解項である。S∈0,1であるため、2つの他の項のうちの一方のみが0に等しくない。このため、Θ及びSを変更することによって、外部入力が時点Θにおいて活性化される(S=1の場合)か、又は時点Θにおいて不活性化される(S=0の場合)。事象を用いてS及びΘを変更することによって、外部入力の活動を制御することができる。
【0031】
本発明の根底にある数学的概念及び方法論は、クーン(Kuhn)他の刊行物(2009)においても詳細に説明されている。クーン他の刊行物(2009)、この文脈では特に、クーン他(2009)の「方法」と題するセクションが参照によりすべて援用される。
【0032】
本発明の方法の好ましい実施の形態では、前記反応速度定数の前記残りの部分は確率分布から選択され、それとは独立して、前記開始濃度の前記残りの部分は確率分布から選択される。それぞれの確率分布は同じであるか又は異なっている場合がある。
【0033】
用語「確率分布」又は「確率密度関数」は、当該技術分野において既知である。この用語は、特定の事象、この事例では反応速度定数の特定の値又は開始濃度の特定の値を、その発生確率に関連付ける。反応速度定数ごとに、それぞれに利用可能な知識の度合いを反映して選択された確率分布から反応速度定数がランダムにサンプリングされる。それとは独立して、開始濃度に関して選択された確率分布から開始濃度がランダムにサンプリングされる。2つ以上の生物学的実体の開始濃度が選択される事例では、開始濃度を選択するのに同じ確率分布又は異なる確率分布を用いることができ、ここでもまた利用可能な部分的(又は完全な)知識を反映する。さらに説明すると、事前知識が一切ないとき、確率分布は好ましくは全ての未知のパラメーターについて同じである。用語「パラメーター」は反応速度定数及び開始濃度を含む。限られた知識若しくは大まかな知識、又は利用可能な類似の相互作用からの知識しか存在しない事例では、スケーリング係数、及び/又はそのようなタイプの情報を反映するための分布関数の幅の変更によって、その種の知識を考慮に入れることができる。さらに、下記でさらに詳述するように、反応速度則はランダムに、好ましくはここでもまた利用可能な知識に依拠する確率を用いて選択することができる。一般的に言えば、本発明のこの前進法はパラメーター推定の既知の過程とは異なる。パラメーター推定では、モデルパラメーターは、観測に適合する最適パラメーターセットを求める目的で、数学的方法によって推定される。提案される手法では、パラメーターは反復的にランダムに選択され、生成された観測値の有意性が統計学的方法を用いて判定される。
【0034】
好ましい実施の形態では、前記分布は、対数正規分布である。対数正規分布において、反応速度定数の対数は正規分布であり、すなわちガウス分布に従う。確率分布の妥当性は、分野における用途及び事前知識に依拠する。さらなる適切な確率分布には、一様分布、指数分布、ポアソン分布、二項分布、コーシー分布、ベータ分布、及びガウス分布が含まれる。
【0035】
本発明による、生物学的実体の濃度を予測する方法のさらなる好ましい実施の形態では、前記反応速度定数の前記残りの部分をランダムに選択し、前記開始濃度の前記残りの部分をランダムに選択し、その後前記連立微分方程式を解くことが反復的に実行される。
【0036】
この実施の形態は、反応速度定数の異なるセットに対する生物ネットワークの応答の評価を可能にする。反応速度定数の異なるセットは、それらの少なくとも一部が次々とランダムに選択される。驚くべきことに、生物ネットワークの反応速度挙動は、限られた数のパラメーターに依拠し、かつ、ほとんどの反応速度定数の異なるランダムな選択は、生物ネットワークの時間依存的挙動に或る特定の影響を及ぼすものの、上記時間依存的挙動、特に定常状態又は平衡状態を根本的に変えることがないことがわかり、本明細書に添付した実施例において立証されている。
【0037】
実験的に求められた反応速度定数が完全に欠如しているとき、かつ/又は実験的に求められた開始濃度が完全に欠如しているときであっても、十分なモンテカルロモデル化によって再現可能な予測を達成することができ、これらの予測を既存の知識によって検証することができることは、特に驚くべきである。換言すれば、実験データである反応速度定数の数がゼロであり、生物ネットワークにおける相互作用を規定する全ての反応速度定数がランダムに選択されるという結果になることが意図的に予想される。特定の理論に拘束されることなく、生物ネットワークの特定のトポロジー、すなわち、ノード及びノード間の接続が(ノード間の接続に関連付けられる反応速度定数を無視して)、潜在的な結果の空間を、反復的なモンテカルロモデル化によって系の時間依存的結果の空間の効果的な調査が可能になる範囲で制限すると見なされる。
【0038】
好ましい実施の形態では、反応速度定数のランダム選択及び常微分方程式系の解決は、複数回行われ、理想的には、利用可能な計算ハードウェアによって制限される実行可能な回数行われる。本発明者らの経験では、現在の計算上の制限に基づいて、10回〜1000回の実行が最適であることがわかっているが、10,000回以上実行すると、精度が増大するという形でさらなる増分値が生成され続けることもわかっている。
【0039】
好ましい実施の形態では、全ての利用可能な実験的に導出された反応速度定数及び開始濃度がシミュレーションにおいて用いられ、本発明者らの経験では概して、既知の反応速度定数が一般的にシミュレーションの精度を改善する。しかしながら、別の好ましい実施の形態では、或る特定の既知の反応速度定数は、確率分布から選択された反応速度定数と選択的に置き換えられる。これは、実験的に導出された定数の精度に関して懸念がある事例において、又は実験的に導出された定数が系が定常状態に達するのを妨げる場合に行うことができる。
【0040】
さらに好ましい実施の形態において、特定の生物システムのシミュレーションのための生物学的実体の反応速度定数及び開始濃度が、同様の系を用いた以前のシミュレーションから導出される。同様の系とは、分析中の系に近い系、たとえば同じであるが特異摂動を伴う生物システムを指す。用語「摂動」は、下記でさらに定義される。定常状態に達するために、以前のシミュレーションが複数のパラメーターセットを用いて実行される。次に、生物学的知識に従って、これらの定常状態のサブセットが選択される。生物学的知識とは、モデル及び既知のパラメーター値によって再現することができる既知のモデル予測を指す。次に、パラメーターセットのこれらのサブセットがシミュレーションに用いられる。
【0041】
さらに好ましい実施の形態では、反応速度定数はシミュレーションの前に適切な方法によって推定される。
【0042】
好ましい実施の形態では、上記反応速度定数の少なくとも10%、少なくとも20%、少なくとも30%、少なくとも40%、又は少なくとも50%、及びそれとは独立して、上記開始濃度の少なくとも10%、少なくとも20%、少なくとも30%、少なくとも40%、又は少なくとも50%が実験データである。明らかに、実験データである、少なくとも60%、少なくとも70%、少なくとも80%、又は少なくとも90%の反応速度定数を用いることも予想される。
【0043】
本発明の方法のさらに好ましい実施の形態では、前記反応速度則のさらに一部が知られており、該反応速度則の残りの部分はランダムに選択される。
【0044】
この実施の形態は、未知の反応速度定数及び/又は未知の開始濃度に加えて、モデルの生物学的実体間の相互作用を支配する反応速度則が部分的に未知である状況にまで及ぶ。反応速度則のランダムな選択は、離散確率分布から実行することができる。確率分布は、反応速度則が離散変数であることの結果として離散している。生物学的実体Aの生物学的実体Bへの転換の反応速度則が未知である事例においてさらに説明すると、反応速度則は、一次反応速度則に関して50%の確率をもたらし、二次反応速度則に関して50%の確率をもたらす確率分布から選択することができる。当然ながら、上述したように、大まかであるか又は類似の相互作用から導出される知識から利点を得ることができる。換言すれば、類似の事例において相互作用の90%が一次反応速度則によって支配され、10%が二次反応速度則によって支配されることが知られている(it is know)場合、分布は一次反応速度則に関して90%の確率をもたらし、二次反応速度則に関して10%の確率をもたらすことができる。
【0045】
より好ましい実施の形態では、上記反応速度則の少なくとも10%、少なくとも20%、少なくとも30%、少なくとも40%、又は少なくとも50%が実験データから導出される。明らかに、実験データから導出される、少なくとも60%、少なくとも70%、少なくとも80%、又は少なくとも90%の反応速度則を用いることも予想される。
【0046】
本発明による、生物学的実体の濃度を予測する方法のさらに好ましい実施の形態では、 (e)前記方法は1つ又は複数のさらなる生物ネットワーク上で付随して実行され、
(f)生物学的実体の濃度は、選択された時点において前記生物ネットワーク間で交換される。1つの生物学的実体の量若しくは濃度、又はより多くの若しくは全ての生物学的実体の量若しくは濃度は交換することができる。その量又は濃度が交換される好ましい生物学的実体には、成長因子、サイトカイン、及びホルモンのような細胞間シグナル伝達分子が含まれる。「量の交換」及び「濃度の交換」は、上記さらなる生物ネットワークのうちの1つ、複数、又は全てに対し、上記量又は濃度を利用可能にすることを指す。上記量がさらなる生物ネットワークに利用可能にされると、該さらなるネットワークにおける相互作用を支配する反応速度則に依拠して、それらの量を、該さらなる生物ネットワークにおける入力として用いることができる。
【0047】
この手法によって、中でも、生物システムの予期される挙動に対応したネットワーク、反応速度則、並びに反応速度定数及び開始濃度の値の範囲の特定の変形を選択することが可能になる。
【0048】
本発明のこの実施の形態によって、ネットワーク間の相互作用のシミュレーションが可能になる。たとえば、1つのシミュレートされたネットワークは細胞を表すことができ、第2のシミュレートされたネットワークは第2の細胞を表し、2つの細胞は情報及び/又は生物学的実体を交換することが可能である。明らかに、2つのみでなく、5つ、10個、数百個、又は数千個の生物ネットワークを付随してシミュレートすることができる。そのようなシミュレーションは、多細胞集合体、組織、器官、生物全体又は相互作用する生物群のシミュレーションである。
【0049】
生物学的実体は、シミュレーションの時間ステップにおいて毎回、又はより大きな間隔、たとえば1つおき、10個おき、又は100個おきの時間ステップにおいて交換することができる。
【0050】
本発明による方法のさらに好ましい実施の形態では、前記生物学的実体の濃度は、野生型と比較して摂動されている。
【0051】
本明細書において用いられる際、用語「摂動される」又は「摂動」は、野生型からの偏差を指す。対応する実験的設定において、そのような偏差は、所与の生物学的実体の過小発現若しくは過剰発現、又は突然変異に起因して生じる場合がある。他の予想される摂動は、薬剤又は他の物質の投与によって生じる。
【0052】
好ましい実施の形態では、摂動系における生物学的実体の濃度は実験的に求められる。モデル化される摂動系は、細胞、組織、器官、生物、又は相互作用する生物のグループとすることができ、摂動は、ノックダウン実験、突然変異、病態、又は薬剤の投与によって引き起こすことができる。
【0053】
本発明によるノックダウン摂動の好ましい実施の形態では、ノックダウンされている1つ又は複数の生物学的実体の濃度(複数可)は、それらの開始濃度の或る特定のパーセンテージに固定される。10%及び0%が好ましいパーセンテージである。加えて、ノックダウン実体の濃度を増加又は減少させる反応が無効にされるため、生物学的実体はシミュレーション全体を通じて固定の濃度に留まる。摂動される実体の開始濃度は、対数正規分布から選択されるか、又は遺伝子発現、RT−PCR、定量的プロテオミクス技術及びメタボロミクス技術のような実験データから選択される。
【0054】
本発明による突然変異摂動の好ましい実施の形態では、生物学的実体に対する突然変異の効果は、文献から既知であるようにモデル化されるか、又は突然変異の効果が未知である場合には、突然変異は、生物情報学技術からの推論を用いてモデル化される。たとえば、サイレント突然変異は、野生型の生物学的実体によって効果的にモデル化され、ミスセンス突然変異は多くの場合に生物学的実体の完全なノックダウン(0%)によってモデル化することができ、既知の機能ドメインを損傷する突然変異は、モデル化された生物学的実体と、損傷されたドメインが相互作用することを意図していた生物学的実体との間の適切なエッジを除去することによってモデル化することができ、構成的に活性化する突然変異は、生物学的実体の不活性形態を活性形態に転換する人工的な非可逆反応(エッジ)を付加することによってモデル化することができ、最後に、酵素生物学的実体の酵素効率を変化させることが知られている突然変異は、反応速度定数に効率の変化の既知の因子を乗算することによってモデル化される。これらの全ての事例において、反応速度定数は実験的に求められるか又は対数正規分布から選択される。
【0055】
本発明による、病態の摂動の好ましい実施の形態では、疾患に関与するとして知られる突然変異の摂動は、上記のセクションにおいて説明したようにモデル化される。加えて、遺伝子発現、タンパク質及びリンタンパク質濃度、代謝産物及びマイクロRNAレベルにおいて体現される活性病態データは、適切な生物学的実体の初期濃度を説明されたレベルに経験的に設定することによってモデルに直接適用される。
【0056】
本発明による薬剤投与摂動の好ましい実施の形態では、薬剤が相互作用することが知られている生物学的実体に対する該薬剤の適用の効果は、いくつかの異なる形のうちの1つでモデル化することができる。
【0057】
a)薬剤が、酵素、たとえばキナーゼである1つ又は複数の生物学的実体の活性を阻害することによって作用する場合、その活性は、キナーゼごとに実験的に規定されたIC50を取り、それを、薬剤の既知のIC50濃度を薬剤のモデル化された細胞濃度で除算してその結果に50%を乗算することによりキナーゼの反応速度定数に適用することによってモデル化される。モデル化された細胞濃度は、一般に、適用濃度とみなされる。たとえば、500nMの薬剤適用をモデル化するために、細胞濃度は一般に500nMであると想定される。しかしながら、モデル化される薬剤に関する因子、たとえば透過性若しくは溶解性、又は薬剤の薬理作用に影響を及ぼす、モデル化される細胞に関する因子、たとえば上方調節されたPGP若しくはMDR−1が知られている場合、モデル化される細胞濃度は、適用される濃度の数分の1に設定することができる。これが行われる場合、実験によるデータに基づくことが好ましい。
【0058】
b)薬剤が、タンパク質間相互作用を阻害することによって作用する場合、その相互作用、すなわちエッジの反応速度定数は、a)におけるような薬剤の既知のIC50によって変更される。
【0059】
c)薬剤が、ほとんどの従来の抗腫瘍薬のように非機構的である場合、薬剤の適用の効果は、薬剤に対する感度及び耐性を担うことが知られる、ネットワーク内の生物学的実体を作動させることによってモデル化される。たとえば、シスプラチン及びカルボプラチンのような、白金を用いた化学療法薬は、DNAをキレートすることによって作用し、この結果、細胞DNA修復ネットワークが過度に刺激されることになる。加えて(Iin addition)、細胞は、PGP及びMDR−1を過剰発現することによって、又はDNA損傷感知及び修復並びにアポトーシスネットワークにおける突然変異を獲得することによって、これらの薬剤に対して耐性を有するようになる。したがって、DNA損傷感知及び修復経路における選択された生物学的実体は、構成的に活性化されたものとして系内でモデル化することができる。加えて、非機構的薬剤の効果は、RNA及びタンパク質発現パターンにおける変化に基づいて間接的にモデル化することができる。このようにしてモデル化される実体(entitites)は、実際の薬剤の実際の細胞への適用の遺伝子発現又はプロテオメトリック実験から取られることが好ましいが、文献から薬剤応答に関して知られていることからモデル化することもできる。
【0060】
本発明の方法のさらに好ましい実施の形態では、初期条件は、
(a)生物学的実体の実験的に求められた濃度、及び/又は
(b)実験的に求められた突然変異データ
を含む。
【0061】
本明細書において用いられる際、用語「初期条件」は、シミュレーションの開始時の、上記生物学的実体の性質、数、状態、及び濃度を指す。
【0062】
本発明は、本発明による生物学的実体の濃度を予測する方法の統計学的有意性を求める、コンピューター実装される方法であって、該統計学的有意性を求める方法は
(a)本発明による生物学的実体の濃度を予測する方法を実行すること、
(b)ステップ(a)において取得された生物学的実体の濃度と、同じ生物学的実体について実験的に求められた濃度との間の一致度を求めること、
(c)前記生物ネットワークのトポロジーをランダム化すること、
(d)ステップ(b)において取得した前記ランダム化された生物ネットワークに対し、本発明による生物学的実体の濃度を予測する方法を実行すること、
(e)ステップ(d)において取得した生物学的実体の濃度と、同じ生物学的実体について実験的に求められた濃度との一致度を求めること、
(f)ステップ(b)において取得した結果をステップ(e)において取得した結果と比較することであって、ステップ(b)における一致度の方が高いことは、本発明による生物学的実体の濃度を予測する方法が、実験的に求められた濃度を偶然よりも良好に予測可能であることを示す、比較すること、
を含む、方法をさらに提供する。
【0063】
本発明のこの態様は、本発明による、生物学的実体の濃度を時間の関数として予測する方法を検証する方法に関する。該方法は、確率分布から反応速度定数をランダムに選択する効果を生物ネットワークのトポロジーのランダム化と比較する。好ましくは、生物ネットワークを前記ランダム化することは、前記ネットワークのエッジをスワップすることによって有効にされる。
【0064】
用語「上記ネットワークのエッジをスワップすること」は、該ネットワークにおける接続の変更を指す。たとえば、生物学的実体の濃度を時間の関数として予測するのに用いられる生物ネットワークにおいて、エッジ1がノードA及びノードBを接続し、エッジ2がノードC及びノードDを接続する場合、スワップされたエッジを有するネットワークの例は、エッジ1がノードA及びノードDを接続し(connections)、エッジ2がノードB及びノードCを接続するネットワークである。
【0065】
上記で言及した「一致度を求めること」は、利用可能な場合には経時変化を用いて、又は定常状態濃度若しくは平衡濃度のような最終濃度を用いて実行することができる。
【0066】
好ましい実施の形態では、エッジの10%、20%、30%、40%、50%、60%、70%、80%、90%、又は全てのエッジがスワップされる。
【0067】
本発明による方法のさらに好ましい実施の形態では、前記実体は生体分子であり、好ましくは遺伝子を含む核酸、タンパク質を含む(ポリ)ペプチド、小分子、並びに生体分子の複合体及び代謝産物から選択される。小分子には、サッカリド、アミノ酸、脂質、ヌクレオチド、ヌクレオシド、並びにそれらの代謝産物及び誘導体が含まれる。
【0068】
本発明による方法のさらに好ましい実施の形態では、上記モデルは、境界条件、好ましくは生理学的状態を表す境界条件を含む。境界条件は、反応速度定数に対して、かつ/又はグラフの接続性に対して影響を有する場合がある。
【0069】
好ましい境界条件は、生物学的実体又は刺激の存在又は恒常性を含む。好ましい境界条件は、所与の量の1つ又は複数の薬剤の存在を含むことができる。他の好ましい境界条件は、細胞外シグナル伝達勾配、並びに細胞間の伝達及び生理学的シグナルによって課された境界条件である。
【0070】
本発明は、先行する請求項のいずれか1項に記載の方法を実行するように構成されたコンピュータープログラムをさらに提供する。
【0071】
本発明によるプログラムを含むコンピューター可読データキャリアがさらに提供される。本発明による方法を実行する手段を含むか、又は本発明によるプログラムがインストールされたデータ処理装置も提供される。
【0072】
本発明は、生物ネットワークの部分的に未知のパラメーターを求める、コンピューター実装される方法であって、該パラメーターは、ネットワークトポロジー、反応速度則、反応速度定数、及び/又は開始濃度から選択され、該方法は、観測される特性と予測される特性との差異を最小にすることであって、該予測される特性は、本発明による生物学的実体の濃度を予測する方法によって予測されるような濃度を含む、最小にすることを含む、方法をさらに提供する。たとえば、或る特定の生物学的実体の定常状態濃度又は平衡濃度は、容易に実験的に求めることができる。本発明による生物学的実体の濃度を予測するin silicoの方法に関して、これらの生物学的実体の予測される定常状態濃度又は平衡濃度は、多様な程度で、未知のパラメーターに割り当てられる値に依拠する。観測される特性と予測される特性との間の差異を最小にすることによって、未知のパラメーターは、取得される反応速度モデルが、実験とシミュレーションとの間の差異が最小にされた特性を最も良好に再現するモデルであるという意味で最適化される。この最小化は、連続最適化、すなわち、本発明による生物学的実体の濃度を予測する方法の実行中の、上記差異の最小化を伴う場合がある。用語「最適化」は、非線形回帰型手法の使用等による、上記で規定したパラメーター、すなわち、たとえば反応速度定数及び/又は開始濃度の最適化に関する。さらに、最適化は、代替的に又は付加的に、ネットワークのトポロジー及び/又は反応速度則の変更のような不連続ステップを伴う場合がある。好ましくは、予測値と実験値との間の一致は、全ての利用可能な測定値に関して最適化される。観測値と予測値との間の差異、又は観測されるパラメーターと予測されるパラメーターとの間の差異を最小にする汎用計算手段及び方法は当該技術分野において既知である。好ましい実施の形態では、上記ネットワークトポロジー及び/又は上記反応速度則は完全に既知である。
【0073】
本発明は、1つ又は複数の実験を選択する、コンピューター実装される方法であって、
(a)本発明による生物学的実体の濃度を予測する方法を実行することによって複数の実験をin silicoで実行することであって、該実行することは未知のパラメーターの異なる選択を用いて実験ごとに反復して行われ、該パラメーターは、ネットワークトポロジー、反応速度則、反応速度定数、及び/又は開始濃度から選択される、実行すること、並びに、
(b)前記複数の実験から、前記異なる選択に依拠して、前記生物学的実体の濃度を予測する方法が予測濃度の最大の分散をもたらす1つ又は複数の実験を選択すること、
を含む、方法をさらに提供する。
【0074】
用語「実験」は、特記されない限り、実行される実世界の実験に関する。この実験は、生物システムにおいて又は生物システムに対してそれぞれ、1つ若しくは複数の遺伝子のノックダウン、1つ若しくは複数の遺伝子に特異的な1つ若しくは複数のsiRNA(低分子干渉RNA)の投与、及び/又は酵素のような1つ若しくは複数のポリペプチドの活性の1つ若しくは複数の修飾因子の投与を実行することを含む実験とすることができる。好ましくは、上記生物システムはin vitro系である。好ましい修飾因子は、薬剤又は薬剤の開発に適したリード化合物である。
【0075】
換言すれば、これらの実験を選択する、コンピューター実装される方法であって、モデルの異なる代替形式間を最も効率的に区別し、全ての可能なモデルを用いて実験をモデル化することによって未知のパラメーターを求めるのに最も効果的となる方法が提供される。実験的に求めることができる、濃度のような特性の予測値間の差異が最大である実験は、異なる形式のモデルのうちのいずれが正確である可能性がより高いかに関するほとんどの情報を提供し、したがって、可能であるがおそらく与える情報がより少ない実験のより大きなセットの中から選択されることができる。
【0076】
1つ又は複数の実験を選択する方法の好ましい実施の形態において、上記ネットワークトポロジー及び/又は上記反応速度則は完全に既知である。好ましくは、上記選択は、上記で詳述したようにランダムな選択である。さらに好ましい実施の形態では、選択された実験(複数可)が実行される。
【図面の簡単な説明】
【0077】
【図1】デビッドソン(Davidson)他[14、4]から再現される内胚葉系中胚葉ネットワークトポロジーを示す図である。トポロジーは、バイオタペストリ(Biotapestry)[15]を用いて生成される。
【図2】alx1−mRNA及びotx−mRNAのシミュレートされた経時変化を示す図である。mRNA量は任意の単位で与えられ、時間は受精後時間で与えられる。800回のシミュレーション実行の平均が各時点においてプロットされる。
【図3】pmar1発現の摂動によって影響を受ける遺伝子の発現を表す図である。mRNA量は任意の単位で与えられる。検出される摂動の効果は、[25]に記載されている実験データに従う。
【図4】無摂動のpmar1−MASO条件(上)及び無摂動のalx1−MOE条件間の、hesc−mRNA(左)及びalx1−mRNA(右)の量を比較する散布図である。x軸は摂動条件のmRNA量を表し、y軸は無摂動条件下のmRNA量を表す。量は、受精後25時間(hpf)の時点でa.u.単位で測定される。800個の異なるパラメーターセットの値がプロットされる。対角線からの偏差は、摂動条件の発現がより高い(対角線の下)か、又は摂動状態の発現がより低い(対角線の上)ことを示す。
【図5】アポトーシスの誘導の前(左)及び後(右)の切断カスパーゼ9の量を示す図である。緑は80%カスパーゼ3ノックダウンに対応し、青は対照に対応し、赤は2つの間の差異に対応する。バー(y軸)の長さは、複数のモンテカルロ−シミュレーションをそれぞれの結果(x軸)と共に示す。初期タンパク質濃度はランダムに選択され、400個全てのシミュレーションに適用される(apply)。
【図6】上枠(A)は、アポトーシスのシミュレートされた誘導の前(左)及び後(右)のカスパーゼ9(未切断)の量を示す図である。下枠(B)左は、スタウロスポリン(staurosporin)によるアポトーシスの誘導を用いない場合(0h+)及び用いた場合(6h+)の、カスパーゼ3siRNAを96時間形質移入されたWI38細胞におけるカスパーゼ9の発現を、アポトーシスの誘導を用いない場合(0h−)及び用いた場合(6h−)の形質移入されていない細胞と比較したウエスタンブロットである。下枠(B)右は、AIDAを用いたウエスタンブロットの定量化を示す図である。値はアクチンを用いて標準化される。
【図7】注釈及びネットワークモデル過程を示す図である。(A)は、癌治療の最小ネットワークを含むシグナル伝達経路の一部の図である。(B)は、PyBioS系内のコンピューター可読反応系におけるこれらの機能的相互作用の説明の図である。
【図8】提案されるモンテカルロ手法のフローチャートである。定常状態シミュレーションが正常状態及び摂動状態(薬剤による阻害)に関して実行される。阻害される成分を除く全てのモデル成分に関して、正常状態及び摂動状態が、反応速度パラメーターの同じセット及び初期値を用いて初期化される。その後、モデルは該モデルの定常状態に入るようにシミュレートされる。この手順は、所与のランダム分布からサンプリングされる異なるパラメーターベクトルを用いてk回反復される。グラフィック検査に関して、それぞれの対照及び処理シミュレーション実行の定常状態の結果は、モデルの全ての成分についてヒストグラムを用いてプロットされる。さらに、対照及び処理間の結果の分布は、対照/処理値のセットごとにコルモゴロフ−スミルノフ検定によって計算されるP値によって定量化される。
【図9】全体統計を示す図である。(A)は、阻害実験ごとに有意に変化したモデル成分の数によって、処理の個々の感度を示す図である。(B)は、ヒストグラムによって、阻害実験の集団(panel)に対するモデル成分の感度を示す図である。このヒストグラムは、この研究において分析される24個の阻害実験のうちの所与の数(X軸)によって有意に変更されるモデル成分の数(Y軸)を示す。(C)24個の異なる阻害実験の対数比の主成分分析(PCA)を示す図である。PCAはJ−Express Pro(ベルゲン、Molmine社)を用いて行われた。(767個のうちの)2つの第1の主成分によって全分散の67%が明らかになる。
【図10】AKTシグナル伝達を標的にする特定の阻害実験を示す図である。AKTによって調節される複合ネットワークを標的にする4つの薬剤/小分子阻害剤(ペリホシン、ウォルトマンニン、メトホルミン、及びラパマイシン)を示す図である。変化は増殖、成長、及びアポトーシスに影響を与える。阻害は、点線(blunted line)によって示される。IRSはインスリン受容体基質であり、AKTはタンパク質キナーゼBであり、Pl3Kはホスファチジルイノシトール3−キナーゼであり、mTORはラパマイシンの哺乳類標的である。
【図11】薬剤標的阻害の直接的及び間接的効果を示す図である。選択された対照状態(上側バー)対処理状態(下側バー)の定常状態頻度のヒストグラムを示す。選択結果は、図15に示される特異的阻害実験の直接的な効果(左枠)及び間接的な効果(右枠)を示す。(A)は、阻害されたAKT2の、「ホスホ−GSK3β」の濃度に対する直接的な効果を示す図である。(B)は、AKT2阻害の間接的な効果が、「S6Kリン酸化」の阻害を示すことの図である。(C)は、PDK1が、「PIP3:リン酸化PKB」の阻害における直接的な効果を示すことの図である。加えて、(D)は、「ホスホ−GSK3β」の阻害が、治療ドメインにおいて間接的な効果とみなされることを示す図である。(E)は、IRSが「ホスホ−IRS:Pl3K」の形態でPl3Kの阻害を通じて直接的な効果を与えることを示す図である。(F)は、活性化IRSの阻害の結果、「TSC2−P」リン酸化の阻害が生じることを示す図である。(G)は、治療ドメインAKTPI3Kが、「ホスホ−GSK3β」に対する直接的な効果を示すことの図である。(H)は、治療ドメインAKTPI3Kが、「S6K1−P」に対する間接的な効果を示すことの図である。合計250回のシミュレーション実行が、阻害実験ごとに行われた。実際には、或る特定の量のシミュレーション実行を超えた実行時間制約に起因して複数のシミュレーション実行が失敗することにより、結果としてのシミュレーション実行回数はこれより小さい(平均220回)。
【図12】AKT阻害を伴う薬剤作用の共発現クラスターを示す図である。異なる阻害実験(列)及びモデル成分(行)に対する、処理状態及び対照状態の平均定常状態値の対数比(底2)を示す図である。赤色は阻害に対するモデル成分の上方調節を示し、緑色は下方調節を示す。対類似度測定としてピアソン相関が用いられ、更新規則として平均連結法が用いられた。J−Express Pro(ベルゲン、Molmine社)を用いてクラスター化が実行された。
【図13】AKT(RACβセリン)に対して高い感度を示すモデル成分のクラスターを示す図である。この研究に用いられた29個の異なる単一阻害実験にわたって同様の感度パターンを示すモデル成分のクラスターを示す図である。薬剤標的ごとの感度値(列)が全てのモデル成分(行)について計算された。−20〜20の段階の感度値がプロットされる。色の段階で示されるように、高い感度成分は赤又は緑である。J−Express Pro(ベルゲン、Molmine社)を用いてクラスター化が実行された。以下の実施例は、本発明を説明するが、限定的なものとして解釈されるべきでない。
【図14】摂動シミュレーションの定性的結果を示す図である。緑色の細胞は、列の摂動に起因して行内の遺伝子の発現が増加したことを示し、赤は発現が低減したことを示す。細胞は、或る特定の領域に含まれた効果を強調するために網掛けにされる。影響を受ける領域は個々の細胞において言及される(E:内胚葉、M:中胚葉、P:PMC、T:全ての領域の合計)。一時的な発現情報は表から除外された。それによって、遺伝子が何らかの所与の時点において有効にされた場合、その効果は遺伝子一般に割り当てられる。
【図15】シミュレーションと実験データとの比較の定性的概観を示す図である。表を短くするために、異なる時点の比較結果が結合されている。符号化は以下の通りである。緑は摂動及び有効にされた遺伝子の全ての実験データが整合していることを示し、緑及び黒のパターンは、所与の摂動及び遺伝子に関して不整合よりも整合の方が多く存在することを示し、灰色は、同程度の数の整合及び不整合が存在することを示し、赤及び黒のパターンは整合よりも不整合の方が多く存在することを示し、赤は不整合のみを示す。最終行は、各列の平均を示す。
【発明を実施するための形態】
【実施例】
【0078】
実施例1:内胚葉系中胚葉ネットワーク
背景
受精卵からの成体の発生は、基本過程と同様に複雑である。発生は、周囲の細胞からのシグナルによって駆動される個々の細胞の特異化、並びに細胞運動性及び卵割事象を含む。主な調節的入力は多数の細胞によって生成されるが、これらの巨視的効果を生成する微視的事象は、細胞レベルで精密に調節されなくてはならない。このため、発生メカニズムは、シグナル伝達事象及びタンパク質相互作用、並びに遺伝子調節及び細胞間相互作用を含む。これらの発生メカニズム、及び異なる種間のこれらのメカニズムの差異を理解することによって、進化メカニズムを見通すことができる[1]。さらに、発生中のいかなる摂動も、何らかの形で生物において顕在化する可能性が高い。
【0079】
発生の科学的分析は、ウニを用いて1890年代に始まった[2]。ウニは歴史的理由からのみでなく、その興味深い進化位置からもモデル生物である。進化プログラムの基本過程は、哺乳類の発生にも同様のものがあると予期される[3]。
【0080】
今日、ウニの発生の理解における最も野心的な取り組みは、アメリカムラサキウニ(S.pur.)における内胚葉系中胚葉ネットワークである[4]。このネットワークは、初期のウニ胚における発生を制御及び調節する複雑な遺伝子調節ネットワーク(GRN)、主に内胚葉、中胚葉、及び第一次間充織(PMC)領域の特異化の構造を捉えようと試みている。遺伝子間の相互作用は、摂動実験、主に、モルホリノ基置換アンチセンスオリゴヌクレオチド(MASO)注入が、特定の遺伝子のmRNA翻訳[5]、mRNA過剰発現、及びエングレイルドリプレッサードメインとの融合を効果的にノックダウンすることから推測される。このネットワーク(図1)は静的であり、相互作用のタイプ、タイミング、及び強さに関する規則を規定する必要がある。本発明の方法を用いて、内胚葉系中胚葉ネットワークが常微分方程式(ODE)を用いてモデル化された。実験データは主に摂動実験[4]に基づき、詳細な研究は限られた数の遺伝子に関してのみ実行された[6、7、8、9、10、11、12、13]ため、実験データは結果としてのモデルを完全にパラメーター化するには不十分である。
【0081】
データがわずかであるにも関わらず時間依存的挙動を予測するために、ランダムにサンプリングされたパラメーターを用いてシミュレーションを実行した。パラメーター変動に対する各ネットワーク成分の感度を検出することに加えて、ネットワークが基づく実験データに関する摂動から、パラメーター変動に対してロバストな効果が検出される。これは同じパラメーターセットを用いて、正常条件及び摂動条件の下でモデルをシミュレートすることによって達成される。結果としてのシミュレーション対の比較によって、用いられたパラメーター値とは独立した効果が示される。
【0082】
比較結果は、モデルトポロジーがどの程度実験データを再現することができるかを示す。ネットワーク内の全ての相互作用がパラメーター変化に対して不変であることを一見して予期することはできない。
【0083】
モデルのランダム化されたバージョンを用いて検証が行われる。元のモデルを用いて実行したのと同じ分析が該モデルのランダム化されたバージョンに適用される。元のモデルが実験データを再現する程度を、ランダム化されたバージョンが該データを再現する程度と比較することによって、元のモデルの妥当性が推論される。
【0084】
方法
ネットワーク妥当性の推論
ネットワークトポロジーの妥当性を推論するために、ネットワーク構造を含む入力ファイルがPyBioS[29、30]においてOEDモデルに転換された。結果としてのモデルは、Python[31]で構築され、容易にシミュレートすることができる。
【0085】
問題となっている異なる胚領域のより現実的なシミュレーションを作成するために、全て互いから独立している、同じ元モデルの3つのコピーを含むモデルが作成された。異なる外部入力を系に適用するため、これらの入力によって異なる領域を判別することができる。PyBioS出力を用いて、モデルのシミュレーションのために、転写調節のためのランダムパラメーターのセットがサンプリングされる。モデルは、70個の時間ステップにわたってシミュレートされた。ここで、1つのタイムステップは受精後1時間(hpf)に対応する。外部で設定された入力を用いて胚領域を判別するため、これらの入力はタイマーとしての役割を果たし、実験データに従って発現の時間フレームを確立する。
【0086】
ノックダウン実験をモデル化するために、問題となっているmRNA生成の速度則を設定することによってPyBioSシンタックスにおけるモデルが、vtranscriptionからvtranscription×0.05に変更された。過剰発現実験の場合、元のモデルにおいて到達した最大発現率が任意の単位で約2/hであるため、vtranscriptionはvtranscription+2に変更された。ランダムにサンプリングされたパラメーター値を用いると、この値は任意となり、MOE実験の分析における誤りの原因となり得る。
【0087】
上述した方法で変更された任意の数のモデル及び元のモデルは、同じパラメーターセットを用いてシミュレートすることができる。この研究において用いられるパラメーターセットは、σ=1.5、μ=0.5の対数正規分布からサンプリングされた。これは、本発明による好ましい分布である。この分布は、一つには、ODE系の数値的安定性を妨げ得るあまりに極端なパラメーター値を回避するために選択された。この研究では、シミュレーションのために800個の異なるパラメーターセットを用いた。
【0088】
全ての可能なパラメーターセットSの下で、有効にされた可能性がある遺伝子Xのそれぞれにつき1つの特異摂動pに関するシミュレーション結果のセットごとに、次のステップが計算される。
【0089】
T=(14,19,25,33,45,66)となるような、[19]において与えられる時間間隔を反映する時点のリストTが規定される。
【0090】
条件p及び無摂動条件C下の全てのパラメーターセットSに関する結果が対で整列される。
【0091】
t∈Tごとに、最も近いシミュレーション時点が求められ、双方の条件に関するシミュレーション結果、及びパラメーターセットs∈Sごとの比rx,s,t=[x]C,s,t/[x]p,s,tが記憶される。これらの値を用いて、摂動及び遺伝子ごとに、全てのパラメーターセットSにわたる平均rx,S,t及び分散varx,S,tを計算する。これが3つ全ての領域及び3つの領域の和について計算された。
【0092】
1つの遺伝子xに関する値が、視覚定位に関して図15に散布図で示されている。摂動の効果は、或る特定の閾値を上回っているか又は下回っている発現比の数として測定される。具体的には、或る特定の閾値zを上回っている比又は或る特定の閾値zを下回っている比を数える。
【0093】
1つの所与の時点tについて、T個より多くのSに関してrx,s,t>z(ただしz≧1)である場合、遺伝子xは時点tにおいて摂動pによって有意に下方調節されるといわれている。他方で、T個より多くのS及び1つの特定のtに関してrx,s,t<z(ただしz<1)である場合、xは時点tにおいて摂動pの下で有意に上方調節されるといわれている。この分析において、z=z=1、及び閾値T=90%を設定する。有意でない変化から判別するために、SのうちのT個について0.8≦rx,s,t≧1.25である場合の効果は有意でないものとみなした。計算された全ての摂動に関する全てのシミュレーション結果を用いて、それらを利用可能な実験データと比較すると、すなわち、実験結果にマッチするシミュレーション結果における上方調節/下方調節/無差別の割合として、モデルがどの程度実験データを再現するか、すなわちモデルの妥当性を計算することが可能である。異なる領域の独立したシミュレーションから全体的な妥当性を推論するために、領域のうちの少なくとも1つにおける実験データとの整合を、モデルと実験データとの間の有意な整合として数えることを選択した。
【0094】
シミュレーションのセットが実験データの所与のセットを再現する程度は、z、z、Tの選択及び有意でない効果に関する閾値に依拠する。系の知識が限られていることに起因して、比較的緩い閾値を設定した(set chose)。モデルのサイズがより小さいか、又は詳細で包括的なモデルが利用可能であるサブセットが用いられる場合、これらの閾値を厳しくすることを選択してもよい。
【0095】
ロバスト性の推論
内胚葉系中胚葉ネットワークモデルの異なる成分のロバスト性の推論は、無摂動モデルのシミュレーション結果のみを用いる。全ての利用可能なパラメーターセットからのシミュレーション結果を比較することによって、遺伝子ごとのシミュレーション結果がパラメーター値に依拠して異なる程度が抽出される。これらの変動の程度は、パラメーター値に対する遺伝子の発現のロバスト性である。
【0096】
遺伝子の発現の1つの真の経時変化が存在すると想定され、これはサンプリングされたパラメーターセットによって必ずしも発見されない。それにもかかわらず、異なるパラメーターセットに関するシミュレーション結果は、経時変化の個々の時点を分析するとき、真の経時変化の周りに正規分布することがさらに想定される。
【0097】
この想定によって、異なる時点において、遺伝子ごとにmRNA濃度の平均μ及び分散σを計算することが可能になる。平均に対する分散の程度を測定するために、絶対平均値又は絶対分散値とは独立して、時点と遺伝子との間でσ/μで比較することができる相対変動varrel=を計算することができる。これは基本的には、対応する平均に対する分散を与える。
【0098】
遺伝子ごとに十分な数の時点を選択するとき、varrelのリストを用いて、遺伝子の発現の全体的なロバスト性を得ることができる。ここでは、遺伝子ごとに最大varrelを選択することによって行われる。結果としての値が実質的に異なり、異なるレベルのロバスト性を示す場合がある。varrel及びそれらの平均は、ロバストな遺伝子と脆弱な遺伝子との区切りを求める絶対値としては解釈されない。なぜなら、これらの区切りは未知である実際のパラメーター値に大きく依拠するためである。
【0099】
ネットワークのランダム化
PyBioS形式のODEモデル。これは、閾値モデルの場合、2つの遺伝子をランダムに選択し、領域ごとに2つのランダムに選択された入力をこれらの遺伝子と交換することによって行われる。このランダム化手順を3000回反復する2つのランダムネットワーク、及びランダム化ステップを30000回反復する1つのモデルを作成した。3つの領域間の境界が維持されることに留意されたい。このランダム化アルゴリズムを用いて、ノード及びエッジの数、平均ノード次数及び次数分布のようなネットワークの全体特徴は、個々の配線が変更される一方で保持されることに留意されたい。ランダム化されたネットワークが構築された後、上述したようにその妥当性を判断することができる。この分析において、3000個のランダム化ステップを用いて2つのランダム化モデルを構築し、それぞれを100個の異なるパラメーターセットを用いてシミュレートした。より強力なランダム化の効果を推論するために、30000個のランダム化ステップを用いる1つのモデルも構築し、該モデルを20個のパラメーターセットを用いてシミュレートした。
【0100】
結果及び検討
内胚葉系中胚葉ネットワークの動的モデル
内胚葉系中胚葉ネットワーク[4]は、ウニS.purにおける内胚葉、中胚葉、及びPMCの分化を駆動する遺伝子間の予測される調節的相互作用を示している。このネットワークの精緻化されたバージョンが基本データ[19]と共に利用可能である[14]。このデータに加えて、いくつかの遺伝子の調節的相互作用が詳細に研究された[6、7、8、9、10、11、12、13]。
【0101】
摂動及び利用可能な経時変化データ(言及したソースにおけるもの)は一般に、胚全体について測定されるが、定性的データが利用可能であり、或る特定の遺伝子が或る特定の領域のみで発現されるか、又はさらには複雑な時空間発現パターンに従うことが示される[7、6、20]。この差次的発現は、胚の細胞間の直接的相互作用及び間接的相互作用によって駆動及び実行される[21、22]。
【0102】
各領域は、遺伝子の発現、TF及びシグナル伝達分子の存在量に関して異なるものの、各細胞が同じ遺伝子情報を含む、すなわち、発生の初期段階においてヒストン修飾が起こらないと想定する。さらに、各領域が、一様な数の細胞から構成される、すなわち、同じ領域内の細胞がTFの同じ組み合わせを含み、同じ遺伝子を発現すると想定する。
【0103】
これらの想定によって、1つの細胞のみをモデル化することにより各領域(内胚葉、中胚葉、及びPMC)をモデル化することが可能になる。このため、同じメカニズムの3つの複製を含むモデルを構築する。このため、モデル化された細胞間の差次的発現は、細胞間シグナル伝達における差異及び異なる開始条件からのみ生じることができる。内胚葉、中胚葉、及びPMCは胚全体を構成するものではなく、細胞間シグナル伝達の正確なメカニズム及び胚内の拡散率は未知であるため、これらのメカニズムは除外することにした。基本メカニズムを除外したにも関わらずこれらの異なるシグナル伝達事象の効果を含めるために、[19]に含まれるデータを用いて、モデル自体から出現しない各領域に対する一時入力を模倣する。これらの適切な初期条件、及び細胞外過程を模倣する人工的入力を所与とすると、モデルは、内胚葉細胞、中胚葉細胞、及びPMC細胞の挙動を再現することが可能であるはずである。
【0104】
この系に対し、単一の遺伝子のノックダウン(KD)又は過剰発現の形態で摂動を加えると、各単一の領域に対する効果を求めることができる。効果を結合して実験データにおいて測定されるような全体効果にするための、所与の時点における胚の正確な組成を知らないため、領域のうちの少なくとも1つにおける摂動の効果の検出によって、それに従う胚全体の挙動を示すことを提案する。そのために、このモデルを用いて、既存のトポロジーを検査することができる一方、新たな実験データと、代替的な仮説の検査とを正規に統合することも容易になる。
【0105】
そのような用途において、常微分方程式(ODE)のモデルが広く用いられている。それらのモデルは定常状態の分析及び経時変化の詳細なシミュレーションを可能にする[23]。それらのモデルは、より単純なモデル化の枠組みよりも詳細な結果を生成するため、ODEのモデルは、モデル化された系に関する、より詳細な情報も必要とする。
【0106】
本明細書において上記で説明した反応速度則を用いて、内胚葉系中胚葉ネットワークのトポロジーを用いるOEDモデルを構築した。図1を参照されたい。このモデルは、細胞タイプごとに合わせて54個の遺伝子及び140個の変異種から構成される。mRNA及びタンパク質の転写、翻訳、分解、並びに複合体の会合及び解離を含めて、モデルは、細胞タイプごとに、287個のパラメータによって制御される合計278個の反応を含む(翻訳、mRNA分解、及びタンパク質分解は1つの偏在性パラメーターによって制御される)。ここでもまた、これらの数はパラメーター推定がこのネットワークに関して現在実現不可能である理由を示している。
【0107】
図2は、異なる遺伝子及び領域に関してシミュレートされた経時変化を示している。これらの経時変化は実験による経時変化を再現しない場合があるが、これはこの分析の目標では決してなく、実際にパラメーター推定を必要とした。それにもかかわらず、これらの経時変化は、本発明者らのモデルが差次的発現を生成することが可能であることを実証する。alx1−mRNA量の経時変化は、PMC条件下でのみ発現されることを明確に示している一方、otxは3つ全ての領域において発現されるが、各領域において異なる程度で発現されることを明白に示している。
【0108】
方法の評価
本明細書において説明される方法は、12個の遺伝子から構成される内胚葉系中胚葉ネットワークのサブモデルに対して検査された[24]。該サブモデルはわずかに変更を加えられ、該サブモデルに関して、摂動データではなく既知の経時変化データを再現するためのパラメーターが推定された。推定されたパラメーターが真のパラメーターであり、かつネットワークの動的挙動が正確であると想定すると、これを、ランダムにサンプリングされたパラメーターを適用してネットワークのトポロジー特徴を抽出するためのベンチマークとして用いることができる。
【0109】
このサブネットワークに関して、推定されたパラメーターを用いて取得されたシミュレートされた摂動効果が、ランダムにサンプリングされたパラメーターを用いたシミュレーションから取得されたシミュレートされた摂動効果と比較された。推定されたパラメーターを用いてモデルをシミュレートすることによって達成される効果は、真の実験データの代わりに基準として用いられた。
【0110】
このサブネットワークの分析によって、ランダムにサンプリングされたパラメーターを有するモデルを用いた摂動効果の検出の受入可能な精度が示される(双方の設定において摂動効果の73.8%が等しく検出された)。偽陰性(推定されたパラメーターの使用において効果が検出され、サンプリングされたパラメーターセットの使用において効果が検出されない)が8.1%を構成し、偽陽性が15.7%を構成し、2.1%が2つの設定において効果が逆の符号を有して検出されたことを示す。
【0111】
これは、サンプリングされたパラメーターセットを用いてOEDモデルのトポロジー特徴を推論することによって、少なくとも小さなネットワークに関して満足のいく結果がもたらされることを示している。多数のシミュレーションに必要とされる計算時間は、選択されるネットワークトポロジー及び反応速度を含む複数の要因に依拠するため、推定されたパラメーターの使用とサンプリングされたパラメーターとの間の正確な比較は、この単一の比較に基づいて達成することができない。それにもかかわらず、パラメーター推定は、パラメーター相互作用及び利用可能なデータに依拠してサンプリングされたパラメーターを用いて効果的にシミュレートすることができる系に関して実現不可能であり得る。
【0112】
ランダムなパラメーター変動に対する遺伝子発現のロバスト性
方法のセクションにおいて詳述したように、遺伝子の発現のロバスト性は、相対変動(varrel=σ/μ)を計算することによって、対応する変動を用いて様々な時点において異なるパラメーターセットを用いて発現の平均を比較することによって測定される。無摂動条件下で800個のパラメーターセットを用いてモデル内の遺伝子ごとにvarrelを計算した。遺伝子ごとに、全ての時点の最も高いvarrel(最も低いロバスト性に関する)が遺伝子のロバスト性を示すものとして用いられる。一般に、ロバスト性は、より初期の測定点において、より高い。これは、問題となっている遺伝子に到達し該遺伝子を生み出すのに時間がかかる、分析される遺伝子の上流の遺伝子の発現における変動に起因する。
【0113】
表1は、ネットワークにおける各遺伝子のvarrelに対する概観を与えている。varrelは用いられる800個のパラメーターセットに関して0.21〜25.69の範囲にわたり、ネットワークにおける遺伝子が、実際は、ランダムなパラメーター変化に対するそのロバスト性において実質的に異なることを示している。
【0114】
【表1】

【0115】
表1:ロバスト性によってソートされた、異なる遺伝子のロバスト性の概観。各遺伝子のロバスト性に関連付けられるスコアは、方法セクションにおいて説明されたように相対変動(varrel)である。800個の異なるパラメーターセットからの結果が用いられる。最小スコア(最大のロバスト性を示す)に関してソートする以外はクラスター化又はグループ化はこの表に適用されない。第3の列は遺伝子の入次数(入来する相互作用の数)を与え、第4の列はノード出次数(遺伝子が有する調節的相互作用の数)を示す。
【0116】
ノード入次数とロバスト性との間の相関(発現のロバスト性が調節因子が増えると共に増加することを示す)を検出することもできないし、ノード出次数とロバスト性との間の相関(重要な調節的遺伝子が他の遺伝子よりもロバストであることを示す)を検出することもできなかった。ネットワークにおける遺伝子の位置(初期内中胚葉、初期PMC、後期内胚葉及び中胚葉、後期PMC)とその(it's)ロバスト性との間の相関性を見出すこともできなかった。
【0117】
これは、現在のアーキテクチャがランダム摂動に対するロバスト性においていかなる特定の遺伝子のセットを優遇することがないか、又はネットワークが、同様のロバスト性を用いて表現される遺伝子のセットを検出するには小さすぎることを示す。
【0118】
異なる摂動のシミュレーション
内胚葉系中胚葉ネットワークの妥当性を推論することを可能にするために、異なる摂動実験をシミュレートする必要があった。ノックダウン実験及び過剰発現(OE)実験のモデルを効率的に作成するために、転写の速度則、たとえばvtranscriptionを、ノックダウンの場合にはvtranscription・0.05に、過剰発現の場合はvtranscription+2に変更した。
【0119】
エングレイルドリプレッサードメインを用いた遺伝子の融合の結果発生する効果をシミュレートしなかった。なぜなら、これは全ての影響を及ぼされる速度則を注意深く調整することが必要とするためである。
【0120】
1000個のパラメーターセットを用いて、元のモデル及び摂動モデルの双方をシミュレートした。非常に硬いODEを生成したサンプリングされたパラメーターセットから生じる数値の問題に起因して、801個のシミュレーション実行のみが終了した。
【0121】
たとえば、alx1、ets1、tbr、及びhescのようないくつかの遺伝子に対するpmar1摂動の効果を提示する。異なる摂動条件下のmRNA量の平均を経時変化図3としてプロットする。効果を視覚化する別の方法が図4に与えられている。ここでは、或る特定の時点における1つのmRNAの無摂動及び摂動条件下の量が、すべてのパラメーターセットに関してプロットされる。視覚的効果は、hesc発現に対するPmar1の阻害的役割を示すと共に、alx1、ets1、及びtbrに対するHesCの阻害的役割を示す、[25]に説明される所見に従う。本発明者らのデータは、この阻害が取り除かれる場合に、alx1、ets1、及びtbrがモデルに従ってより高い率で発現されることをさらに示している。
【0122】
全てのシミュレートされた摂動実験を視覚化したものが図14に与えられている。結果は、摂動された遺伝子の標的を示す摂動によって分析することができるか、遺伝子の調節因子を示す有効にされた遺伝子によって分析することができる。列方向の分析によって、実際に異なる胚領域に関連付けられる同時調節された遺伝子のグループが明確に同定される。たとえば、alx1、dri、ets1、hnf6、及びpmar1のノックダウンはPMC条件下のみにおいて検出可能な効果を有する。調節的効果の一致パターンを有する行は、同様の領域/調節因子に関連付けられる遺伝子のグループを示す。図1の左下部分からの針状のマトリックス遺伝子のグループは、異なる遺伝子の摂動に対するそれらの感度を考慮すると、右下部分の二次間充織細胞遺伝子のグループと大幅に異なる。GataE及びHoxは概して反対の効果を有し、これはHoxがgatae発現に対して有する阻害的役割によって説明することができる。hox−KDの甚大な効果を、ネットワークトポロジーにおけるその比較的限定された役割と比較すると、多数の検出された効果がgatae発現の阻害に起因することが明らかである。ネットワークトポロジーから明らかであるが、これはシミュレーション結果のみから判別することがむしろ困難である。シミュレーション結果のみを考慮すると、Hox及びGataEは共に全ての有効にされた遺伝子を並行に調節していると誤解される場合がある。
【0123】
内胚葉系中胚葉ネットワークトポロジーの妥当性
内胚葉系中胚葉ネットワーク(図1)の妥当性は、シミュレーション結果において正確に再現される実験結果の比に関して評価される。計算分析の結果は定性的結果しか生じないため、実験データは定量的データから定性的データに変換されることを必要とした。既存の実験データのうちのいくつかは、データが不明瞭であるために除外された。たとえば、gatae−MASOに応答したfoxbの発現は、−2.8/−2.4/−3.4、−2.4/NS、+3.7である(ただし、スラッシュは異なる実験からの測定値を示し、カンマは反復測定を示し、NSは有意でないと考えられる値を示す)[19]。このため、妥当性は不明瞭でない実験結果のみを用いて計算される。元の内胚葉系中胚葉ネットワークのモデル、並びに[19]に示されるMASO及び過剰発現実験に関連するモデルを用いて、元のモデル及び摂動モデルにおける遺伝子の発現が比較される。本発明者らのモデルは、内胚葉特異的濃度、中胚葉特異的濃度、及びPMC特異的濃度を含み、胚の正確な組成は未知であるため、実験データとシミュレーション結果のうちの任意のものとの間の合致するものを考察して、再現された摂動を示す。
【0124】
実験データに従う信頼可能な計算に必要な複数のパラメーターセットを試験した。モデルは多数のパラメーターを含むため、全ての可能なパラメーターセットを含むパラメーター空間は巨大である。驚くべきことに、20個ほどの少ないパラメーターセットであっても、800個のパラメーターセットを用いて得られた結果と同様な結果を生じるのに十分であった。これらの発見は、表3に要約され、本発明者らの分析において検出されるネットワークのトポロジーの、すなわちパラメーター不変の特徴に対し示唆を与える。
【0125】
実験データと比較した定量的結果にわたる概観が表2に与えられ、定性的概観が図15に与えられる。いくつかの摂動の効果が非常に良好に再現されている(gatae−KD、alx1−KD、及びbra−KD)。
【0126】
【表2】

【0127】
表2:異なる摂動実験に関するシミュレーション結果と実験データとの間の比較の概要。実験データとシミュレーション結果との間の合致するものの数が、列1において与えられ、全ての可能な合致するものの割合としての合致するものの数が第2の列に与えられている。表の底部に全ての摂動の平均が与えられている。
【0128】
図15は、これらの観測を裏付けるが、多くの場合に実験データに従って摂動に反応する遺伝子(pks、nrl、fvmo、alx1、及びbra)、及び、sm50、sm27、及びficolinのような、その予期されない挙動が、異なるパラメーターセットにより変化する多数のアップストリーム相互作用によって生じる可能性がある、実験データに従って反応することがほとんどない遺伝子をさらに示す。foxb、foxa、及びeveのような、重要な調節的遺伝子である他の遺伝子も、非常に多くの場合に実験データの再現に失敗する。遺伝子の双方のセットが実験データを不十分に再生するが、第2のセット内の遺伝子は、上述したものとは対照的に重要な調節的役割を有するため、これらの遺伝子をより差し迫って精緻さを欠いていると見なす。
【0129】
【表3】

【0130】
表3:用いられるパラメーターセット数及び実験データに対する全体的な一致度
要約すると、結果はets1及びsoxb1を伴う調節的相互作用の詳細な分析の必要性を指示し、一方で他の遺伝子の調節も同様にチェックする必要がある(foxb、foxa、及びeve)。
【0131】
全体として、モデルは実験データの42%を再生する。これはモデルが実験データとの一致度を増大するために精緻化を必要とすることを示す。この精緻化は、より多くの実験データに大きく依拠しなくてはならない。なぜなら、分析された遺伝子に対するモデル化されたMASO摂動の5920個の可能な効果のうちの小さな部分しか実験データに関連付けられないためである。
【0132】
ランダムネットワークとの比較
方法セクションにおいて説明したように、元のODEモデルにおけるエッジをランダムにスワップすることによって、比較可能な特徴を有する2つのランダムなネットワークが構築された。
【0133】
ランダム化されたネットワークは、サンプリングされたパラメーターセットを用いて、対照条件及び摂動条件の下でシミュレートされた。これらのパラメーターセットは元のモデルと同様であるが、前のセクション及び表3において説明された所見に起因して、ランダム化されたモデルの場合、元のモデルの場合の800個ではなく100個のパラメーターセットのみを用いた点が異なる。ランダム化されたモデルは、それらの一時入力においてのみ異なる3つの同一のサブモデルも含んでいた。ここで分析された2つのランダム化されたモデルは、実験データの20.15%及び23.5%のみを再現することが可能であった。3000個のエッジの代わりに30000個のエッジを同様の結果と切り換えた1つのランダム化されたモデルも検査した。
【0134】
これは元のモデルの内胚葉部分の範囲内であるが、いずれのランダム化されたモデルも、完全なモデルのPMC部分ほど有意にこの一致度を超えない(表4を参照のこと)。また、8個のランダム化されたネットワークのうちのいかなる3つのから成る組み合わせも、実験データとの全体一致度が25%を超えることはなかった。したがって、ランダム化されたネットワークの計算される一致度は、モデルの一般的なトポロジー特徴(サイズ及び接続性等)、実験データ、及び一時入力のみに依拠すると想定する。これは、特定の一時入力を用いてここで説明されたモデルのような比較可能な一般的な特徴のモデルと、ここで用いられた実験データとの間の一致度が、偶然のみで25%未満になると予期することができることを示す。
【0135】
【表4】

【0136】
表4:異なる胚領域に関するシミュレーション結果を用いて、実験データとの達成された一致度が示されている。この表に関して800個全てのパラメーターセットが用いられている。
【0137】
このため、偶然により予期される実験データとの一致度は、本発明者らのモデルを用いて観測された一致度の約半分であり、ランダムネットワークと比較して内胚葉系中胚葉ネットワークの妥当性が有意に改善されることを示す。
【0138】
本発明の方法を用いてウニ内胚葉系中胚葉ネットワークにおいて得られる結果が、刊行物クーン他(2009)においてさらに詳細に説明されている。クーン他の刊行物(2009)、及び本文脈では特にクーン他(2009)の「結果(Reults)及び検討」と題するセクションが、参照によりその全体が援用される。
【0139】
実施例2
アポトーシス
アポトーシスネットワークのin silico表現は、97個の微分方程式及び113個(全て未知)の反応速度パラメーターを含む。予測される濃度は、ノックダウン実験から得られた実験データと比較された。カスパーゼC3はsiRNAを用いてWi38細胞においてノックダウンされた。シミュレーションにおいてカスパーゼ3の実験的なノックダウンを反映するために、シミュレーションにおけるカスパーゼ3の初期濃度は、対照状況においてカスパーゼ3の濃度の20%に設定された。いずれの事例においても、すなわち対照及びノックダウンにおいて、アポトーシス経路の時間依存的挙動の400個のシミュレーションが実行された。結果が以下の表5に示される。
【0140】
【表5】

【0141】
表5:モデルが平衡状態に到達した後の、アポトーシス正常細胞及びカスパーゼ3ノックダウン(20%)細胞におけるアポトーシス経路のタンパク質及びタンパク質複合体の濃度。値は400個のシミュレーションの平均に対応する。
【0142】
カスパーゼ9及び切断カスパーゼ9の実験結果との比較。実験の設定において、カスパーゼ3のノックダウンによって、切断カスパーゼ9の量の有意な減少、及び(切断されていない)カスパーゼ9のわずかな増加が生じた。これらの結果は、本発明による方法を用いて行われた予測と良好に一致する。図5及び図6に示すように、カスパーゼ9から切断されたもの(from)の濃度が減少する一方、カスパーゼ9はそれに付随して増加する。
【0143】
実施例3
大型癌ネットワークに対するノックダウン効果の計算予測の生成
癌治療の効果を予測するための最小ネットワーク
癌治療のための最小相互作用ネットワークの開発及び注釈における第1のステップとして、癌に関連する遺伝子及びそれらの機能相互作用に関する大規模な情報を利用してきた。癌はおそらく、複数の遺伝子及び経路を伴う最も複雑な疾病のうちの1つであり(ビルド(Bild)他(2006)、ハナハン(Hanahan)及びワインバーグ(Weinberg)(2000)、ワインバーグ(2007))、たとえばアポトーシスの回避及び反成長シグナルに対する不感受性につながる、細胞生理学における激しい機能変化の症状と見なされる。これらの機能変化は、癌発病及び進行に伴う主要な分子及び経路と関連する。ほとんどの癌研究は、癌遺伝子及び腫瘍抑制遺伝子の突然変異に起因するこれらの経路の異常な活性の結果に集中してきた(キンズラー(Kinzler)及びフォーゲルシュタイン(Vogelstein)(1996))。細胞増殖及びアポトーシスの調節に重要なのは、細胞シグナル伝達ネットワーク、すなわち、広大なクロストークを示す複合ネットワークによる、成長シグナル及び死シグナルの認識及び統合である。経路間の正のフィードバックループが、不活性から永久に活性化された状態への遷移を引き起こす場合があり、これは連続した細胞増殖へとつながり、このため癌の発症に寄与する(キム(Kim)他(2007))。
【0144】
【表6】

【0145】
表6:癌における標的にされた治療薬剤。細胞表面受容体又は開始された経路の下流成分を標的にする異なる抗癌剤の選択。阻害実験は、阻害された単一の又は複数のモデル成分に関する。これらの成分は表9において説明される。
【0146】
現在の抗癌剤の大部分は、少なくとも名目的に、単一の標的に向けられるが、多くの場合、一緒に多くのオフターゲット効果を発生する。これらの薬剤のほとんどは、想定される標的及び関連する経路に対して直接的な阻害性効果を有する(表6)。たとえば、重要な標的経路はPI3K/AKTシグナル伝達経路である。なぜなら、該経路は細胞の成長及び存続の多くの局面にとって重大であるためである。該経路はゲノム異常によって影響を及ぼされ、これは多くの癌において経路の活性化へとつながる(ヘネシー(Hennessy)他(2005))。したがって、いくつかの癌実体、たとえば前立腺癌に関する前臨床及び臨床試験においてAKT阻害剤ペリホシンが用いられた(ヴァン アンメルセン(Van Ummersen)他(2004))。この薬剤は、おそらくAKTのPHドメインと相互作用することによって膜への局在を防ぐ。ウォルトマンニン及びLY294002は、in vitro及びin vivoでの抗腫瘍活性を示す(ヘネシー他(2005))。メトホルミン塩酸塩は、筋肉細胞及び脂肪細胞上のインスリン受容体の数及び/又は親和性を増大させ、受容体及び後受容体結合部位におけるインスリンへの感度を増大させる(NCI薬剤辞書http://www.cancer.gov)。ラパマイシン(WyethAyerst社ラパミューン)はAKTの下流で機能するmTOR(ラパマイシンの哺乳類標的)の特異的阻害剤である(ヘイ(Hay)及びソネンバーグ(Sonenberg)(2004))。mTOR阻害剤は、臨床試験において、乳癌及び他の固形腫瘍を有する患者に関して試験されている(チャン(Chan)他(2005)、イダルゴ(Hidalgo)及びロウィンスキー(Rowinsky)(2000)、ナガタ(Nagata)他(2004))。細胞において、エベロリムスはイムノフィリンFK結合タンパク質−12(FKBP−12)に結合し、mTOR、すなわち主要な調節キナーゼに結合してその活性化を阻害する免疫抑制複合体を生成する。mTOR阻害は、下方調節されたPTENによって生じるトラスツズマブ耐性を克服する試みにおいて研究されている。テムシロリムス(Wyeth社Torisel)は、キナーゼmTORの阻害剤であり、S6K1のリン酸化を阻止し(フェーヴル(Faivre)他(2006))、進行性腎細胞癌の治療に用いられる。ソラフェニブ(ネクサバール)は、RAF−キナーゼ、VEGFR−2、PDGFR、FLT−2、及びc−KITに対する経口マルチキナーゼ阻害剤であり(ストランバーグ(Strumberg)(2005))、血管形成及び腫瘍増殖を標的にする。これはインターフェロンα又はインターロイキン2に基づく治療に耐性を示す進行性腎細胞癌又は腎臓癌を有する患者の治療に適している。MEKは癌細胞の成長及び存続に伴うMAPK経路の重要な構成要素である。PD−325901は、この経路を遮断し癌細胞を破壊するように設計された新たな薬剤である。PD−0332991は、サイクリン依存性キナーゼ、特にCDK4を選択的に阻害し、網膜芽細胞腫(Rb)タンパク質のリン酸化を阻害することができる。Rbリン酸化の阻害は、Rb陽性の腫瘍細胞が細胞周期のS期に入るのを防ぎ、その結果DNA複製が抑制され、腫瘍細胞増殖が減少する。AEG35156はX連鎖アポトーシス阻害タンパク質(XIAP)、すなわち多くの細胞において過剰発現されるアポトーシスの阻害剤の細胞発現を選択的に阻止する。この化学物質は、細胞毒性薬剤と相乗的に作用して腫瘍細胞におけるXIAPの総レベルを低減し、アポトーシスに対する腫瘍細胞耐性を克服する。別の化合物、FJ9は、Frizzed−7 Wnt受容体とDishevelledのPDZドメインとの間の相互作用を妨害し、古典的(canonical)Wntシグナル伝達を下方調節し、腫瘍細胞成長を抑制する(フジイ(Fujii)他(2007))。エンザスタウリン塩酸塩は、ATP結合部位に結合して、プロテインキナーゼCβを選択的に阻害する。
【0147】
細胞の成長及び生存に不可欠な重要なシグナル伝達経路は、突然変異、増幅、転位を含むゲノム異常に起因してヒト癌において頻繁に活性化される。RAS−RAF−MEK−ERK、PI3K−AKT−mTOR、又はJAK−STATシグナル伝達経路のような成長及び生存経路を対象とする、ますます多くの合理的に設計された小分子阻害剤が、今や癌の治療の臨床試験に入ってきている(ヘネシー他(2005)、マッキュブリー(McCubrey)他(2008)、ヴァン アンメルセン他(2004))。それにもかかわらず、多くの阻害剤は、事前に知られていない「オフターゲット」によって生じる予期されない毒性に起因して、又は薬剤標的自体が脱制御に感受性のある複数の機能相互作用に関与するため、臨床試験において失敗する。加えて、標的にされた薬剤の臨床的失敗は、予期されないフィードバックループの存在、代替シグナル伝達経路の補完的上方調節、又は薬剤耐性突然変異によっても生じ、それらの全ては標的阻害の効果を回避し、腫瘍細胞を生存及び増殖させる(ブルヘルト他(2005))。
【0148】
したがって、これらの実施例が示すように、予測モデルは、複数の標的の複雑度及び経路間のクロストークに対処するために、関連する機能相互作用の多く(又は理想的には全て)を含むべきである。そのようなモデルは、患者の予期されない副作用又は不感受性を明らかにすることによって、新規な標的薬剤の開発の重要なサポートを提供することができる。しかしながら、これまで、癌過程の計算的モデリングは主に、RAF(キム他(2007))、AKT(アラウホ(Araujo)他(2007))、又はWNTシグナル伝達(キム他(2007))のような個々の下位経路に焦点を当ててきた。しかしながら、いくつかの分離した経路をより大きなフレームワークに統合することは、経路間のクロストークも捕捉し、薬剤作用の予測に不可欠である場合がある。この観点において、本発明者らは、次の段落で説明される癌に関連する過程の「最小」相互作用ネットワークに対する薬剤作用の効果をシミュレートすることによって計算による予測力を検査した。薬剤、それらの分子標的又は標的のセット、及びそれらが機能する細胞相互作用ネットワークに関する情報を凝集し、次のステップは阻害性効果をコンピューターにおいて翻訳することである。これは、薬剤をそれらの意図される薬剤標的タンパク質に関連させ、薬剤標的に関連付けられる単一の又は複数のモデル成分の阻害によって該薬剤標的の阻害の効果をシミュレートすること(to simulate)によって行われる(表6)。
【0149】
モデル集団及び注釈ワークフローの自動化
がん治療のための最小予測相互作用ネットワークを構成する主成分の選択は、主に生物学的機能に基づく。薬剤効果の計算モデル化は、経路スキーマ(図7A)を、モデル成分の濃度及びこれらの成分が関与する反応の反応速度に関する情報を保有することができるコンピューターモデルに変換することを必要とする。この過程は、適切なコンピューターオブジェクトの設計、反応の実装、これらの反応に対する反応速度則の割り当て、及びモデル分析を含む(図7B)。
【0150】
モデル集団、シミュレーション、及び分析をサポートする計算ツールがシステムズ生物学の基礎を構築する(ビリング(Wierling)他(2007))。いくつかのモデル化ツール、たとえばCellDesigner(フナハシ(Funahashi)他(2003))、E−Cell(トミタ(Tomita)他(1999))、Gepasi(メンデス(Mendes)(1993)、メンデス(1997))、及びその後継機Copasi(ホープス(Hoops)他(2006))、ProMoT/Diva(ギンケル(Ginkel)他(2003))、Virtual Cell(レーブ(Loew)及びシャーフ(Schaff)(2001)、スレプチェンコ(Slepchenko)他(2003))、及びSystems Biology Workbenchにおいて統合されるツール(ハッカ(Hucka)他(2002))が最近提案されている。これらのシステムのほとんどは、わずかな反応を用いて小さなモデルを手動で分析するように設計されている。しかしながら、関連する予測モデルの開発は、多数の成分、反応、及び反応速度パラメーターに関する情報を必要とする。このため、モデル化過程の自動化は、いくつかのステップにおいて、たとえば反応ネットワークの注釈、モデルの流れの可視化、微分方程式の数値積分、及びモデル分析において、大きなネットワークを扱う際、サポートを必要とする。
【0151】
この研究において、モデル化及びシミュレーションシステムPyBioS(ビリング他(2007))を用いてそのような手法を実証する。PyBioS(http://pybios.molgen.mpg.de)は、パスウェイデータベースへのインターフェースを提供することによってモデルの自動生成をサポートし、これによって関連する反応系への高速で自動化されたアクセスが可能になる。癌に関連する反応ネットワークに対する既存の知識の多くは、PyBioSシステムに直接取り込むことを可能である、BioCyc(カープ(Karp)他(2005))、KEGG(カネヒサ(Kanehisa)他(2006))、及びReactome(ジョシートープ(Joshi-Tope)他(2005)、ワストリック(Vastrik)他(2007))等のパスウェイデータベースに凝集される。
【0152】
本発明者らの分析は、3つの部分、すなわち癌に関連する経路を含むモデルの確立、たとえば薬剤の効果に起因するモデル成分の阻害の導入、及びモデル分析から構成される。薬剤標的阻害の効果をシミュレートする予測ネットワークの原型は、文献からの情報、及び癌の特質に焦点をあてた20個の異なるシグナル伝達経路に関する公共のパスウェイデータベースからの情報を組み合わせることによって編集された。現在のところ、ネットワークは1913個の個々の反応を有する767個の成分を含む(概要は表7を参照されたい)。
【0153】
【表7】

【0154】
表7:注釈付きヒト癌ネットワークに含まれるモデル成分のリスト
Cancer Gene Census:http://www.sanger.ac.uk/genetics/CGP/Census/
**ラス(Russ)及びランペル(Lampel)(2005)に基づく創薬可能遺伝子
【0155】
経路注釈は大方依然として手動で実行されるが、モデル化システムへの静的モデルのアップロードを記憶し容易にする自動化された注釈のためのツールが利用可能である。本発明者らの原型ネットワークの経路注釈は、シグナル伝達経路の情報を収集及び記憶する過程を自動化するReactome Curator Toolを用いて実行された(http://www.reactome.org)。全体ネットワークは、20個の異なる経路から構成される。それらの経路は、増殖因子(EGF、NGF、IGF−1、TGF−β)、細胞増殖(Wnt、Rb、ノッチ受容体、ヘッジホッグ)、サイトカイン(インターロイキン2、STAT−JAK)、炎症(トール様受容体)、アポトーシス(TNFα、FAS、TRAIL)、及び代謝調節(Gタンパク質共役型受容体)のような刺激によって活性化されるシグナル伝達カスケードを構成する。反応はPyBioSモデル化システムに取り込まれ、対応する数学的ODEモデルが反応系から自動的に作成された。
【0156】
モンテカルロ手法:不確実性に直面したモデル化
モデル化は、モデルサイズと予測精度とのトレードオフである。高い精度を有するモデルは、或る程度少ない数のモデル成分に基づいて非常に詳細な(of large detail)計算予測を生成する。しかしながら、これらの予測は、多くの場合にin vivo条件下で関連パラメーターを測定することが困難であることによって支障をきたし、これは関与する異なる経路間のクロストークを無視することによって悪化する。しかしながら、パラメーター適合戦略は、任意のそのような手法の一般的な難点を被るものの、不正確なモデルであっても、適合を生成するために十分なパラメーターを変形することができれば、一般的に極めて良好に適合することができる。特に、医学関連のモデルは、モデル精度に関する結果をもたらす多数のモデル成分を伴う可能性が高い。この観点において提案される戦略はこの目的のために設計される。モデルに関与する反応は、合成反応、複合体及び生成物形成、並びに不可逆質量作用反応速度
【0157】
【数6】

【0158】
によって説明される分解反応のような少ない数の異なる反応タイプから構成される。ただし、kは反応速度定数であり、Sはi番目の基質の濃度である。
【0159】
可逆反応はそれぞれ不可逆質量作用反応速度を用いる別個の前方向反応及び後方向反応によって説明される。複雑な形成及び解離の場合、100[a.u.]の解離定数を有する可逆質量作用反応速度が適用された。合成及び減衰反応が適宜導入された。適用されたモデルの速度則が図8に要約される。
【0160】
【表8】

【0161】
表8:異なるタイプの反応速度
本発明者らのモデル化手法において、同じネットワーク異なる状態、たとえば薬剤標的の異なるセットを不活性化したか又は不活性化していない生物ネットワーク間の差異の予測に焦点を当て、薬剤処理の効果の単純化されたモデルを表す。パラメーターのうちの多くにおける不確実性を補償するために、パラメーターベクトルの成分が適切な確率分布から選択され、利用可能な知識が反映される。ここで説明されるシミュレーション実行のセットでは特に、未知の反応速度パラメーターが対数正規分布から平均μ=2.5及び標準偏差σ=0.5でサンプリングされた。各パラメーターベクトル及び入力値の各ベクトルを用いて、通常対照状態及び阻害実験に対応する「処理された」状態を、該2つの状態の初期の差異が1つのパラメーター又は(parameteror)モデルパラメーターの小さなセットのみの変化に起因するようにモデル化した。異なる阻害実験のシミュレーションに関して、表9に列挙されているモデル成分が、標的タンパク質の完全な除去に対応する0.0[a.u.]の固定濃度を用いて初期化された。全ての成分の定常状態レベルが到達されるまで、異なるパラメーターベクトルを用いて対照及び処理シミュレーションが250回反復された(図8)。
【0162】
「処理された」状態と「処理されていない」状態との間のモデル挙動の差異が、個々のモデル成分ごとの定常状態濃度を比較することによって分析された。成功したシミュレーション実行にわたるモデル成分の2つの状態の最終定常状態値における差異を、コルモゴロフ−スミルノフ検定を用いて統計学的有意性に関して分析し(コノバー(Conover)(1971))、特別な治療によって影響を受けるモデル成分を特定する。
【0163】
薬剤標的阻害の系統的分析
提案されるフレームワークの検査として、表9に列挙されている異なる阻害実験に従って、無摂動状態におけるモデル成分の挙動を、処理された状態と比較した。薬剤標的のコンピューターによる阻害によって、いくつかのレベルで結果を予測する機会が与えられる。一方で、これは、定常状態濃度に従って統計学的に有意な変化のセットを計算することによって処理の感度の全体評価を可能にする。他方で、これは直接効果(治療における所望の効果)及び間接効果(治療の潜在的な望ましくない副作用)のような、経路成分における特別な変化の同定を可能にする。
【0164】
【表9−1】

【0165】
【表9−2】

【0166】
表9:抗癌剤の効果をシミュレートする阻害実験(表6を比較されたい)及びそれぞれの実験において阻害された関連付けられるモデル成分
【0167】
大域分析によって薬剤作用の差異をモニタリングする
図9は、全体統計を要約している。異なる阻害実験に関する複数の有意な変化によって表される摂動感度は或る程度変動する(図9A)。いくつかの阻害ドメイン、たとえば、AKT2酵素(モデル成分PIP3:リン酸化PKB)の活性体の阻害が、異なる阻害実験において単一の標的(IRS)又は複数の標的(mTOR、IRS、AKT、及びPI3K)の阻害によって60個より多くの異なるモデル成分に影響を及ぼすのに対し、他の阻害ドメインは非常に特異的であり、たとえばSTAT阻害は767個のモデル成分のうちの10個未満に影響を及ぼす。他方で、標的感度は非常に高く、ほとんどのモデル成分は阻害性効果に対してロバストである(図9B)。モデル成分の最も大きな部分(767個のうちの520個)は、阻害実験のいずれによっても影響を及ぼされない。最も有意な変化は、単一の阻害実験(73個)又は5つ未満の異なる阻害実験(2つ以上5個以下について138個)においてのみ観測される。モデル成分の小さな部分(35個)が、阻害実験のうちの多数(8個以上)における効果を示している。多数の薬剤標的阻害によって影響を及ぼされる成分は、シグナル伝達経路において中心的役割を果たすもの、たとえばGSK3β、及びそのリン酸化形式、若しくはWntシグナル伝達経路からのAPC及びアキシンに結合するGSKβ、又はmTORシグナル伝達の中心成分である。一般に、選択された薬剤及び阻害ドメイン(表9)は主成分分析(図9C)によって示されるように3つの異なるグループに入る。特に興味深いのはIRS、AKT2、PI3K、PDK1のグループ及びそれらの組み合わせである。なぜなら、これらの阻害実験はほとんどのモデル成分に非常に影響を及ぼすためである。
【0168】
AKTシグナル伝達を標的にすることによって、直接的な下流効果及び間接的な下流効果が明らかにする
大域分析を補完するものとして、モデル化手法は主要処理に関して詳細な予測を生成する。いくつかの阻害実験は、リン酸化GSK3β(ホスホ−GSK3β)、及びTSC1:阻害TSC2−1−P複合体のようなAKTシグナル伝達の直接的な不活性化又は間接的な不活性化を標的にする。これらの成分は、たとえば増殖及びアポトーシスとして癌関連細胞の読み出しに影響を与える。これらの介入点及び読出し点の概略図が図10に示されている。それぞれのモデル成分の阻害によって、薬剤副作用の経路クロストークの重大性を示す、直接的な効果及び間接的な効果が明らかになる。
【0169】
図11は、阻害実験の直接的な効果(左枠)又は間接的な効果(右枠)を示す選択された結果を示している。たとえば、リン酸化GSKβ(ホスホ−GSK3β)の定常状態濃度の低減は、AKT2阻害単体(AKT2)の直接的効果として、又はPI3K阻害と組み合わせて(AKTPI3K)見て取ることができる(図11A、図11G)。リン酸化GSKβ(ホスホ−GSK3β)の濃度における直接的な低減は、活性AKT2及びPI3K(AKT2、PI3K)の阻害に起因する。リン酸化GSKβ(ホスホ−GSK3β)の非リン酸化形式は、Wntシグナル伝達において重要な役割を果たす。PDK阻害(図11C)は、AKT2のリン酸化において直接的な効果を有し、結果としてPIP3:リン酸化−PKB複合体の下方調節をもたらす。PIP3がセリン/トレオニンキナーゼPDK1及びAKT2を原形質膜に入れることが既知である。ここでAKT2はPDK1媒介リン酸化によって活性化される。IRS阻害実験において、PI3Kの阻害は、インスリン受容体基質(IRS)のリン酸化(活性)形態の阻害に起因して直接的な効果と見なされ、PI3K(図11E)の主要活性剤は、それに続く、ホスホ−IRS:PI3K複合体の定常状態濃度の低減をもたらす。
【0170】
補間的に、モデル化戦略は多くの間接的な効果を同定する。S6K1は、小さな(40−S)リボソームサブユニットの成分であり、このサブユニットがリボソーム形成及びしたがってタンパク質合成に関与することを可能にする。S6K1のリン酸化体が阻害実験AKT及びAKTP3Kにおいて下方調節されることがわかった(図11B、図11H)。S6Kリン酸化の阻害は、AKTシグナル伝達経路におけるPI3Kカスケードの下流成分及びmTORシグナル伝達に対する下流効果に起因する。GSK3βのリン酸化(図11D)は、PDK(阻害実験PDK1)によって間接的に制御される。なぜなら、AKTの阻害に起因して、PDKの阻害がホスホ−GSK3βの下方調節につながるためである。さらなる間接的な効果として、活性化されたIRSの阻害によって(阻害実験IRS)、TSC2のリン酸化の低減につながる(図11F)。
【0171】
阻害実験の予測効果の共発現クラスターは文献による知識との一致を示す
多くのモデル成分が、阻害実験の集団全体にわたって同様の発現パターンを示しており、薬剤作用のこれらの共発現クラスターのうちのいくつかを、文献からの以前のデータによって説明することができる。図12は、AKTを伴う阻害実験によって影響を及ぼされるモデル成分の特殊な例を示している。
【0172】
AKT活性化は、プレクストリン相同ドメイン(PHD)のホスファチジルイノシトール−3,4,5−三リン酸(PIP3)への結合、及びその後に続く調節アミノ酸セリン473及びトレオニン308のリン酸化によって開始される膜への局在によって駆動される(ヴィヴァンコ(Vivanco)及びソイヤーズ(Sawyers)(2002))。AKTと原形質膜との病的会合はAKTを癌につなげる共通関連性(common thread)である。AKTの活性体の阻害に基づく薬剤効果(阻害実験AKT2、DvlAKT2、AKTPI3K、IRSAKTPI3K、mTORIRSAKTPI3K、表9を参照されたい)に関して、下流成分の下方調節を同定することができる(図12)。さらに、阻害実験PI3K及びPDK1は、AKT阻害に基づく阻害実験と比較して同様の阻害された成分を示している。この観測は、AKTがPI3K及びPDK1の効果の一次下流仲介物質として同定される文献データ(ヘネシー他(2005))と一致する。
【0173】
IRS阻害実験における成分のうちのいくつかは、AKT2と同様の挙動を示す。しかしながら、IRS阻害実験は、3つの特異成分(RAC−βセリン/トレオニンプロテインキナーゼ[AKT2]、そのPKB調節因子[PKB:PKB調節因子、AKT2はPKBの別名である]との複合体、及びPI3K)の上方調節によって特徴付けられるのに対し、ホスホ−IRS:PI3K複合体は下方調節される。成分eEF2K−P、elF4G−P、リン酸化40S小リボソームタンパク質、elF4B−P、TSC1:阻害されたTSC2−1−P、S6K1−P、活性化されたmTORC1、阻害されたTSC2−1−P、及びリン酸化AKT複合体の阻害は、AKT阻害に起因する。AKT1によるTSC1:TSC2複合体のリン酸化の結果、TSC1:阻害されたTSC2−1−P複合体、及びプロテオソーム経路を通じたその後の分解がもたらされる。PDE3Bリン酸化の下方調節はAKT阻害に起因し、PI3K阻害実験においても注目される。これによって、ウォルトマンニンを通じたPI3K阻害が、PDE3Bのリン酸化及び活性化、並びにインスリンの抗脂肪分解効果を阻害するという報告(ラーン(Rahn)他(1994))が裏付けられる。さらに、PI3Kによる活性化に起因して、mTORがS6Kをリン酸化及び活性化し、mRNAの翻訳を増加させるため、S6Kのリン酸化におけるPI3Kの影響が知られている(ローデス(Rhodes)及びホワイト(White)(2002)、サルティエル(Saltiel)及びカーン(Kahn)(2001))。
【0174】
感度分析を用いたネットワーク依存性のモニタリング
感度分析を用いて、異なるモデル成分間の関係及び相互作用を理解する。これらの相互作用を定量化するためのいくつかの方法が提案された。たとえば、代謝制御分析(MCA)は基本系における無限小の摂動に対するモデル成分の感度を研究する(クリップ(Klipp)他(2005))。チョー(Cho)他(2003)において、マルチパラメトリック感度分析(MPSA)が紹介された。この方法によって、数学モデルがモンテカルロシミュレーションに基づく手法に関して非常に感受性があるパラメーターの同定が可能になる。チアン(Jiang)他(2007)において、時間依存定量的制御値を計算することによって系の局所挙動が分析された。時間におけるパラメーター感度軌跡を用いて、反応速度パラメーターにおける無限小の変化に対する代謝産物の感度が求められた。
【0175】
特定のモデル成分の100%の阻害に加えて、より小さい入力変動におけるモデル成分の挙動を調査するために、感度分析を実行した。一度に29個のモデル成分を用いて小さな阻害(10%)が実行され、以前のセクションで説明した異なる阻害実験において分析された。これらのモデル成分のそれぞれについて、その定常状態値の90%に固定された選択された成分を除いた対照状態の定常状態値を用いて阻害実験のODEモデルを初期化した。続いて、10%阻害に関する定常状態の変化が計算された。10%低減の場合の各モデル成分の定量化された感度が、
【0176】
【数7】

【0177】
によって与えられる。ただし、i=1,...,nは全ての成分にわたるインデックスであり、j=1,...,mは阻害された標的のインデックスである。したがって、
【0178】
【数8】

【0179】
は制御実行におけるインデックスiを有する成分の定常状態濃度を表し、
【0180】
【数9】

【0181】
はインデックスjを有する薬剤標的の10%阻害モデルにおける成分iの定常状態濃度を表す。
【0182】
結果としての感度値のクラスター化によって、29個の阻害実験のセットにわたって同様の感度パターンを示すモデル成分のグループが明らかになる。
【0183】
図13は、AKT2(別名RACβセリン)に対し高い感度を示すモデル成分の選択されたクラスターを示している。非活性AKT2をわずかに低減することによって活性AKT2(PIP3:リン酸化PKB複合体)、並びにその後の標的、すなわちTSC1:阻害されたTSC2及びホスホ−GSK3βにおける有意な変化がもたらされる。感度分析によって、GSK3βのリン酸化体(ホスホ−GSK3β)も、様々な他の薬剤標的(たとえばSRC、GSK3、PIP3複合体、PDK1、PI3K)における小さな変化に感受性があることが明らかになっている。
【0184】
実施例1で引用される引用文献
【0185】
【表10−1】

【0186】
【表10−2】

【0187】
【表10−3】

【0188】
実施例3で引用される引用文献
【0189】
【表11−1】

【0190】
【表11−2】

【0191】
【表11−3】

【0192】
【表11−4】

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

【0195】
[EI]が酵素/阻害剤複合体の濃度であり、[E]及び[I]がそれぞれ遊離型酵素及び阻害剤の濃度(the concentrations of the concentrations of)である所与の反応に関して、定常状態における阻害された酵素と遊離型酵素との間の比は、解離定数によって以下のように与えられる。
【0196】
【数11】

【0197】
ただし、k及びk−1は、それぞれ解離反応の反応速度定数及び会合反応の反応速度定数(the kinetic constants of the dissociation respective association reaction)であり、Kは解離定数である。
【0198】
GEが遊離型酵素及び阻害された酵素の総キナーゼ濃度である場合、遊離型酵素の濃度yを
【0199】
【数12】

【0200】
によって計算することができる。
【0201】
阻害剤の濃度が各キナーゼの濃度と比較して高い場合、この式は、薬剤によって阻害され得る任意のキナーゼに適用することができる。
【0202】
阻害剤の全体濃度が細胞内のキナーゼの濃度と同じ範囲内にある場合、個々のキナーゼに結合される阻害剤の量も考慮に入れなくてはならない。

【特許請求の範囲】
【請求項1】
生物ネットワークの反応速度モデルを生成する、コンピューター実装される方法であって、
(a)ネットワークトポロジーを選択することであって、該トポロジーのノードは生物学的実体を表し、該トポロジーのエッジは該実体間の相互作用を表す、ネットワークトポロジーを選択すること、
(b)前記相互作用に反応速度則及び反応速度定数を割り当てること、並びに
(c)前記生物学的実体に開始濃度を割り当てること、
を含み、
(i)前記反応速度定数の一部分、及びそれとは独立して前記開始濃度の一部分が実験データであり、
(ii)前記反応速度定数の残りの部分、及びそれとは独立して前記開始濃度の残りの部分がランダムに選択される、方法。
【請求項2】
生物ネットワークにおいて、生物学的実体の濃度を時間の関数として予測する方法であって、
請求項1に記載の方法によって生物ネットワークのモデルを生成すること、及び
(d)連立微分方程式を解くことであって、該微分方程式は前記生物学的実体の濃度の時間依存性を規定する、連立微分方程式を解くこと、
を含み、それによって前記濃度を時間の関数として取得する、方法。
【請求項3】
前記反応速度定数の前記残りの部分は確率分布から選択され、それとは独立して、前記開始濃度の前記残りの部分は確率分布から選択される、請求項1又は2に記載の方法。
【請求項4】
前記分布は、対数正規分布、一様分布、指数分布、ポアソン分布、二項分布、コーシー分布、ベータ分布、又はガウス分布である、請求項3に記載の方法。
【請求項5】
前記反応速度定数の前記残りの部分をランダムに選択し、前記開始濃度の前記残りの部分をランダムに選択すると共に、その後前記連立微分方程式を解くことが反復的に実行される、請求項2〜4のいずれか1項に記載の方法。
【請求項6】
(e)前記方法は1つ又は複数のさらなる生物ネットワーク上で付随して実行され、
(f)生物学的実体の前記濃度は、選択された時点において前記生物ネットワーク間で交換される、
請求項2〜5のいずれか1項に記載の方法。
【請求項7】
更に前記反応速度則の一部分が知られており、該反応速度則の残りの部分はランダムに選択される、請求項1〜6のいずれか1項に記載の方法。
【請求項8】
前記生物学的実体の前記濃度は、野生型と比較して摂動されている、請求項1〜7のいずれか1項に記載の方法。
【請求項9】
前記摂動された濃度は実験的に求められた濃度である、請求項8に記載の方法。
【請求項10】
前記濃度は、ノックダウン若しくは過剰発現の実験、又は病的状態、又はin vitro若しくは細胞の薬剤阻害実験において求められる、請求項9に記載の方法。
【請求項11】
初期状態は、
(a)生物学的実体の実験的に求められた濃度、及び/又は
(b)実験的に求められた突然変異データ
を含む、請求項1〜10のいずれか1項に記載の方法。
【請求項12】
請求項2〜11のいずれか1項に記載の方法の統計学的有意性を求める、コンピューター実装される方法であって、該統計学的有意性を求める方法は
(a)請求項2〜11のいずれか1項に記載の方法を実行すること、
(b)ステップ(a)において取得された生物学的実体の濃度と、同じ生物学的実体について実験的に求められた濃度との間の一致度を求めること、
(c)前記生物ネットワークのトポロジーをランダム化すること、
(d)ステップ(b)において取得した前記ランダム化された生物ネットワークに対し、請求項2〜11のいずれか1項に記載の方法を実行すること、
(e)ステップ(d)において取得した生物学的実体の濃度と、同じ生物学的実体について実験的に求められた濃度との一致度を求めること、
(f)ステップ(b)において取得した結果をステップ(e)において取得した結果と比較することであって、ステップ(b)における一致度の方が高いことは、請求項2〜11のいずれか1項に記載の前記方法が、実験的に求められた濃度を偶然よりも良好に予測可能であることを示す、比較すること、
を含む、方法。
【請求項13】
(c)に従ってランダム化することは、前記ネットワークのエッジをスワップすることによって有効にされる、請求項12に記載の方法。
【請求項14】
前記実体は生体分子であり、好ましくは遺伝子を含む核酸、タンパク質を含む(ポリ)ペプチド、小分子、並びにそれらの複合体及び代謝産物から選択される、請求項1〜13のいずれか1項に記載の方法。
【請求項15】
前記モデルは、境界条件、好ましくは生理学的状態を表す境界条件を含む、請求項1〜14のいずれか1項に記載の方法。
【請求項16】
生物ネットワークの部分的に未知のパラメーターを求める、コンピューター実装される方法であって、該パラメーターは、ネットワークトポロジー、反応速度則、反応速度定数、及び/又は開始濃度から選択され、該方法は、観測される特性と予測される特性との差異を最小にすることであって、該予測される特性は、請求項2〜11のいずれか1項に記載の方法によって予測されるような濃度を含む、最小することを含む、方法。
【請求項17】
1つ又は複数の実験を選択する、コンピューター実装される方法であって、
(a)請求項2〜11のいずれか1項に記載の生物学的実体の濃度を予測する方法を実行することによって複数の実験をin silicoで実行することであって、該実行することは未知のパラメーターの異なる選択を用いて実験ごとに反復して行われ、該パラメーターは、ネットワークトポロジー、反応速度則、反応速度定数、及び/又は開始濃度から選択される、実行すること、並びに、
(b)前記複数の実験から、前記異なる選択に依拠して、前記生物学的実体の濃度を予測する方法が予測濃度の最大の分散をもたらす1つ又は複数の実験を選択すること、
を含む、方法。
【請求項18】
請求項1〜17のいずれか1項に記載の方法を実行するように構成されるコンピュータープログラム。
【請求項19】
請求項18に記載のプログラムを含むコンピューター可読データキャリア。
【請求項20】
データ処理装置であって、請求項1〜17のいずれか1項に記載の方法を実行する手段を備えるか、又は該データ処理装置にインストールされる請求項18に記載の前記プログラムを有する、データ処理装置。

【図1】
image rotate

【図2】
image rotate

【図3】
image rotate

【図4】
image rotate

【図5】
image rotate

【図6】
image rotate

【図7】
image rotate

【図8】
image rotate

【図9】
image rotate

【図10】
image rotate

【図11−1】
image rotate

【図11−2】
image rotate

【図12】
image rotate

【図13】
image rotate

【図14】
image rotate

【図15】
image rotate


【公表番号】特表2012−502335(P2012−502335A)
【公表日】平成24年1月26日(2012.1.26)
【国際特許分類】
【出願番号】特願2011−525474(P2011−525474)
【出願日】平成21年9月3日(2009.9.3)
【国際出願番号】PCT/EP2009/007223
【国際公開番号】WO2010/025961
【国際公開日】平成22年3月11日(2010.3.11)
【出願人】(301078836)マックス−プランク−ゲゼルシャフト・ツア・フェルデルンク・デア・ヴィッセンシャフテン・エー.ファウ.・ベルリン (1)