Transcript 1 - 京都大学
集中講義（九州大学数理学研究院）
バイオ構造データに対する数理モデルと
アルゴリズム（４）
ブーリアンネットワーク
阿久津 達也
京都大学 化学研究所
バイオインフォマティクスセンター
内容
ブーリアンネットワーク
定常状態の高速検出アルゴリズム
確率ブーリアンネットワーク
制御系列計算のアルゴリズム
代謝ネットワークのブーリアンモデルと
頑健性解析
計算量
情報科学では、入力データのサイズ(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
A’ = B
B’ = A and C
C’ = not A
time ｔ
time ｔ＋１
A B C
A’ B’ C’
0
0
0
0
1
1
1
1
0
0
1
1
0
0
1
1
0
0
1
1
0
0
1
1
0
1
0
1
0
1
0
1
INPUT
0
0
0
0
0
1
0
1
1
1
1
1
0
0
0
0
OUTPUT
状態変化の例：
１１１ ⇒ １１０ ⇒ １００ ⇒ ０００ ⇒ ００１ ⇒ ００１ ⇒ ００１ ⇒ 。。。
ブーリアンネットワーク研究の意義
単純化しすぎとの批判
でも、単純化しないと理論解析、推定、制御などが困難
定性的理解には有用（？）
最も単純な非線形システムの一つ
シミュレーションなら複雑なモデルでも可能
ブーリアンネットワークで不可能なら、他の（多くの）非
線形システムでも不可能
本質的にデジタル回路と同じ
計算機科学（特にブール関数、アルゴリズム理論）にお
ける成果や技法が利用可能
定常状態の高速計算アルゴリズム
アトラクター (1)
定常状態
細胞の種類に対応？
例
011 ⇒ 101 ⇒ 010 ⇒
101 ⇒ 010 ⇒…
111 ⇒ 110 ⇒ 100 ⇒
000 ⇒ 001 ⇒ 001 ⇒
001 ⇒ …
time ｔ
time ｔ＋１
A B C
A’ B’ C’
0
0
0
0
1
1
1
1
0
0
1
1
0
0
1
1
0
0
1
1
0
0
1
1
0
1
0
1
0
1
0
1
INPUT
0
0
0
0
0
1
0
1
1
1
1
1
0
0
0
0
OUTPUT
アトラクター (2)
time ｔ
time ｔ＋１
A B C
A’ B’ C’
0
0
0
0
1
1
1
1
0
0
1
1
0
0
1
1
0
0
1
1
0
0
1
1
0
1
0
1
0
1
0
1
INPUT
0
0
0
0
0
1
0
1
1
1
1
1
0
0
0
0
OUTPUT
111
011
110
100
000
001
101
010
点アトラクター
アトラクターの生物学的解釈
異なるアトラクター ⇔ 異なる種類の細胞
点アトラクター
周期が１のアトラクター
（静的）定常状態に対応
以下を満たす状態
v(t ) (v1 (t ),, vn (t ))
v(t 1) v(t )
（もしくは、
v(1) v(0)
）
アトラクターの分布に関する研究
平均入次数を定数（K）とし、ＢＮをランダムに生
成した場合のアトラクターのサイズや個数
多くの研究があるが、詳しいことはわかっていない
点アトラクター（非周期的）であれば平均、1個であるこ
とが簡単にわかる
計算機実験により、n頂点の場合、 O( n ) 個との予
想[Kauffman, The Origin of Order, 1993]があるが、正しくなさそう？
[Samuelsson & Troein, PRL, 2003], [Drossel et al., PRL, 2005]
入次数＝２
A
入次数＝３
A
定理：各頂点に対してランダムに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]
アトラクターが平均１個の理由
B
Prob( A(1) =1 ) = 0.5
A
ランダムに
ブール関数を選ぶ
C
000 000 確率 0.53
t=0 t=1
A
B
C
0
0
0
0
0
0
Prob( A(1) =0 ) = 0.5
確率 0.5
確率 0.5
確率 0.5
001 001 確率 0.53
111 111 確率 0.53
8 0.53 1
アトラクターの高速計算
n
頂点の場合、2n
n
個の状態があるので、 O(2 )
程度の時間かければ計算可能
小さい n にしか対応できない
ヒューリスティックな高速アルゴリズムはあるが、計算
時間に理論的保証が無い
[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個目の頂点までへと
順次試していき、矛盾が検出されたら後戻りする
ProcedureIdentAttractor( v, m)
if m n 1 the noutput v; re tu rn;
for all b {0,1} do
v m b;
if fi ( v) v i is found for some i m the ncon tin ue
;
e lseIdentAttractor( v, m 1);
re tu rn.
[Zhang et al., EURASIP JBSB 2007]
再帰型アトラクター検出アルゴリズム（２）
0
0
0
0
Output
1
0
0
0
1
0
0
1
再帰型アトラクター検出アルゴリズム（３）
0-1割り当てを最初の頂点から順々に探索。点アトラクター
にならないとわかった時点でバックトラック
これは大幅な高速化 （測定値ではなく、理論値）
０
００ －＞ Ｘ
０１
０１０ －＞ Ｘ
０１１ －＞ Ｘ
１０
230≒109
2100=1030
1.230≒240
1.2100≒80,000,000
小さい周期のアトラクター検出にも対応可能
K
2
3
4
5
6
平均計算量
1.19n
1.27n
1.34n
1.41n
1.45n
点アトラクターの高速計算：計算量解析
全部でn個の頂点のうちm個の頂点まで0-1割り当
てを行った時、vi(0)≠vi(1) が検出できる確率は
C
m
0.5 m K 0.5
n
n CK
よって、m個すべての頂点がアトラクターの条件を
満たし、m+1個目の頂点の0-1割り当てがチェック
される確率は
K m
m
1
0
.
5
n
よって、m個の頂点に対する0-1割り当て2m通りの
うち、次のステップに再帰呼び出しが進む回数の
期待値は
m
m
2 1 0.5
n
m
t=0 t=1
K
K
よって、平均計算量は上の式において m を1からn
まで動かした場合の最大値のオーダーとなる
v1
vm-1
vm
vm+1
K
最悪時の時間計算量
入次数 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]
アトラクター検出問題からSATへの変換
入次数 K 以下のBNの点アトラクター検出
(K+1)-SAT
vj
vk
vi
ブーリアンネットワークの制御
システム生物学
北野宏明氏が提唱 e.g., [Kitano, Nature, 420:206-210, 2002]
生命システムの制御理論の構築（主要目標の一つ）
線形制御理論は確立しているが、生命システムは非線形
生命システムに対する制御理論を構築することが必要
創薬やがん治療などにつながる可能性がある
iPS細胞では４個の転写因子を導入することにより、多
能性幹細胞を構成
先行研究 （Datta et al.）
確率ブーリアンネットワーク （PBN）の制御
いくつかの頂点（遺伝子）を制御可能と仮定
目標状態、初期状態を与え、目標状態にできるだけ近づ
き、かつ、コストが低い制御系列を計算
古典的制御理論に基づく動的計画法アルゴリズムを提案
⇒ 2n×2n サイズの行列を扱うことが必要
⇒ 指数時間アルゴリズム
素朴な疑問
多項式時間で制御系列を見つけることは可能か？
ＮＯ!
[Akutsu et al., 2007]
（ただし、P≠NPの仮定のもとで）
[Datta et al., Machine Learning, 2003]
確率ブーリアンネットワーク
(PBN: Probabilistic Boolean Network)
確率ブーリアンネットワーク (PBN)
一つの頂点に対する制御規則が複数
どの制御規則が選ばれるかは（事前に定められた確率
分布に従い）ランダムに決まる
利点:ノイズなどの影響が扱える。マルコフ過程に関す
る理論が適用できる
欠点：ほとんどの場合、O(2n) 時間以上の計算時間が
かかるため、大規模なシステムを扱えない
B
確率0.6で A(t+1) = B(t) AND C(t)
A
C
確率0.4で A(t+1) = B(t) OR (NOT C(t))
ブーリアンネットワークの行列表現
状態遷移表
A
0
0
0
0
1
1
1
1
B
0
0
1
1
0
0
1
1
C
0
1
0
1
0
1
0
1
A’
0
0
1
1
0
0
1
1
B’
0
0
0
0
0
1
0
1
C’
1
1
1
1
0
0
0
0
0 0
1 1
0 0
0 0
0 0
0 0
0 0
0 0
0
1
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
1
0
0
1
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
1
0
0
0
0 1
0 0
0 0
0 0
0 0
0 0
1 0
0 0
000
001
010
011
100
101
110
111
この行列を A とし、 v(t ) (v1 (t ),, vn (t ))
とすると v(t 1)T A v(t )T
PBNの行列表現
0
0
0 0.7 0 0.1 0 1
0 0
0
0
0
0
0 0
0.6 0.6 0.7 0
0 0
0
0 0.2 0 0.8 0
0 0
0
0
0 0.3 0
0.4 0.4 0 0.1 0
0 0 0.3 0
0
0
0
0
0
.
9
0
0 0
0 0.8 0.8 0 0.2 0
0 0
0
0
0
0
0
0
.
2
0
0
0
.
6
0
0 0
0
0
0
.
1
0
0
.
1
0
0
0
.
1
000
001
010
011
100
101
110
111
この行列の意味は以下のとおり
Aw, u Pr(v(t 1) u | v(t ) w )
Dattaらのモデル
確率ブーリアンネットワークに基づく
x(t ) ( x1 (t ),, xn (t )) {0,1}n
時刻 t における全体の状態:
制御規則： 2n×2n 行列 A
A は m 個のブール変数により変化 u(t ) (u1 (t ),, um (t ))
Pr(x(t 1) b | x(t ) a ) A(u(t ))
コスト関数
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)
M
E C (x(M )) C (x(t ), u(t ))
t 0
s.t.
Pr(x(t 1) | x(t )) A(u(t ))x (t ),x (t 1)
動的計画法アルゴリズム
J T (x( M )) C (x( M ))
111
J t (x(t )) min m C (x(t ), u(t )) A(u(t ))x (t ),y J t 1 (y )
u ( t ){0 ,1}
y 000
ブーリアンネットワーク制御モデルの定式化
入力
内部頂点： 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テーブル D[b1 , b2 ,, bn , t ]
状態 [b1 , b2 ,, bn ] から目標状態に至る制御系列
があれば 1、それ以外は０
アルゴリズム（再帰式）
1, if [b1 ,, bn ] v M
D[b1 , b2 ,, bn , M ]
0, otherwise
1, if thereis (c, u) such that
D[b1 , b2 ,, bn , t 1]
D[c1, ,cn ,t] 1 and c f (b, u)
0, otherwise
動的計画法アルゴリズムの説明
D[b1 , b2 ,, bn , t 1]
D[1,1,1, 2] =1
D[0,0,0, 2] = 0
1, if thereis (c, u) such that
D[c1, ,cn ,t] 1 and c f (b, u)
0, otherwise
u1=1, u2=1
DP計算
D[0,1,1, 3] = 1
ただし、DP表は
指数サイズ
ネットワークが木構造の場合
効率よく（多項式時間で）制御系列が計算可能
木構造： ループ無しの連結グラフ
木構造
木でない構造
木構造に対するアルゴリズム （１）
基本アイデア: 以下の S[vi,t,b] という表を、 t=0
から t=M という順に、かつ、木の葉から根へと
いう順に、動的計画法で計算
S[vi,t,b]
1, if thereexist u(0), u(1),, u(t ) such thatvi (t ) 1
S [vi , t ,1]
0, otherwise
1, if thereexist u(0), u(1),, u(t ) such thatvi (t ) 0
S [vi , t ,0]
0, otherwise
木構造に対するアルゴリズム （２）
S[vi,t,b] の計算:
1, if thereexist [bi1 , , bik ] such that f (bi1 , , bik ) 1 holds
S [vi , t 1,1]
and S [vi j , t , bi j ] 1 holds for all j 1, , k
0, otherwise
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への変換
確率ブーリアンネットワークの制御
頂点⇔遺伝子、制御規則⇔ブール関数
動的計画法による指数時間アルゴリズム
ブーリアンネットワークの制御
動的計画法による指数時間アルゴリズム
木構造に対する多項式時間アルゴリズム