集中講義（九州大学数理学研究院） バイオ構造データに対する数理モデルと アルゴリズム（４） ブーリアンネットワーク

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

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

3. 計算量 • 情報科学では、入力データのサイズ(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くらいなので、実用アルゴリズムを開発できる可能性はある

4. ブーリアンネットワーク

5. ブーリアンネットワーク • 遺伝子ネットワークの数理モデル • 頂点⇔遺伝子 • 状態：　1 (発現) と 0 (非発現)のみ • 制御規則 • ブール関数(AND, OR, NOT …) • もし、y が x を直接制御していれば、y から x に辺が存在 • 状態は（クロックに）同期して変化 • 基本的にデジタル回路と同じ [Kauffman, The Origin of Order, 1993]

6. A time ｔ time ｔ＋１ A ’ B’ C’ A B C 0 0 0 0 0 1 0 0 1 0 0 1 0 1 0 1 0 1 0 1 1 1 0 1 B C 0 1 0 0 0 0 0 1 0 0 A ’ = B 1 0 1 0 1 1 1 0 1 1 0 B ’ = A and C 1 1 1 INPUT OUTPUT C ’ = not A ブーリアンネットワークの例 ブーリアンネットワーク 状態遷移表 状態変化の例： １１１　⇒　１１０　⇒　１００　⇒　０００　⇒　００１　⇒　００１　⇒　００１　⇒　。。。

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

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

9. time ｔ time ｔ＋１ A ’ B’ C’ A B C 0 0 0 0 0 1 0 0 1 0 0 1 0 1 0 1 0 1 0 1 1 1 0 1 0 1 0 0 0 0 0 1 0 0 1 0 1 0 1 1 1 0 1 1 0 1 1 1 INPUT OUTPUT アトラクター(1) • 定常状態 • 細胞の種類に対応？ • 例 • 011 ⇒ 101 ⇒ 010 ⇒ 101 ⇒ 010 ⇒… • 111 ⇒ 110 ⇒ 100 ⇒ 000 ⇒ 001 ⇒ 001 ⇒001 ⇒ …

10. 111 010 000 100 110 011 001 101 time ｔ time ｔ＋１ A ’ B’ C’ A B C 0 0 0 0 0 1 0 0 1 0 0 1 0 1 0 1 0 1 0 1 1 1 0 1 0 1 0 0 0 0 0 1 0 0 1 0 1 0 1 1 1 0 1 1 0 1 1 1 INPUT OUTPUT アトラクター (2)

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

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

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

14. アトラクターが平均１個の理由 Prob( A(1) =1 ) = 0.5 Prob( A(1) =0 ) = 0.5 ランダムに ブール関数を選ぶ t=0 t=1 0 0 A 確率 0.5 0 0 B 確率 0.5 0 0 C 確率 0.5 C B A

15. アトラクターの高速計算 • n頂点の場合、2n個の状態があるので、　　　　　程度の時間かければ計算可能 • 小さい nにしか対応できない • ヒューリスティックな高速アルゴリズムはあるが、計算時間に理論的保証が無い [Irons, Pysica D, 2006], [Devloo et al., Bull. Math. Biol. 2003], … • 点アトラクターの有無の判定はNP完全 [Akutsu et al., GIW 1998] • 平均的に高速なアルゴリズムを開発 [Shuqin et al., EURASIP JBSB 2007]

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

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

18. 再帰型アトラクター検出アルゴリズム（３） • 0-1割り当てを最初の頂点から順々に探索。点アトラクターにならないとわかった時点でバックトラック • ０ • ００　－＞　Ｘ • ０１ • ０１０　－＞　Ｘ • ０１１　－＞　Ｘ • １０ • これは大幅な高速化　（測定値ではなく、理論値） • 230≒1091.230≒240 • 2100=10301.2100≒80,000,000 • 小さい周期のアトラクター検出にも対応可能

19. 点アトラクターの高速計算：計算量解析 t=0 t=1 v1 • 全部でn個の頂点のうちm個の頂点まで0-1割り当てを行った時、vi(0)≠vi(1) が検出できる確率は • よって、m個すべての頂点がアトラクターの条件を満たし、m+1個目の頂点の0-1割り当てがチェックされる確率は • よって、m個の頂点に対する0-1割り当て2m通りのうち、次のステップに再帰呼び出しが進む回数の期待値は • よって、平均計算量は上の式においてmを1からnまで動かした場合の最大値のオーダーとなる vm-1 vm vm+1 K

20. 最悪時の時間計算量 • 入次数 K以下のBNの点アトラクター検出  (K+1)-SAT に変換可能 • K=2 の場合、O(1.322n)時間 • さらに工夫をすれば、O((1.322-δ)n) 時間にすることも可能 • 点アトラクター検出はK=2でもNP困難 • ブール関数を変数とその否定のAND/ORに限った場合にはKに関わらずO(1.587n) 時間で検出可[Melkman, Tamura & Akutsu, 2010]

21. アトラクター検出問題からSATへの変換 • 入次数 K以下のBNの点アトラクター検出  (K+1)-SAT vj vk vi

22. ブーリアンネットワークの制御

23. システム生物学 • 北野宏明氏が提唱e.g., [Kitano, Nature, 420:206-210, 2002] • 生命システムの制御理論の構築（主要目標の一つ） • 線形制御理論は確立しているが、生命システムは非線形 • 生命システムに対する制御理論を構築することが必要 • 創薬やがん治療などにつながる可能性がある • iPS細胞では４個の転写因子を導入することにより、多能性幹細胞を構成

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

25. 確率ブーリアンネットワーク(PBN: Probabilistic Boolean Network) • 確率ブーリアンネットワーク (PBN) • 一つの頂点に対する制御規則が複数 • どの制御規則が選ばれるかは（事前に定められた確率分布に従い）ランダムに決まる • 利点:ノイズなどの影響が扱える。マルコフ過程に関する理論が適用できる • 欠点：ほとんどの場合、O(2n) 時間以上の計算時間がかかるため、大規模なシステムを扱えない 確率0.6で A(t+1) = B(t) AND C(t) A B C 確率0.4で A(t+1) = B(t) OR (NOT C(t))

26. A ’ B’ C’ A B C 0 0 0 0 0 1 0 0 1 0 0 1 0 1 0 1 0 1 0 1 1 1 0 1 0 1 0 0 0 0 0 1 0 0 1 0 1 0 1 1 1 0 1 1 0 1 1 1 ブーリアンネットワークの行列表現 状態遷移表 この行列を　A　とし、 とすると

27. PBNの行列表現 この行列の意味は以下のとおり

28. Dattaらのモデル • 確率ブーリアンネットワークに基づく • 時刻 tにおける全体の状態: • 制御規則：2n×2n行列A • A はm個のブール変数により変化 • コスト関数 • C(x(t), u(t)): u(t)という制御を状態 x(t)に適用した際のコスト • C(x(t)): 最終状態v(t) のコスト

29. Dattaらの制御手法 • 入力：　初期状態x(0), 制御規則A(u(t)), 目的時刻M , コスト関数 • 出力：　コストの合計を最小化する制御系列u(0), u(1), …, u(M-1) • 動的計画法アルゴリズム

30. ブーリアンネットワーク制御モデルの定式化 • 入力 • 内部頂点：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]

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

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

33. 動的計画法アルゴリズムの説明 D[1,1,1, 2] =1 D[0,0,0, 2] = 0 u1=1, u2=1 DP計算 D[0,1,1, 3] = 1 ただし、DP表は 指数サイズ

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

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

36. 木構造に対するアルゴリズム　（２） • S[vi,t,b] の計算: • S[vi,t,b] の計算例 • S[v3,t+1,1]=1 ifS[v1,t,1]=1 andS[v2,t,1]=1 • S[v3,t+1,0]=0 ifS[v1,t,0]=1 orS[v2,t,0]=1 • Note:このアルゴリズムは Datta et al. の 　　動的計画法アルゴリズムと本質的に異なる

37. まとめ • ブーリアンネットワーク • 頂点⇔遺伝子、制御規則⇔ブール関数 • 点アトラクター検出 • 再帰的アルゴリズム（O(2n)時間より高速） • SATへの変換　　 • 確率ブーリアンネットワークの制御 • 動的計画法による指数時間アルゴリズム • ブーリアンネットワークの制御 • 動的計画法による指数時間アルゴリズム • 木構造に対する多項式時間アルゴリズム