集中講義（九州大学数理学研究院） バイオ構造データに対する数理モデルと アルゴリズム（４） ブーリアンネットワーク 阿久津　達也 京都大学　化学研究所 バイオインフォマティクスセンター

内容 ブーリアンネットワーク 定常状態の高速検出アルゴリズム 確率ブーリアンネットワーク 制御系列計算のアルゴリズム 代謝ネットワークのブーリアンモデルと頑健性解析

計算量 情報科学では、入力データのサイズ(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 ｔ ｔ＋１ B C 1 INPUT OUTPUT 状態変化の例： １１１　⇒　１１０　⇒　１００　⇒　０００　⇒　００１　⇒　００１　⇒　００１　⇒　。。。

ブーリアンネットワーク研究の意義 単純化しすぎとの批判 最も単純な非線形システムの一つ 本質的にデジタル回路と同じ でも、単純化しないと理論解析、推定、制御などが困難 シミュレーションなら複雑なモデルでも可能 定性的理解には有用（？） 最も単純な非線形システムの一つ ブーリアンネットワークで不可能なら、他の（多くの）非線形システムでも不可能 本質的にデジタル回路と同じ 計算機科学（特にブール関数、アルゴリズム理論）における成果や技法が利用可能

定常状態の高速計算アルゴリズム

アトラクター (1) 定常状態 細胞の種類に対応？ 例 011 ⇒ 101 ⇒ 010 ⇒ 101 ⇒ 010 ⇒… アトラクター　(1) 定常状態 細胞の種類に対応？ 例 011 ⇒ 101 ⇒ 010 ⇒ 101 ⇒ 010 ⇒… 111 ⇒ 110 ⇒ 100 ⇒ 000 ⇒ 001 ⇒ 001 ⇒　001 ⇒ … A ’ B’ C’ time ｔ ｔ＋１ 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 ｔ ｔ＋１ B C 1 INPUT OUTPUT 000 010 001 101 100 110 011 111

点アトラクター （もしくは、 ） アトラクターの生物学的解釈 異なるアトラクター ⇔ 異なる種類の細胞 点アトラクター 異なるアトラクター　⇔　異なる種類の細胞 点アトラクター 周期が１のアトラクター （静的）定常状態に対応 以下を満たす状態 （もしくは、　　　　　　　　　　　　　　　　）

アトラクターの分布に関する研究 平均入次数を定数（K）とし、ＢＮをランダムに生成した場合のアトラクターのサイズや個数 多くの研究があるが、詳しいことはわかっていない 点アトラクター（非周期的）であれば平均、1個であることが簡単にわかる 計算機実験により、n頂点の場合、　　　　　　個との予想[Kauffman, The Origin of Order, 1993]があるが、正しくなさそう？ [Samuelsson & Troein, PRL, 2003], [Drossel et al., PRL, 2005]

定理：各頂点に対してランダムにn変数ブール関数を選んでブーリアンネットワークを構成した場合、周期１のアトラクターの個数は平均的に１個 証明 時刻 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]

アトラクターが平均１個の理由 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]

再帰型アトラクター検出アルゴリズム（１） 0-１割り当てを１個目の頂点からn個目の頂点までへと順次試していき、矛盾が検出されたら後戻りする [Zhang et al., EURASIP JBSB 2007]

再帰型アトラクター検出アルゴリズム（２） 1 1 1 Output

再帰型アトラクター検出アルゴリズム（３） K 2 3 4 5 6 1.19n 1.27n 1.34n 1.41n 1.45n 0-1割り当てを最初の頂点から順々に探索。点アトラクターにならないとわかった時点でバックトラック ０ ００　－＞　Ｘ ０１ ０１０　－＞　Ｘ ０１１　－＞　Ｘ １０ これは大幅な高速化　（測定値ではなく、理論値） 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細胞では４個の転写因子を導入することにより、多能性幹細胞を構成

先行研究 （Datta et al.） ＮＯ! [Akutsu et al., 2007] （ただし、P≠NPの仮定のもとで） 確率ブーリアンネットワーク （PBN）の制御 いくつかの頂点（遺伝子）を制御可能と仮定 目標状態、初期状態を与え、目標状態にできるだけ近づき、かつ、コストが低い制御系列を計算 古典的制御理論に基づく動的計画法アルゴリズムを提案 ⇒　 2n×2n サイズの行列を扱うことが必要 ⇒　指数時間アルゴリズム 素朴な疑問 多項式時間で制御系列を見つけることは可能か？ ＮＯ!　[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 　　　（時刻０で初期状態、時刻Ｍで目標状態） [Akutsu et al., J. Theo. Biol. 2007]

計算困難性に関する結果 制御系列を見つけるのは一般に計算困難（ＮＰ困難） ネットワークが以下のような単純な構造をしていても計算困難 　　　　　　　⇒　Dattaらの開発したような 　　　　　　　　　　非効率的アルゴリズムを使うのも仕方がない

BN制御の動的計画法アルゴリズム DattaらのアルゴリズムをBNに合わせ簡単化 DPテーブル アルゴリズム（再帰式） 状態　　　　　　　　　　から目標状態に至る制御系列があれば 1、それ以外は０ アルゴリズム（再帰式）

動的計画法アルゴリズムの説明 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表は 指数サイズ

ネットワークが木構造の場合 効率よく（多項式時間で）制御系列が計算可能 木構造：　ループ無しの連結グラフ

木構造に対するアルゴリズム　（１） 基本アイデア: 以下の S[vi,t,b] という表を、 t=0 から t=M という順に、かつ、木の葉から根へという順に、動的計画法で計算 S[vi,t,b]

木構造に対するアルゴリズム （２） S[vi,t,b] の計算: S[vi,t,b] の計算例 木構造に対するアルゴリズム　（２） 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への変換　　 確率ブーリアンネットワークの制御 動的計画法による指数時間アルゴリズム ブーリアンネットワークの制御 木構造に対する多項式時間アルゴリズム