集中講義(九州大学数理学研究院) バイオ構造データに対する数理モデルと アルゴリズム(4) ブーリアンネットワーク 阿久津 達也 京都大学 化学研究所 バイオインフォマティクスセンター
内容 ブーリアンネットワーク 定常状態の高速検出アルゴリズム 確率ブーリアンネットワーク 制御系列計算のアルゴリズム 代謝ネットワークのブーリアンモデルと頑健性解析
計算量 情報科学では、入力データのサイズ(n)に対して、計算時間がどのように変化するかを理論的に解明することが重要 O(n): かなり速い(文字列検索など) O(n log n): 結構速い(ソートなど) O(n2): まあまあ速い(アライメントなど) O(n3): ちょっと遅い(RNA二次構造予測など) O(n4): 結構遅い(Pseudo-knotつきRNA二次構造予測など) NP困難: すごく遅い (マルチプルアライメント、スレッディングなど) P=NP は理論計算機科学における最大の難問 P≠NPならば、NP困難問題に対する理論的に効率的なアルゴリズム(多項式時間アルゴリズム)は存在しない しかし、タンパク質配列などは n ≦ 1000 くらいなので、実用アルゴリズムを開発できる可能性はある
ブーリアンネットワーク
ブーリアンネットワーク 遺伝子ネットワークの数理モデル 頂点⇔遺伝子 制御規則 状態は(クロックに)同期して変化 状態: 1 (発現) と 0 (非発現)のみ 制御規則 ブール関数 (AND, OR, NOT …) もし、y が x を直接制御していれば、y から x に辺が存在 状態は(クロックに)同期して変化 基本的にデジタル回路と同じ [Kauffman, The Origin of Order, 1993]
ブーリアンネットワークの例 A B C ブーリアンネットワーク 状態遷移表 ’ = B ’ = A and C ’ = not A A time t t+1 B C 1 INPUT OUTPUT 状態変化の例: 111 ⇒ 110 ⇒ 100 ⇒ 000 ⇒ 001 ⇒ 001 ⇒ 001 ⇒ 。。。
ブーリアンネットワーク研究の意義 単純化しすぎとの批判 最も単純な非線形システムの一つ 本質的にデジタル回路と同じ でも、単純化しないと理論解析、推定、制御などが困難 シミュレーションなら複雑なモデルでも可能 定性的理解には有用(?) 最も単純な非線形システムの一つ ブーリアンネットワークで不可能なら、他の(多くの)非線形システムでも不可能 本質的にデジタル回路と同じ 計算機科学(特にブール関数、アルゴリズム理論)における成果や技法が利用可能
定常状態の高速計算アルゴリズム
アトラクター (1) 定常状態 細胞の種類に対応? 例 011 ⇒ 101 ⇒ 010 ⇒ 101 ⇒ 010 ⇒… アトラクター (1) 定常状態 細胞の種類に対応? 例 011 ⇒ 101 ⇒ 010 ⇒ 101 ⇒ 010 ⇒… 111 ⇒ 110 ⇒ 100 ⇒ 000 ⇒ 001 ⇒ 001 ⇒ 001 ⇒ … A ’ B’ C’ time t t+1 B C 1 INPUT OUTPUT
アトラクター (2) A ’ B’ C’ B C 1 time INPUT OUTPUT 000 010 001 101 100 110 アトラクター (2) A ’ B’ C’ time t t+1 B C 1 INPUT OUTPUT 000 010 001 101 100 110 011 111
点アトラクター (もしくは、 ) アトラクターの生物学的解釈 異なるアトラクター ⇔ 異なる種類の細胞 点アトラクター 異なるアトラクター ⇔ 異なる種類の細胞 点アトラクター 周期が1のアトラクター (静的)定常状態に対応 以下を満たす状態 (もしくは、 )
アトラクターの分布に関する研究 平均入次数を定数(K)とし、BNをランダムに生成した場合のアトラクターのサイズや個数 多くの研究があるが、詳しいことはわかっていない 点アトラクター(非周期的)であれば平均、1個であることが簡単にわかる 計算機実験により、n頂点の場合、 個との予想[Kauffman, The Origin of Order, 1993]があるが、正しくなさそう? [Samuelsson & Troein, PRL, 2003], [Drossel et al., PRL, 2005]
定理:各頂点に対してランダムにn変数ブール関数を選んでブーリアンネットワークを構成した場合、周期1のアトラクターの個数は平均的に1個 証明 時刻 t=0 においてすべての頂点の状態が0であったとする(すべての i について vi(0)=0 であったとする) 任意の頂点 vi について、時刻 t=1 に vi の状態が 0 となる(vi(1)=0となる)確率は 1/2 よって、すべての i について vi(1)=0 となる確率は 1/2n(つまり、000…0 がアトラクターとなる確率が1/2n) 時刻 t=0 におけるすべての状態(000..00から111…1まで)についての和をとることにより、周期1のアトラクターの個数の期待値は 2n×(1/2n) = 1 参照:[Harvey et al., ECAL, 1997; Mochizuki, JTB, 2005]
アトラクターが平均1個の理由 t=0 t=1 A B C B A C Prob( A(1) =1 ) = 0.5 ランダムに ブール関数を選ぶ A B C t=0 t=1 確率 0.5
アトラクターの高速計算 n 頂点の場合、2n 個の状態があるので、 程度の時間かければ計算可能 点アトラクターの有無の判定はNP完全 ヒューリスティックな高速アルゴリズムはあるが、計算時間に理論的保証が無い [Irons, Pysica D, 2006], [Devloo et al., Bull. Math. Biol. 2003], … 点アトラクターの有無の判定はNP完全 [Akutsu et al., GIW 1998] 平均的に高速なアルゴリズムを開発 [Shuqin et al., EURASIP JBSB 2007]
再帰型アトラクター検出アルゴリズム(1) 0-1割り当てを1個目の頂点からn個目の頂点までへと順次試していき、矛盾が検出されたら後戻りする [Zhang et al., EURASIP JBSB 2007]
再帰型アトラクター検出アルゴリズム(2) 1 1 1 Output
再帰型アトラクター検出アルゴリズム(3) K 2 3 4 5 6 1.19n 1.27n 1.34n 1.41n 1.45n 0-1割り当てを最初の頂点から順々に探索。点アトラクターにならないとわかった時点でバックトラック 0 00 -> X 01 010 -> X 011 -> X 10 これは大幅な高速化 (測定値ではなく、理論値) 230≒109 1.230≒240 2100=1030 1.2100≒80,000,000 小さい周期のアトラクター検出にも対応可能 K 2 3 4 5 6 平均計算量 1.19n 1.27n 1.34n 1.41n 1.45n
v1 vm-1 vm vm+1 点アトラクターの高速計算:計算量解析 t=0 t=1 K 全部でn個の頂点のうちm個の頂点まで0-1割り当てを行った時、vi(0)≠vi(1) が検出できる確率は よって、m個すべての頂点がアトラクターの条件を満たし、m+1個目の頂点の0-1割り当てがチェックされる確率は よって、m個の頂点に対する0-1割り当て2m通りのうち、次のステップに再帰呼び出しが進む回数の期待値は よって、平均計算量は上の式において m を1からnまで動かした場合の最大値のオーダーとなる v1 vm-1 vm vm+1 t=0 t=1 K
最悪時の時間計算量 入次数 K 以下のBNの点アトラクター検出 (K+1)-SAT に変換可能 点アトラクター検出はK=2でもNP困難 K=2 の場合、O(1.322n) 時間 さらに工夫をすれば、O((1.322-δ)n) 時間にすることも可能 点アトラクター検出はK=2でもNP困難 ブール関数を変数とその否定のAND/ORに限った場合にはKに関わらず O(1.587n) 時間で検出可 [Melkman, Tamura & Akutsu, 2010]
アトラクター検出問題からSATへの変換 入次数 K 以下のBNの点アトラクター検出 (K+1)-SAT vj vk vi
ブーリアンネットワークの制御
システム生物学 北野宏明氏が提唱 e.g., [Kitano, Nature, 420:206-210, 2002] 生命システムの制御理論の構築(主要目標の一つ) 線形制御理論は確立しているが、生命システムは非線形 生命システムに対する制御理論を構築することが必要 創薬やがん治療などにつながる可能性がある iPS細胞では4個の転写因子を導入することにより、多能性幹細胞を構成
先行研究 (Datta et al.) NO! [Akutsu et al., 2007] (ただし、P≠NPの仮定のもとで) 確率ブーリアンネットワーク (PBN)の制御 いくつかの頂点(遺伝子)を制御可能と仮定 目標状態、初期状態を与え、目標状態にできるだけ近づき、かつ、コストが低い制御系列を計算 古典的制御理論に基づく動的計画法アルゴリズムを提案 ⇒ 2n×2n サイズの行列を扱うことが必要 ⇒ 指数時間アルゴリズム 素朴な疑問 多項式時間で制御系列を見つけることは可能か? NO! [Akutsu et al., 2007] (ただし、P≠NPの仮定のもとで) [Datta et al., Machine Learning, 2003]
確率ブーリアンネットワーク (PBN: Probabilistic Boolean Network) 一つの頂点に対する制御規則が複数 どの制御規則が選ばれるかは(事前に定められた確率分布に従い)ランダムに決まる 利点:ノイズなどの影響が扱える。マルコフ過程に関する理論が適用できる 欠点:ほとんどの場合、O(2n) 時間以上の計算時間がかかるため、大規模なシステムを扱えない A B C A(t+1) = B(t) AND C(t) A(t+1) = B(t) OR (NOT C(t)) 確率0.6で 確率0.4で
ブーリアンネットワークの行列表現 状態遷移表 A ’ B’ C’ B C 1 とすると この行列を A とし、
PBNの行列表現 この行列の意味は以下のとおり
Dattaらのモデル コスト関数 確率ブーリアンネットワークに基づく 時刻 t における全体の状態: 制御規則: 2n×2n 行列 A A は m 個のブール変数により変化 コスト関数 C(x(t), u(t)): u(t) という制御を状態 x(t) に適用した際のコスト C(x(t)): 最終状態 v(t) のコスト
Dattaらの制御手法 入力: 初期状態 x(0), 制御規則 A(u(t)), 目的時刻 M , コスト関数 出力: コストの合計を最小化する制御系列 u(0), u(1), …, u(M-1) 動的計画法アルゴリズム
ブーリアンネットワーク制御モデルの定式化 入力 内部頂点: v1 ,…, vn , 制御頂点: u1 ,…, um 初期状態: v0 目標状態: vM ブーリアンネットワーク 出力 以下を満たす制御頂点の状態系列: u(0), u(1), …, u(M) v(0)=v0, v(M)=vM (時刻0で初期状態、時刻Mで目標状態) [Akutsu et al., J. Theo. Biol. 2007]
計算困難性に関する結果 制御系列を見つけるのは一般に計算困難(NP困難) ネットワークが以下のような単純な構造をしていても計算困難 ⇒ Dattaらの開発したような 非効率的アルゴリズムを使うのも仕方がない
BN制御の動的計画法アルゴリズム DattaらのアルゴリズムをBNに合わせ簡単化 DPテーブル アルゴリズム(再帰式) 状態 から目標状態に至る制御系列があれば 1、それ以外は0 アルゴリズム(再帰式)
動的計画法アルゴリズムの説明 D[1,1,1, 2] =1 D[0,0,0, 2] = 0 D[0,1,1, 3] = 1 u1=1, u2=1 DP計算 D[0,1,1, 3] = 1 ただし、DP表は 指数サイズ
ネットワークが木構造の場合 効率よく(多項式時間で)制御系列が計算可能 木構造: ループ無しの連結グラフ
木構造に対するアルゴリズム (1) 基本アイデア: 以下の S[vi,t,b] という表を、 t=0 から t=M という順に、かつ、木の葉から根へという順に、動的計画法で計算 S[vi,t,b]
木構造に対するアルゴリズム (2) S[vi,t,b] の計算: S[vi,t,b] の計算例 木構造に対するアルゴリズム (2) S[vi,t,b] の計算: S[vi,t,b] の計算例 S[v3,t+1,1]=1 if S[v1,t,1]=1 and S[v2,t,1]=1 S[v3,t+1,0]=0 if S[v1,t,0]=1 or S[v2,t,0]=1 Note: このアルゴリズムは Datta et al. の 動的計画法アルゴリズムと本質的に異なる
まとめ ブーリアンネットワーク 点アトラクター検出 確率ブーリアンネットワークの制御 ブーリアンネットワークの制御 頂点⇔遺伝子、制御規則⇔ブール関数 点アトラクター検出 再帰的アルゴリズム(O(2n)時間より高速) SATへの変換 確率ブーリアンネットワークの制御 動的計画法による指数時間アルゴリズム ブーリアンネットワークの制御 木構造に対する多項式時間アルゴリズム