JP3443039B2 - ネットワーク推定方法および装置 - Google Patents

Description

【発明の詳細な説明】

【０００１】

【発明の属する技術分野】本発明は、互いに関連し、干
渉し合う別々の要素のネットワークとしてモデル化でき
るシステムを、そのネットワークの挙動や状態を表す情
報から逆に推定することに関する。特に、生物学的現象
における遺伝子の発現データ、タンパク質の濃度データ
等から、その現象を引き起こしている遺伝子の制御ネッ
トワークや酵素やタンパク質の反応を示す代謝ネットワ
ークの構造等を推定することに関する。

【０００２】

【技術の背景】最近の分子生物学の急激な進歩によっ
て、主要な生物に関する多くの情報が蓄積されつつあ
る。その典型例が、酵母、線虫、ショウジョウバエであ
る。例えば、線虫（C.elegans)では、その全細胞系譜の
推定、神経系の結合の同定などがすでに行われており、
昨年、ＤＮＡの全塩基配列の決定が発表された。酵母、
大腸菌、マイコプラズマ等のモデル生物においても、Ｄ
ＮＡの全塩基配列がすでに決定されており、究極の目標
であるヒトのＤＮＡの全塩基配列も数年以内には終了す
ると思われる。しかしながら、配列が分かるということ
と、遺伝子やその機能、さらには遺伝子間の相互作用が
分かるという間には、大きなギャップがある。

【０００３】このような遺伝子の機能や遺伝子間の相互
作用の決定は非常に困難な仕事である。現在、分子生物
学の多くの研究は、注目する現象に関わる遺伝子および
その転写産物の特定と、それらのカスケードの決定に重
点を置いている。このような遺伝子間の相互作用の理解
を行う段階において、複雑な相互作用を人間の直感のみ
を用いて行うことは極めて難しい。これは、データが大
量にあることやこれに関わる相互作用等の関係も色々と
考えられるからである。我々は、この困難をコンピュー
タという計算能力をもった道具を導入することで乗り越
えようとしている。

【０００４】この考えを図式したのが図１である。これ
は、コンピュータを利用する部分（バーチャル）と、実
際に生物学的に実験する部分（リアル）との関係を示す
図である。図１において、まず、解明しようとする現象
を特定すると、現在知られている生物学的知識（図１
フェーズ−I Ａ）から、モデルをコンピュータ上に実
装する（図１ フェーズ−I Ｂ）。また、すでに仮説
が提示されている場合には、その仮説を実装する。実装
されたモデルを使ってシミュレーションを行い（図１
フェーズ−I Ｃ）、観測されているデータとの整合性
を検証する（図１ フェーズ−I Ｄ）。シミュレーシ
ョンの結果と、観測データと一致しない場合は、２つの
可能性が考えられる。まず、シミュレーションが正確で
ない場合である。これは、シミュレーションの精度をあ
げて、確実なものとすることで解決する。次に考えられ
るのは、モデルが不完全または間違っている場合であ
る。この場合は知られている実験データに整合するモデ
ルを構築する必要がある。その後に、できるならば、こ
のモデルを用いて、コンピュータ上で実験を行う（図１
フェーズ−II Ｆ）。フェーズIIの後、又はフェーズ
IIができないなら、直接、仮説的に決めたモデル等が正
しいか実際に生物学的な実験を開始する（図１フェーズ
−IIIおよびIV）。このようにして、「生物のリバース
・エンジニアリング」を行うことができる。このような
コンピュータ上におけるシミュレーションにより、仮説
的に求めたモデルの絞り込みを必要としている分野は、
上述の生物の分野に限らない。

【０００５】

【発明が解決しようとする課題】本発明の目的は、コン
ピュータにより、仮説的に求めたモデルの絞り込みを効
率良く行うことである。例えば分子生物学においては、
ＤＮＡチップやマイクロアレイ、さらには、ＰＣＲなど
を利用して得られた遺伝子産物やタンパク質の量の時間
的変化に関するデータ、電気生理学的実験や電圧感受性
色素などを用いた神経回路網の活動に関するデータ等か
ら、その背後にある遺伝子やタンパク質、酵素の相互作
用のネットワーク、および、神経細胞のネットワークで
表されるモデルを、コンピュータを用いて推定すること
である。与えられたデータから考えられるモデルの総数
は莫大な数である。従来は、これを人間の勘に頼ってあ
たりを付け、生物学的な実験により確かめていた。本発
明は、与えられたデータから、類似したデータを生成す
る可能性のあるモデル（ネットワーク）を、コンピュー
タを用いて推定することにより、確からしいモデルの数
をある程度まで絞り込もうとするものである。これは、
上述の図１に示した場合では、フェーズIに関するもの
である。これらの相互作用のモデルの推定（ネットワー
クの推定）は、疾病の原因遺伝子の特定や、投薬の効果
の推定等に必要である。また、未知の遺伝子や遺伝子産
物の推定にも用いることができる。

【０００６】

【課題を解決するための手段】上記目的を達成するため
に、本発明は、あるデータを生成する要素を少なくとも
１つ用いて、互いに関連している要素間の関係をネット
ワークとして記述でき、前記ネットワークにより生成さ
れたデータが与えられたときにおける、前記与えられた
データから、前記データを再現できるネットワークをコ
ンピュータ・システムを用いて推定する方法であって、
前記ネットワークは、ネットワーク構造と、それに対応
するパラメータ・セットと、ネットワーク構造およびパ
ラメータ・セットから生成されるデータと前記与えられ
たデータとの適応度とからなるトリプレットとして表現
され、前記与えられたデータを実現することが想定され
る、部分的に既知のネットワーク構造を基にネットワー
ク構造を生成し、それに対応するパラメータ・セットお
よび適応度を作成し、前記適応度を利用して最適化した
複数の候補トリプレットを第１記憶手段に格納するステ
ップと、前記与えられたデータとは異なる、ネットワー
ク構造の一部が変異したネットワーク構造から生成でき
る他のデータを用いて、前記第１記憶手段に格納された
候補トリプレットから適切な候補トリプレットを絞り込
み、第２記憶手段に格納するステップとを備えることを
特徴とする。

【０００７】前記候補トリプレットを第１記憶手段に格
納するステップの前記適応度を利用した最適化は、作成
したネットワーク構造からＮ個を選択して、各々に最適
化したパラメータ・セットをＭ個適合してＮ×Ｍ個のト
リプレットを作成し、さらに、該Ｎ個のネットワーク構
造から新たなＮ個のネットワーク構造を生成して、各々
に最適化したパラメータ・セットをＭ個適合してＮ×Ｍ
個のトリプレットを作成し、作成した２Ｎ×Ｍ個のトリ
プレットから適応度の高いＰ個のトリプレットを選択す
ることもできる。前記候補トリプレットを第１記憶手段
に格納するステップの前記適応度を利用した最適化は、
さらに、前記Ｐ個のトリプレットの近傍を探索して、よ
り適合度の高いトリプレットを見つけて置き換えること
もできる。このように広域探索をしたり、さらに局所探
索をハイブリッドで用いることにより、効率的に最適な
トリプレットを選択することができる。前記候補トリプ
レットを第１記憶手段に格納するステップの前記適応度
を利用した最適化は、遺伝的アルゴリズム（ＧＡ）、シ
ミュレーテッド・アニーリング（ＳＡ）、及び／又は、
山登り法を用いることができる。これにより、与えられ
たデータを再現するための、ネットワーク構造に対する
最適なパラメータを早く推定することができる。このよ
うにして、最終的に、効率的に比較的少数の最適化され
た候補トリプレットを得ることができる。上述のネット
ワーク推定方法を実行する装置および、上述のネットワ
ーク推定方法をコンピュータ・システムに実現させるこ
とができるプログラムを格納した記録媒体も本発明であ
る。

【０００８】

【０００９】

【００１０】

【発明の実施の形態】本発明の実施形態を、図面を参照
して詳細に説明する。本発明は、例えば、特定の生物学
的現象における遺伝子の発現データ、タンパク質の濃度
データ等のプロファイルから、その現象を引き起こして
いる遺伝子の制御ネットワークや、酵素やタンパク質の
反応を示す代謝ネットワークの構造と反応の強さや速度
などを決定するいろいろな変数を推定する場合に適用す
ることができる。この場合を例に以下説明する。

【００１１】（遺伝子ネットワーク）本発明の実施形態
を、図２に示すような遺伝子の相互関係を示しているネ
ットワークを例に説明する。図２（ａ）は、遺伝子をノ
ードとして、その遺伝子が発現する物質等の発現に対し
て活性を示す関係と、抑制を示す関係を示している。そ
れをネットワーク構造として図示したものが図２（ｂ）
である。さて、与えられたターゲット・プロファイル
（発現プロファイル）に対する遺伝子のネットワーク
は、構成要素（ＤＮＡ，ＲＮＡ，タンパク質等）間の反
応関係を示すネットワーク構造（トポロジー）、ネット
ワーク内部での反応や状態のモデルを記述する際に必要
な変数の集合であるパラメータ・セット、および、ター
ゲット・プロファイルとの類似度を示す適応度の３つを
単位とした、「トリプレット」で表すことができる。

【００１２】まず、遺伝子のネットワーク構造の例につ
いて説明する。構成要素（ノード）間の反応規則は、簡
単化すると、Ｎ×Ｎの結合重み行列Ｗ（パラメータ・セ
ット）および結合行列Ｃ（ネットワーク構造）で表現さ
れる。結合重み行列Ｗの値ｗは、例えば、＋１．０と−
１．０の間の実数をとる。ｗの負の値は抑制を示し、正
の値は活性を示す。結合行列Ｃでは、要素ｃ_pqの０と１
の２値で、ノードｐとノードｑの間の結合関係を示す。
０はノード間に結合関係がないことを意味し、１は結合
関係があることを示す。また、Ｎノードを有する遺伝子
ネットワーク構造のある時間ｔにおける発現状態は、Ｎ
次元空間のベクトルｘ（ｔ）で表現することができる。
ベクトルｘの各要素は、ｘ_i（ｔ）（ｉ＝１，…，Ｎ）
で表される。そして、Ｓ_i（ｔ）を、要素ｉ（ｉ＝１，
…，Ｎ）の発現型の濃度であり、ｈ_iを活性要素ｉの活
性レベルを決定するしきい値であるとすると、各要素
は、例えば、次のように定義することができる。この関
数は、様々に定義することが可能である。

【数１】

【００１３】（実施形態の処理の概略）このような遺伝
子ネットワークの場合は、複数の発現物質の濃度の時間
的変化で示される発現プロファイルが与えられて、この
発現プロファイルを発現する遺伝子ネットワークの候補
を探し出すことが本発明の目的である。例えば、図３
（ａ）に示すような発現プロファイルを与えて、図３
（ｂ）に示すような遺伝子ネットワークの候補を探し出
すのである。なお、この図３（ａ）に示した発現プロフ
ァイルは、野生型の遺伝子構造に対する発現プロファイ
ルである。最初に用いる発現プロファイルは、これに示
されるような野生型であることが望ましい。ネットワー
クの結合（リンク）部分の数字は、上述における結合重
み行列Ｗの値ｗの絶対値を示している。なお、結合重み
行列Ｗの値ｗの符号は、結合の活性か抑制かにより正負
が決められる。また、遺伝子Ｂがない突然変異（ミュー
タント）のプロファイル、遺伝子Ｃのない突然変異（ミ
ュータント）のプロファイルがそれぞれ図３（ｃ）、図
３（ｄ）に示されている。このプロファイルの用いられ
方については後で説明する。

【００１４】本実施形態は、図４の概略フローチャート
に示すようなプロセスにより行われる。図４において、
まず、対象である発現プロファイルを実現することがで
きるネットワークの構造を想定し、そのネットワーク構
造（トポロジー）を生成して、トポロジー・プール３０
０を作成する（Ｓ１０２）。想定されるネットワーク構
造の数が少ない場合は、想定される全てのネットワーク
構造を生成して、トポロジー・プール３００に格納する
ことができるが、想定されるネットワーク構造の数が多
い場合は、トポロジー・プール３００に格納するトポロ
ジーに対して何らかの制限を設ける必要がある。なお、
必ずしも、トポロジー・プール３００にトポロジーを生
成して格納する必要はない。想定されるネットワーク構
造により構成されるトポロジー空間を定義することがで
きれば、この定義されたトポロジー空間から任意のトポ
ロジーを生成するアルゴリズムを用いることにより、ト
ポロジー・プール３００からトポロジーを取り出すのと
同様のことができる。

【００１５】生成したトポロジーのそれぞれに、発現プ
ロファイルからのパラメータを適応度を利用して適合さ
せる（Ｓ１０４）。パラメータは、上述の結合重み行列
Ｗの値としきい値である。また、適応度は、このパラメ
ータの確からしさ（どれだけ合致しているか）である。
この適応度は、評価関数を用いて計算されるが、例え
ば、この評価関数として最小自乗誤差などがある。

【００１６】ネットワーク構造（トポロジー）に対し
て、パラメータと適応度を組み合わせられたものは、前
述したように、トリプレット（トポロジー、パラメー
タ、適応度の３点セット）という。処理結果のトリプレ
ットをトリプレット・プール４００に格納する。トポロ
ジーを選択してパラメータを適合させ、適応度を算出す
るときに、ＧＡ（遺伝的アルゴリズム）やＳＡ（シミュ
レーテッド・アニーリング：焼きなまし法）等の手法を
用いて、与えられている発現プロファイルに対して最適
化する。

【００１７】ここまでは、処理対象である標本（トリプ
レット）を抽出する段階である。これから、目的に適合
した標本（トリプレット）をトリプレット・プール４０
０から絞り込む段階となる。この処理には、例えば、突
然変異アナリシスを用いている（Ｓ１０６）。この処理
により、標本（トリプレット）が絞り込まれて、目的と
する候補トリプレットの集合（候補トリプレット・プー
ル５００）ができる。

【００１８】（ネットワーク抽出過程）ネットワーク構
造の生成（Ｓ１０２）の処理は、例えば、最大でいくつ
のノードを含むネットワークを想定するか定める。その
ネットワークの構造を、予備知識無しから生成するか、
部分的にその構造が分かっているネットワークを基に生
成する。生成したネットワークをトポロジー・プール３
００に格納する。この場合に用いられるネットワークの
分かっている構造には、例えば、ファンアウトの制限が
ある。これを説明するのが図５である。ファンアウトと
は、ノードから他のノードへの働きかけ（活性、抑制）
の数である。図５（ａ）は、ＡからＢ〜Ｅへ、４個のフ
ァンアウトがあることを示している。このファンアウト
の数には制約があると考えられ、この制約は、図５
（ｂ）に示す例のように、確率で示されている。図５
（ｂ）では、ファンアウトの数の確率は、２付近で極大
となり、５以上のファンアウトの数の確率は無視できる
程度であることを示している。なお、この制約も様々に
定義することが可能である。また、すでに同定されてい
る結合や未結合が分かっている場合は、その結合以外の
部分に対してネットワーク構造を発生させる。ネットワ
ーク構造を発生させるときに、ある部分の結合構造を仮
定して、発生させることも可能である。このような色々
な構造に対する制約的な規則を適用して、トポロジーを
作成することもできる。

【００１９】トポロジー・プール３００に格納されてい
るネットワーク構造（トポロジー）の全てにパラメータ
や適応度を適合して、トリプレット・プール４００に格
納するか、抽出して格納するかは、生成され得るネット
ワークの数による。生成され得るネットワークの数が少
ないとき、すなわち、ネットワークのノードの数が少な
い場合は、あり得る全てのネットワークを網羅的に処理
することが可能であるので、生成されたトポロジーに全
てパラメータおよび適応度を適合させてトリプレット・
プール４００に格納する。しかし、生成されるネットワ
ークの数が膨大で、全数を検討することができない場合
は、トポロジーを選択してパラメータや適応度を割り当
て、トリプレット・プール４００に格納する。

【００２０】（トポロジーが多数存在する場合）さて、
トポロジーの数が膨大で、トポロジーの選択が必要な場
合における処理を、図６のフローチャートと図７および
図８を用いて、もっと詳しく説明する。図７および図６
のフローチャートを参照する。まず、最初にいくつのト
ポロジーを選択するか定める。この場合は、Ｎ個のトポ
ロジーをランダムに選択するとする。そして、図４のト
ポロジー・プール３００（又は、トポロジー空間）か
ら、ランダムに１つトポロジーを取り出す（Ｓ２０
２）。これに対して、与えられた発現プロファイルを実
現できる発現パラメータを、ＧＡ（遺伝的アルゴリズ
ム）／ＳＡ（シミュレーティド・アニーリング）などを
用いて適合させる（Ｓ２０４）。このとき、各発現パラ
メータによる発現プロファイルとの適応度の高いものか
ら、異なるパラメータ・セットをＭ個とり、トリプレッ
ト・プール４００に格納する。これをＮ個のトポロジー
について繰り返す。繰り返しにおいて、ランダムにトポ
ロジーを選択するとき（Ｓ２０２）、前に選択したトポ
ロジーとは異なる構造のものを選択する。繰り返しが終
了した段階（Ｓ２０６でＹＥＳ）では、Ｎ×Ｍのトリプ
レットがトリプレット・プール４００に格納されてい
る。なお、上述では、ランダムにトポロジーを選択して
いるが、他の確率的サンプリングを用いてもよい。ま
た、異なるパラメータ・セットをＭ個とるとき、パラメ
ータ・セット同士の類似度を計算して、その類似度が極
めて近いとき（あるしきい値より小さいとき）は、同じ
として選択しないようにしてもよい。これは、同じよう
な解の近傍を複数選択しても意味が無いことが多いの
で、代表値のみを選択していることを意味している。

【００２１】次に、トリプレット・プール４００に格納
されているＮ個のトリプレットのトポロジー組みかえを
したＮ個のネットワーク構造（トポロジー）をＧＡ／Ｓ
Ａ等により得る（Ｓ２０８）。これらのＮ＋１〜２Ｎの
得られた各トポロジーに対して、上述と同様に、パラメ
ータ・セットをＧＡ／ＳＡにより適合させ、発現プロフ
ァイルとの適応度の高いものをＭ個残して、トリプレッ
ト・プール４００に格納する（Ｓ２１０）。これによ
り、２Ｎ×Ｍ個のトリプレットがトリプレット・プール
４００に格納されることになる。この２Ｎ×Ｍ個のトリ
プレットから、適応度が高いものから所定のしきい値以
下の値を有するＰ個を選択してトリプレット・プール４
００に残す（Ｓ２１２）。これで、トリプレット・プー
ル４００に格納されているトリプレットはＰ個となっ
た。ここまでの処理は、トリプレットに対する広域探索
である。

【００２２】（局所探索）トリプレット・プール４００
に格納されているＰ個のトリプレットの近傍を、例えば
ＳＡ（シミュレーテッド・アニーリング）により探索し
て、より適応度の高いトリプレットを探す（Ｓ２１
４）。これを図示したのが、図８である。図８におい
て、Ｐ個の各トリプレットにおけるトポロジー空間のＳ
Ａの対象となる近傍に対して、探索を行っていることを
示している。これにより、より適応度の高いトリプレッ
トを得ることができる。より適応度が高いトリプレット
が見つかった場合は、それもトリプレット・プール４０
０に格納する（Ｓ２１６）。これにより、Ｐ個以上のト
リプレットがトリプレット・プール４００に格納される
ことになる。なお、このとき、より適応度が高いトリプ
レットが、あるトリプレットの近傍に発見された場合、
そのトリプレットに代えて、より適応度が高いトリプレ
ットをトリプレット・プール４００に格納することもで
きる。この場合、トリプレット・プール４００に格納さ
れているトリプレットの個数は、Ｐ個のままである。こ
こで、ＳＡを用いて、近傍の探索を行ったが、他の手
法、例えば山登り法（可能性のある探索点のなかで最も
有望な点を選んで探索を進めていく方法）等により行う
ことが可能である。

【００２３】このようにして、トリプレット・プール４
００に格納されているトリプレットの適応度を高めてい
るので、サンプリング段階を終了した時点で、トリプレ
ット・プール４００に格納されているトリプレット中
に、目的とするトリプレットが入っていることが、高い
確率であり得る。しかし、この段階のトリプレットの数
は、まだ、これに対して生物学的な実験を行うのに、十
分には小さくない。

【００２４】（絞り込み過程）次に、トリプレット・プ
ール４００に格納されているトリップレットから、より
少数の候補トリプレットを選別する処理の段階である。
この処理を、突然変異分析を用いて行っている図９およ
び図１０を例として説明する。図９において、トリプレ
ット・プール４００から各トリプレットに対して、ミュ
ータント（突然変異）・トリプレットを作成する（Ｓ４
０２）。この場合のミュータント・トリプレットは、あ
る遺伝子を除去（ノックアウト）して、その遺伝子から
の結合を全て取り去ることで作成している。このとき、
トリプレットにおけるその他の結合や、パラメータは変
化させない。これを遺伝子（ネットワーク構造のノー
ド）に対して行う。したがって、トリプレットのトポロ
ジーにおけるノードのうち、ｘ個の遺伝子のそれぞれの
ノードをノックアウトすると、ミュータント・トリプレ
ットをそれぞれ格納したミュータント・プール４５２〜
４５６はｘ個できる。そして、与えられたそれぞれのタ
ーゲット・データである突然変異プロファイルとの間の
適応度を各ミュータント・プールごとに評価する（Ｓ４
０４）。そして、各ミュータント・プールごとに評価し
た適応度を統合して、最終的に、ある値以上の適応度を
有する候補群を選出して、候補トリプレット・プール５
００に格納する（Ｓ４０６）。この場合、適応度によ
り、得られた候補トリプレットをソートしておいて、そ
の適応度を参照して、高い順に生物学的な実験を行って
もよい。ここで用いる突然変異としては、上述のノック
アウトの他に、例えば、ヘテロザイゴート（Heterozygo
te)，過剰発現(over expression)、温度感受性突然変異
(temperature sensitive mutation)等があり、いずれも
使用可能である。例えば、へテロザイゴートは、相同染
色体の特定の遺伝子座においての一方の対立遺伝子をつ
ぶしてしまうことにより、全体としてその遺伝子の発現
量が例えば半分などに減少してしまう突然変異である。
また、過剰発現は、ある遺伝子について過剰に発現して
しまう突然変異である。温度感受性突然変異は、温度感
受性のある物質を用いて、計測中のある時点で遺伝子の
発現を止めてしまう、又は、過剰に発現させる突然変異
である。

【００２５】（適応度の統合の計算例）この適応度の統
合の計算例について、図１０を用いて説明する。図１０
では、全ての突然変異ごと（突然変異ではない野生型も
含む）の各トポロジーごとの適応度を合計したものを、
そのトポロジーの適応度として取り扱う例を示してい
る。この計算例では単純な合計をしているが、必要に応
じて、突然変異ごとに重み付けて加重加算等を行っても
よい。なお、突然変異の発現プロファイルは、対象の全
ての突然変異遺伝子ネットワークに対して得られている
とは限らない。この場合は、得られている突然変異の発
現プロファイルのみを用いて、適応度の算出および統合
を行う。このようにして、最終的に、候補トリプレット
を得ることができる。

【００２６】（実施例）具体的な遺伝子のトリプレット
を例に、上述の処理を詳しく説明する。ここでは、図３
（ａ）に示したような時系列データである発現プロファ
イルを実現する遺伝子のネットワークを例とする。ま
ず、トポロジーの結合関係を示す結合行列をもとに染色
体にコーディングを行う。ここでは、図３（ａ）に示し
た発現物質が３つであるので、ノードが３つあるネット
ワークのトポロジーとしては、例えば、以下のような結
合行列Ｃで表される。これは、図３（ｂ）に示されてい
るものである。

【数２】

【００２７】このトポロジーの結合は４カ所あり、各々
に対する重み付けをする必要があるので、４つの遺伝子
型（genotype）を持つ染色体を生成する。さらに各ノー
ドのしきい値も最適化する必要があるので、結局、全体
として４（結合数）＋３（ノードのしきい値数）で７の
遺伝子型を持つ染色体を生成することになる。この場合
の染色体の荷重およびしきい値の取りうる範囲（実数
値）は、図１１（ａ）のようになる。このような図３
（ａ）に示す発現プロファイルを実現可能な複数のトポ
ロジーを発生して、トポロジー・プール３００（図４）
に格納する。この発生した各トポロジーに対して、ＧＡ
を用いてパラメータ（荷重およびしきい値）の最適化を
行う。まず、トポロジーに対して決められた範囲のパラ
メータをＮ組ランダムに発生させ、染色体の初期集団と
する。図１１（ｂ）に、上述の結合行列で定められたト
ポロジーに対する初期集団を示す。このように、トポロ
ジーに対して、リストの形で染色体が定まる。この初期
集団に対してＧＡを適用することにより、最適化を行う
のである。この最適化の指標として、各パラメータの適
応度を用いる。

【００２８】初期集団の適応度は、染色体のパラメータ
をもとに、各個体で発現プロファイルを生成し、ターゲ
ット・プロファイルとの間で、例えば最小自乗誤差ＴＳ
Ｓ（Total Sum Square error)の処理を行うことにより
求める。この最小自乗誤差ＴＳＳは、以下の計算を行
う。図３（ａ）に示すように、ターゲット・プロファイ
ルが濃度の時系列データとして与えられている。このた
め、各発現物質Ａ，Ｂ，Ｃそれぞれのある時間ｔにおけ
る濃度はＬ_TA（ｔ），Ｌ_TB（ｔ），Ｌ_TC（ｔ）と表さ
れ、生成した染色体の各パラメータをもとに、上述の式
（１）により計算される各発現物質の濃度は、Ｌ
_EA（ｔ），Ｌ_EB（ｔ），Ｌ_EC（ｔ）と表すことができ
る。これらの間の最小自乗誤差ＴＳＳは次のようにあら
わすことができる。

【数３】 このようにして求めた最小自乗誤差ＴＳＳを各発生した
ネットワーク（染色体）の適応度とする。

【００２９】さて、このように求められる適応度を利用
することにより、ＧＡを用いてトポロジー毎にパラメー
タ（結合荷重としきい値）の最適化を行うことができ
る。この処理について、図１２，図１３を用いて、ＧＡ
のクロスオーバと突然変異を用いて最適化を行っている
例で説明する。

【００３０】図１２（ａ）で示されているように、トポ
ロジーのパラメータは、結合荷重としきい値のリストで
表されている。まず、親となる個体対の決定を、例えば
トーナメント戦略を用いて行う。これは、トーナメント
・サイズの数の個体（例えば２）をランダムに選択す
る。この選択した個体の中から一番適応度が高い個体を
一方の親とする。これをもう一度行うことにより、親と
なる個体対が決定される。この個体対に対して、図１２
（ｂ）に示すように、交叉（crossover）を行って子孫
のパラメータを作成し、さらに、突然変異を適用して、
生成したパラメータで適応度を計算する。このような処
理を、例えば、定められた数の子孫が作成されるまで行
うと、第２世代の個体集団が得られる。

【００３１】交叉および突然変異の処理を、図１３を用
いて具体的に説明する。図１３（ａ）において、リスト
で表された親１および親２のパラメータの矢印で示され
た部分を入れ替えて、子１および子２が作成される。こ
れが交叉の処理である。図１３（ｂ）は、子１の丸を付
けられたパラメータの１つ（結合荷重０．４）が、別の
値（結合荷重０．６）に変化したことが示されている。
この変化が突然変異で、例えば、固定された確率でラン
ダムにおきるように処理される。この交叉と突然変異の
処理を、例えば、決められた世代数分繰り返すことで、
あるトポロジーに対するパラメータを最適化することが
できる。なお、全世代で、エリート保存戦略を適用して
もよい。エリート保存戦略は、集団中で最も適応度の高
い個体を、交叉や突然変異により変化させず、そのまま
次世代に残す手法である。

【００３２】以上の処理をすべてのトポロジーに対して
終了したときに、その中で、所定の適応度Ｔ⁺（この例
の場合、最小自乗誤差）以下の適応度を持つトポロジー
のサブグループＤ⁺を決定する。このサブグループＤ
⁺は、次の式により定義される。

【数４】 抽出されたトポロジー・サブグループのトリプレット
（トポロジー，パラメータ，適応度）を、トリプレット
・プール４００へ格納する。なお、この実施例ではトポ
ロジー・プール３００に格納されているすべてのトポロ
ジーに対してＧＡにより最適なパラメータを作成してい
るが、トポロジー・プール３００に格納されたトポロジ
ーが莫大な数である場合には、定められた数のトポロジ
ーをランダムに選択して、これに対して最適なパラメー
タをＧＡにより生成してもよい。ここまでが、トリプレ
ットの広域探索である。

【００３３】このトリプレット・プールに格納されてい
るトリプレットに対して、さらに別のＧＡ手法を適用し
て最適化したものを追加する。ここで適応するＧＡは２
次元ＧＡである。この２次元ＧＡについて図１４を用い
て説明する。図１４（ａ）は、結合行列で表した例で、
それを図式化したものを図１４（ｂ）に示す。図１４
（ａ）に示すように、親１および親２の結合行列の十字
形で示される部分（即ちＢを中心とする結合関係：図１
４（ｂ）参照）を交叉させて、子１および子２を生成す
る。このように、図１３（ａ）で説明したリスト（パラ
メータ）上の交叉とは異なり、ネットワークの物理的な
構造自体に着目して交叉を行っている。このような処理
は、２次元構造の最適化に有用である。２次元ＧＡを適
用して得られたトポロジーに対してパラメータを、ＧＡ
を前述と同様に用いることで最適化を行って求める。こ
れを例えば、決められた世代数実行する。そして、上述
と同じように、適応度Ｔ⁺以下であるトリプレットのみ
を、さらにトリプレット・プール４００に格納する。こ
の追加して行うＧＡでもエリート保存戦略を適用するこ
とができる。

【００３４】トリプレット・プール４００に格納されて
いるトリプレットを対象としてＳＡ（シミュレーティド
・アニーリング）を実行して、各トポロジーの近傍トポ
ロジー空間を探索し、得られたトポロジーに対して、前
と同様にＧＡにより最適化されたパラメータを求める。
この近傍の探索を行うのは、ＧＡが局所的な探索能力に
劣っているためである。ここで用いているＳＡの基本ア
ルゴリズムは、疑似コードを用いると、以下のようにな
る。これは、各トポロジーに対して適用される部分であ
る。 上述のACCEPT(E, E', T)中のＥは、ＳＡではエネルギー
と呼ばれる状態を規定する変数で、この場合は、各ネッ
トワークの適応度を使用している。Ｅ’は、「次の状態
生成（）」処理で生成されたトリプレットの適応度であ
る。また、Ｔは、ＳＡのアルゴリズムでは温度パラメー
タで、この場合は、この処理のみで用いられる状態を規
定する変数である。変数Ｔは、処理が始まるときに、初
期値Ｔ_{in it}に設定される。また、停止条件や平衡条件
は、例えば、処理回数を用いることができる。「ACCEPT
(E, E', T)」， 「REDUCE(T)」は、次のように定義され
る。

【数５】 上述の疑似コードで表される処理を説明する。処理対象
のトポロジーが選択されると、変数ＴとＥとが初期化さ
れる。そして、回数で規定されるループの内側の繰り返
し処理で、まず、「次の状態発生（）」という処理を行
う。この「次の状態発生（）」は、以下のような処理を
行う。処理対象のトポロジーを基に、結合関係を変化さ
せたトポロジーを発生させる。この変化は、図１５
（ａ）に示した遷移図に基づき、処理対象のトポロジー
を、確率Ｐにより変化させる。この遷移確率は、上述の
変数Ｔを用いて、次のように定義される。

【数６】 ここで、Ｔは上述の変数、しきい値Thresholdは、ここ
では例えば、０．５に設定されている。しきい値を導入
した理由は、遷移確率をあまり高くすると、対象とした
トポロジーに対して多数の変化したトポロジーを探索す
ることになり、近傍探索とはならないからである。Ｔ
は、外側のループの「REDUCE(T)」で変化するので、内
側のループ内では遷移確率は変化しない。外側のループ
のステップにより、遷移確率がどのように変化するのか
を示したのが、図１５（ｂ）である。この場合、γは例
えば、０．８である。トポロジーの変化を図１５（ｃ）
に示したトポロジーを例に説明する。これは、図３
（ｂ）のトポロジーと同じものである。なお、遷移確率
を０．５とした場合で説明する。

【００３５】図１５（ｃ）において、ＡからＡの結合
は、活性化で結合しているので、図１５（ａ）の状態遷
移図では、活性６０６の状態にある。この状態からは、
未結合６０４への状態遷移のみが可能である。この状態
遷移が起こる確率は０．５である。次に、ＡからＢへの
結合を見るとこれも活性として結合しているので、同様
に、０．５の確率で、未結合の状態に変化する。Ｂから
Ｂへの結合は、未結合なので、抑制と活性とに０．２５
づつ（合計で０．５）の確率で遷移する。以下同様に、
すべての結合関係について、図１５（ａ）に示した状態
遷移図と、図１５（ｂ）の確率に基づいて、変化させ
る。このようにして、対象のトポロジーを構成する各結
合関係を上述で説明した確率で、各遷移図に従って変化
させた結果得られたトポロジーに対して、次に、前述の
ＧＡを用いて、最適なパラメータを決定し、そのパラメ
ータに対する適応度を計算する。これで新しいトリプレ
ットが作成される。これで「次の状態発生（）」は完了
し、次の処理である「ACCEPT(E, E', T)」を実行する。

【００３６】この「ACCEPT(E, E', T)」では、処理対象
であるトポロジーを有するトリプレットの適応度（Ｅ）
と、発生したトポロジーの適応度（Ｅ’）を比較し、発
生させたトポロジーの適応度が低い（より適合してい
る）と、１（真）を返す。また、発生した適応度が高く
ても（適合していなくても）、ある確率で１（真）を返
している。この場合、Ｔが高い場合（あまり処理回数が
多くないとき）に１（真）を返す確率を多くしている。
「ACCEPT(E, E', T)」は１（真）を返すと、「更
新（）」処理が行われて、「次の状態発生（）」で発生
したトリプレットに更新する。このようにして、内側の
ループを平衡状態が得られるまで（この場合はある回数
の処理が行われるまで）続ける。所定回数の処理が行わ
れると、外側の処理（「REDUCE(T)」）が行われて、新
しいＴが定められ、また、このＴを用いて、内側の処理
が行われる。外側の処理も、この場合、例えばある所定
回数行われる。そして、最終的にたどり着いたトリプレ
ットの適応度がＴ⁺であると、トリプレット・プールに
入れる。

【００３７】上述のようにこの段階終了時では、トリプ
レット・プール４００には、各トポロジーの近傍を探索
して、より良い適応度のトリプレットと置き換えられた
ものが格納されている。このトリプレットから、最終的
な候補トリプレットを選択する処理を以下で説明する。
これは、種々の突然変異の発現プロファイルを用いて、
これに対する適応度を計算することにより行う。この処
理を、図１６のフローチャートを用いて説明するこれに
は、まず、図３（ｃ）、図３（ｄ）に示されているよう
な、突然変異によるプロファイルを用意する。図３
（ｃ）、図３（ｄ）は、遺伝子Ｂがない突然変異（ミュ
ータント）のプロファイル、遺伝子Ｃがない突然変異
（ミュータント）のプロファイルである。そして、それ
ぞれに対応したミュータント・プールを用意し、上述の
トリプレット・プール内の全トリプレットを各ミュータ
ント・プールに格納する（Ｓ８０２）。各ミュータント
・プールに格納されたトリプレットに対して、トリプレ
ットのトポロジーをその突然変異のトポロジーとし、パ
ラメータはそのまま用いて、対応する突然変異のプロフ
ァイルに対する適応度を計算する（Ｓ８０４）。各ミュ
ータント・プールにおけるトリプレットのうち、計算し
た適応度が各突然変異ごとの定めた所定のしきい値以下
であるトリプレットを選択する（Ｓ８０６）。選択した
各ミュータント・プール内のトリプレットで、全部のミ
ュータント・プールに共通するトリプレットであり、か
つ、所定のしきい値以下の適応度を有するものに対応す
るトリプレットを、候補トリプレット・プールへ格納す
る（Ｓ８０８）。これで、適切な候補トリプレットを確
定することができる。

【００３８】（他の適用）上述の説明では、遺伝子ネッ
トワークの例で説明したが、例えば遺伝子・代謝回路ネ
ットワークにも適用することが可能である。図１７は、
遺伝子・代謝回路のネットワークについて示している。
図１７（ａ）に示すように、遺伝子と、代謝における酵
素、タンパク質をノードとして、活性、抑制、媒介等の
関係をリンクで表す。そして、これを用いたネットワー
ク例を図１７（ｂ）に示す。この図１７（ｂ）に示した
ようなネットワーク構造に対して、上述で説明した処理
を適用する。上述においては、対象ネットワークとし
て、遺伝子ネットワーク、代謝回路ネットワークを示し
たが、そのほかにも神経回路ネットワーク等に対しても
適用することができる。この場合のデータ（ターゲット
・プロファイル）としては、ニューロン活動電位等を用
いる必要がある。このように、本発明の手法は色々なも
のに対して用いることができる。適用できるターゲット
・プロファイルの色々な例を図１８に示す。図１８
（ａ）は、ターゲット・プロファイルが０次元空間デー
タである場合を示している。このデータの例としては、
時間や周波数に対して変化するようなデータである。例
えば、濃度、活性度、電位等がある。図１８（ｂ）は、
１次元空間データを示している。データは、空間として
ｘに、および、時間又は空間に依存する量である。この
量としては、濃度、活性度、電位等がある。図１８
（ｃ）および（ｄ）として示したターゲット・プロファ
イルは、２次元空間データおよび３次元空間データであ
る。これは、２次元（ｘ，ｙ）および３次元（ｘ，ｙ，
ｚ）の空間、および、時間又は周波数に依存する量のデ
ータである。この量としては、濃度、活性度、電位等が
ある。図１２に示すように、色々なデータ形式で、ター
ゲット・プロファイルを与えることができる。

【００３９】このように、ターゲット・プロファイルが
存在し、そのターゲット・プロファイルを生じるネット
ワーク構成を推定するような場合に、本発明は有効であ
る。本発明は、スタンド・アローンのコンピュータ・シ
ステムばかりではなく、複数のコンピュータから構成さ
れる例えばクライアント・サーバ・システム等に適用し
てもよい。本発明に関するプログラムを格納した記憶媒
体から、プログラムをシステムで読み出して実行するこ
とにより、本発明の構成を実現することができる。この
記録媒体には、ＤＶＤ、ＭＤ、ＭＯ、フロッピー・ディ
スク、ＣＤ−ＲＯＭ、磁気テープ、ＲＯＭカセット等が
ある。

【発明の効果】上記の説明のように、本発明は、データ
としてターゲット・プロファイルが存在し、そのターゲ
ット・プロファイルを生じるネットワーク構成を推定す
るような場合に、極めて有効である。特に、遺伝子ネッ
トワークや代謝ネットワークのように、分子生物学にお
けるモデルの推定に対して有効に適用することができ
る。

【図面の簡単な説明】

【図１】分子生物学におけるコンピュータ利用を説明す
る図である。

【図２】遺伝子ネットワークの例を示す図である。

【図３】発現プロファイルおよびそれに対応するネット
ワークを示す図である。

【図４】実施形態における処理例の概略フローチャート
である。

【図５】ネットワークのファンアウト制限を説明する図
である。

【図６】選択処理のフローチャートである。

【図７】選択処理の手法を説明するための図である。

【図８】ＳＡによる再探索を説明する図である。

【図９】候補の選別処理を説明する図である。

【図１０】適応度の総合処理の例を説明する図である。

【図１１】ＳＡを行う場合の表記例および初期集団を示
す図である。

【図１２】ＧＡの適用例を示す図である。

【図１３】具体的にＧＡ処理（交叉と突然変異）を行っ
ている例を示す図である。

【図１４】２次元ＧＡを適用した例を示す図である。

【図１５】ＳＡによる近傍探索を説明するための例であ
る。

【図１６】突然変異による候補トリプレット選択の処理
を示すフローチャートである。

【図１７】遺伝子・代謝ネットワークの例を示す図であ
る。

【図１８】ターゲット・プロファイルのデータ例を示す
図である。

───────────────────────────────────────────────────── フロントページの続き (73)特許権者 599079344 濱橋 秀互 東京都江東区木場６−８−15 ハウス木 場３Ｃ (72)発明者 北野 宏明 埼玉県川越市西小仙波町２−18−３ (72)発明者 諸橋 峰雄 東京都北区西ケ原１−27−３−1010 (72)発明者 京田 耕司 神奈川県横浜市青葉区奈良町2423−208 (72)発明者 濱橋 秀互 神奈川県川崎市中原区木月1566佐藤ビル 303号室 (56)参考文献 特表2002−518021（ＪＰ，Ａ) 特許2881711（ＪＰ，Ｂ２) 国際公開98／30722（ＷＯ，Ａ１) 富永大介、岡本正宏，「逆問題（ｉｎ ｖｅｒｓｅ ｐｒｏｂｌｅｍ）解決のた めの遺伝的アルゴリズムを用いた多次元 非線形数値最適化手法の開発」，化学工 学論文集，日本，社団法人化学工学会, 1999年 ３月10日，Ｖｏｌ．25，Ｎｏ. ２，ｐｐ．220−225，ＪＳＴ資料番号： Ｓ0110Ｂ 篠原歩、竹田正幸、野田清志，「重み 付きネットワークを用いた遺伝子発現ネ ットワークの発見」，1998年度科学研究 費補助金研究成果報告書，日本，1999年 ３月，ｐｐ．290−293，ＪＳＴ資料番 号：Ｎ19993074 矢田哲士・他，「隠れマルコフモデル と遺伝的アルゴリズムによるＤＮＡ配列 のシグナルパターン抽出」，情報処理学 会論文誌，日本，社団法人情報処理学 会・発行，1996年 ６月15日，Ｖｏｌ. 37，Ｎｏ．６，ｐｐ．1117−1129 矢田哲士・他，「ＤＮＡ配列の複合モ チーフを表現する隠れマルコフモデルの 生成」，情報処理学会論文誌，日本，社 団法人情報処理学会・発行，1999年 ２ 月15日，Ｖｏｌ．40，Ｎｏ．２，ｐｐ. 750−767 阿久津達也，「遺伝子ネットワークの 推定アルゴリズム」，数理科学，日本, 株式会社サイエンス社・発行，1999年 ５月28日，Ｖｏｌ．37，Ｎｏ．６，ｐ ｐ．40−46 北野宏明，「システム・バイオロジー の試み」，数理科学，日本，株式会社サ イエンス社・発行，1999年 ５月28日, Ｖｏｌ．37，Ｎｏ．６，ｐｐ．54−60 北野宏明・他，「遺伝的アルゴリズ ム」，日本，産業図書株式会社，1993年 ６月 ３日，初版，ｐｐ．61−88，Ｉ ＳＢＮ：４−7828−5136−７ 阿久津達也・他，「遺伝子発現量の時 系列データからの遺伝子ネットワーク推 定」，1999年度人工知能学会全国大会 （第13回）論文集，日本，社団法人人工 知能学会，1999年 ６月15日，ｐｐ. 328−329，ＪＳＴ資料番号：Ｘ0580Ａ 阿久津達也・他，「遺伝子制御ネット ワーク解明のための学習理論」，人工知 能学会研究会資料（ＳＩＧ−ＦＡＩ− 9703−６），日本，社団法人人工知能学 会，1998年 ３月26日，ｐｐ．37−42 北野宏明・他，「遺伝的アルゴリズム ２」，日本，産業図書株式会社，1995年 ５月 ８日，初版，ｐｐ．71−87，Ｉ ＳＢＮ：４−7828−5141−３ Ｔｏｍｉｎａｇａ，Ｄ．，ｅｔ．ａ ｌ，”Ｄｅｓｉｇｎ ｏｆ Ｃａｎｏｎ ｉｃａｌ Ｍｏｄｅｌ Ｄｅｓｃｒｉｂ ｉｎｇ Ｃｏｍｐｌｅｘ Ｎｏｎｌｉｎ ｅａｒ Ｄｙｎａｍｉｃｓ”，Ｐｒｏ ｃ．ＩＦＡＣ Ｉｎｔ．Ｃｏｎｆ．Ｃｏ ｍｐｕｔｅｒ Ａｐｐｌｉｃａｔｉｏｎ ｓ ｉｎ Ｂｉｏｔｅｃｈｎｏｌｏｇ ｙ，1998年，ｐｐ．85−90，ＪＳＴ資料 番号：Ｋ19990184 久原哲・他，「ＤＮＡチップの情報科 学的取り扱い」，数理科学，日本，株式 会社サイエンス社，1999年 ６月 １ 日，Ｖｏｌ．37，Ｎｏ．６，ｐｐ．33− 39 (58)調査した分野(Int.Cl.⁷，ＤＢ名) G06N 1/00 - 7/08 G06F 19/00 C12N 15/00 ＪＳＴファイル（ＪＯＩＳ) ＣＳＤＢ（日本国特許庁) ＩＮＳＰＥＣ（ＤＩＡＬＯＧ)

Claims (6)

(57)【特許請求の範囲】
1. 【請求項１】 あるデータを生成する要素を少なくとも
   １つ用いて、互いに関連している要素間の関係をネット
   ワークとして記述でき、前記ネットワークにより生成さ
   れたデータが与えられたときにおける、前記与えられた
   データから、前記データを再現できるネットワークをコ
   ンピュータ・システムを用いて推定する方法であって、前記ネットワークは、 ネットワーク構造と、 それに対応するパラメータ・セットと、 ネットワーク構造およびパラメータ・セットから生成さ
   れるデータと前記与えられたデータとの適応度とからな
   るトリプレットとして表現され、 前記与えられたデータを実現することが想定される、部
   分的に既知のネットワーク構造を基にネットワーク構造
   を生成し、それに対応するパラメータ・セットおよび適
   応度を作成し、前記適応度を利用して最適化した複数の
   候補トリプレットを第１記憶手段に格納するステップ
   と、 前記与えられたデータとは異なる、ネットワーク
   構造の一部が変異したネットワーク構造から生成できる
   他のデータを用いて、前記第１記憶手段に格納された候
   補トリプレットから適切な候補トリプレットを絞り込
   み、第２記憶手段に格納するステップとを備えることを
   特徴とするネットワーク推定方法。
2. 【請求項２】 請求項１に記載のネットワーク推定方法
   において、前記候補トリプレットを第１記憶手段に格納するステッ
   プの前記適応度を利用した最適化は、作成したネットワ
   ーク構造からＮ個を選択して、各々に最適化したパラメ
   ータ・セットをＭ個適合してＮ×Ｍ個のトリプレットを
   作成し、さらに、該Ｎ個のネットワーク構造から新たな
   Ｎ個のネットワーク構造を生成して、各々に最適化した
   パラメータ・セットをＭ個適合してＮ×Ｍ個のトリプレ
   ットを作成し、作成した２Ｎ×Ｍ個のトリプレットから
   適応度の高いＰ個のトリプレットを選択することを特徴
   とするネットワーク推定方法。
3. 【請求項３】 請求項２に記載のネットワーク推定方法
   において、前記候補トリプレットを第１記憶手段に格納するステッ
   プの前記適応度を利用した最適化は、さらに、前記Ｐ個
   のトリプレットの近傍を探索して、より適合度の高いト
   リプレットを見つけて置き換えることを特徴とするネッ
   トワーク推定方法。
4. 【請求項４】 請求項１又は２に記載のネットワーク推
   定方法において、 前記候補トリプレットを第１記憶手段に格納するステッ
   プの前記適応度を利用した最適化は、遺伝的アルゴリズ
   ム（ＧＡ）、シミュレーテッド・アニーリング（Ｓ
   Ａ）、及び／又は、山登り法を用いることを特徴とする
   ネットワーク推定方法。
5. 【請求項５】 あるデータを生成する要素を少なくとも
   １つ用いて、互いに関連している要素間の関係をネット
   ワークとして記述でき、前記ネットワークにより生成さ
   れたデータが与えられたときにおける、前記与えられた
   データから、前記データを再現できるネットワークを推
   定するネットワーク推定装置であって、ネットワークを、 ネットワーク構造とそれに対応するパ
   ラメータ・セットと、ネットワーク構造とパラメータ・
   セットから生成されるデータと前記与えられたデータと
   の適応度とからなる候補トリプレットとして表現して記
   憶する第１の記憶手段と、 最終的な候補としてのトリプレットを格納する第２の記
   憶手段と、 前記与えられたデータを実現することが想定される、部
   分的に既知のネットワーク構造を基にネットワーク構造
   を生成し、それに対応するパラメータ・セットおよび適
   応度を作成し、前記適応度を利用して最適化した複数の
   候補トリプレットを前記第１の記憶手段に格納する手段
   と、 前記与えられたデータとは異なる、ネットワーク
   の一部が変異したネットワークから生成できる他のデー
   タを用いて、前記第１の記憶手段に格納されているトリ
   プレットから適切な候補トリプレットを絞り込み、前記
   第２の記憶手段に格納する手段とを備えることを特徴と
   するネットワーク推定装置。
6. 【請求項６】 請求項１〜４のいずれかに記載のネット
   ワーク推定方法をコンピュータ・システムに実現させる
   プログラムを格納した記録媒体。