九州大学学術情報リポジトリ

(1)

九州大学学術情報リポジトリ

Kyushu University Institutional Repository

システム生物学研究のための並列コンピューティングシステムの設計 : 遺伝子ネットワークの相互作用推定と細胞内情報伝達システムの解析を例にして

岡本, 正宏

九州大学大学院農学研究院生物機能科学部門 | 九州大学大学院システム生命科学府生命情報科学講座 | 九州大学バイオアーキテクチャーセンターシステムデザイン部門

https://doi.org/10.15017/1467693

出版情報：九州大学情報基盤センター広報 : 全国共同利用版. 6 (3), pp.202-210, 2007-03. Computing and Communications Center Kyushu University

バージョン：

権利関係：

(2)

システム生物学研究のための並列コンピューティングシステムの設計：

遺伝子ネットワークの相互作用推定と細胞内情報伝達システムの解析を例にして

岡本正宏

¹

はじめに

分子生物学の発展により、生体内の個々の遺伝子解析からゲノム解析へ、さらにそこから遺伝子セットとしての機能解析が進みつつある。また、同様に細胞内の個々の代謝過程

（酵素反応）から代謝系の解析へ進展し、代謝セット（細胞内タンパク質ネットワーク）

としての機能解析も行われようとしており、生命のシステム論的解析（システム生物学）

あるいはネットワーク解析がますます進むものと期待される。システム生物学は、生物をひとつのシステムとして捉える研究分野と定義される。システム生物学の解析手法としては、１）システム同定・推定(system identification)：構成因子とそのネットワーク構造を推定・同定、２）システム解析(system analysis)：システムのダイナミクスの解析、

３）システム制御(system control)：システムを特定の状態に誘導する制御、４）システム設計(system design)：特定の挙動を再現するシステムの設計、の４つがある。ここでは、

システム同定・推定の例として、大規模遺伝子ネットワークの相互作用推定、システム解析の例として、細胞周期に関わるシグナル情報伝達系を取り上げ、システム生物学研究のために並列コンピューティングシステムがどのように活用されているのか、概説する。

大規模遺伝子ネットワークの相互作用推定

一般に、観測されるシステム要素の動的挙動から、システム要素間の相互作用あるいはネットワークを推定することは、一種の逆問題(inverse problem)である。このことを遺伝子ネットワーク推定問題に当てはめてみる。図

1

に遺伝子ネットワークの概念図を示す。

図

1

において、遺伝子２は、遺伝子

8

の発現を活性化し、同時に、遺伝子

11

の発現を抑制し、また、遺伝子

11

は遺伝子

15

の発現を活性化するといった、遺伝子間の相互作用が既知であるとする。そうすると、もし、遺伝子

2

が強制発現された場合は、遺伝子

8

は活性化、遺伝子

11

と

15

は、発現が抑制され、逆に、遺伝子２が破壊された場合は、遺伝子

8

の発現は抑制、遺伝子

11

と

15

は活性化されることは容易に想像できる。しかし、遺伝子ネットワーク推定問題においては、遺伝子間の相互作用を推定することがタスクであるため、特定の遺伝子を破壊したり、強制発現させたりした時の、他の遺伝子の発現パターンを

DNA

マイクロアレイや

DNA

チップで調べ、遺伝子破壊や強制発現によってどの遺伝子の発現が活性化あるいは抑制されるかを調べ、それに基づいてネットワーク構造を推定しなければならない。つまり、事象結果(resulting event)から相互作用(interaction)を推定・

1

九州大学大学院農学研究院生物機能科学部門，九州大学大学院システム生命科学府生命情報科学講座，九州大学バイオアーキテクチャーセンターシステムデザイン部門

(3)

同定する情報科学的手法が必須となる。遺伝子相互作用推定問題に対してオールマイティーな情報科学的手法はまだ確立されていない。

図 1 遺伝子ネットワーク推定における逆問題

我々は、逆問題解決のための革新的な突破口として、微分方程式の立式に、べき乗表現 (power‑law formalism)に基づいた S‑system モデルを、観測データを再現する多数の実数値パラメータの自動推定法に生物の進化（淘汰、交叉、突然変異）をモデルにした遺伝的アルゴリズムをベースとする手法を導入し、遺伝子ネットワークの相互作用問題に適用してきた(岡本,小野:人工知能学会誌, 18(5), 502‑509 (2003))。ｎ個のシステム構成要素（状態変数）Xi (i=1,2,…,n) の値（濃度、発現量に相当）が時間的に変動し、Xi 同士が相互作用しているネットワークシステムを考える。上記の遺伝子ネットワークの問題であれば、

n 個の遺伝子を考え、Xi は i 番目の遺伝子の発現量を連続値で表している。ここでは Xi の値は非負である。このシステムの各 Xi の速度式（ここでは偏微分項のないｎ元連立常微分方程式を取り扱う）は一般に次式で表される。

(i=1,2,…,n) （１）

式（１）において、右辺第 1 項、第２項はそれぞれ、状態変数 Xi の生成および消費(分解) を表す関数である。通常は、質量作用則に基づいて立式するが、我々が取り扱う系には遺伝子ネットワークのように、反応のメカニズム、相互作用の仕方に関する情報がほとんど得られないものが多い。状態変数の数(n)は既知で、非線形的な相互作用のあるシステムを

(+) (-)

(+)

遺伝子

2

遺伝子8

遺伝子11

遺伝子15 これらの相互作用がわかっていれば

、

● 遺伝子2が過剰発現した場合、

遺伝子8の発現は活性化される遺伝子11の発現は抑制される遺伝子15の発現は抑制される

● 遺伝子2が破壊された場合、

遺伝子8の発現は抑制される遺伝子15の発現は活性化される

遺伝子15の発現は活性化される結果として起こる事象

結果として起こる事象から、相互作用を推定しなければならない (+)

(-)

(+)

遺伝子

2

遺伝子8

遺伝子11

遺伝子15

(+) (-)

(+)

遺伝子

2

遺伝子8

遺伝子11

遺伝子15 これらの相互作用がわかっていれば

、

遺伝子15の発現は活性化される結果として起こる事象これらの相互作用がわかっていれば

、

遺伝子15の発現は活性化される結果として起こる事象

結果として起こる事象から、相互作用を推定しなければならない結果として起こる事象から、相互作用を推定しなければならない

) ,..., , ( )

,..., ,

( ¹ ² ⁿ _i ¹ ² ⁿ

i X X X X X X

dt

dXi ⁼ F ⁺ ^- F ^-

(4)

如何に記述するかが問題である。Savageau はその解決策として、べき乗表現に基づく S‑system 、あるいは BST(biochemical system theory) を提案した (M.A.Savageau,

Biochemical Systems Analysis (Addison‑Wesley, MA, 1976))。それは、式（１）の F

⁺ _i

, F

^‑ _i

を以下のように近似するものである。

(i=1,2,…,n)

（２）

式（２）において、g

_ij

は、状態変数 Xi の生成過程に関与する状態変数 Xj の相互作用係数であり、同様に h

_ij

は、Xi の分解過程（消費過程）に関与する Xj の相互作用係数である。

たとえば、g

_ij

が正の値なら、Xi の生成過程に対し Xj は＋の作用を及ぼし、同様に h

_ij

の値が負なら、Xi の分解過程に対し Xj は−の作用を及ぼすことになる。α

_i

, β

_i

は、それぞれ Xi の生成項、分解項に乗じる係数である。式（２）は、状態変数 Xi の生成過程(右辺第１項)と分解過程(右辺第２項)にシステムを構成しているすべての状態変数 Xj (j=1,2,…,n) が関与していると仮定している。もし、Xi の生成過程（あるいは分解過程）に Xj が関与していない（相互作用がない）場合、g

_ij

（あるいは h

_ij

）の値はゼロということになる。しかし、図２で示すように、S‑system モデルでは、生成過程、分解過程がそれぞれ 1 つの項で表現されているため、生成項、分解項が複数の経路で構成されている場合は、一般質量作用則(generalized mass action law (GMA))を近似した表現になる。逆問題では、観測される状態変数(Xi)の動的挙動（タイムコース）が与えられ、その挙動を説明するような状態変数間の相互作用を表す実数値パラメータを推定しなければならない。つまり、観測データ（Xi のタイムコース）を再現するような式(２)の実数値パラメータ（α

_i

, β

_i

, g

_ij

, h

_ij

）を効率よく（高速かつ高精度に）探索する必要がある。

図２ S‑system モデル

Õ Õ = =

-

= ⁿ

j i j n

j j

i X ^g X ^h

dt

dXi _ij

ij

1 1

a b

• g _ij > 0 : X _j

は

X _i

の合成過程を活性化する

• g _ij < 0 : X _j

は

X _i

の合成過程を抑制する

• g _ij = 0 : X _j

はX

_i

の合成過程に影響を与えない

• h _ij > 0 : X _j

は

X _i

の分解過程を活性化する

• h _ij < 0 : X _j

は

X _i

• h _ij = 0 : X _j

はX

_i

の分解過程に影響を与えない物質Xiの分解過程

1 1

ij ij

n n

g h

i i j i j

j j

d X X X

dt a b

= =

= Õ - Õ

物質Xiの合成過程

X ₁

X ₂

+

X ₃

+ +

遺伝子ネットワーク

１つにまと物質Xiのめられた合成過程

X X _i _i

g _ij h _ij

X X _j _j

S-system

１つにまとめられた物質Xiの分解過程

• g _ij > 0 : X _j

は

X _i

• g _ij < 0 : X _j

は

X _i

• g _ij = 0 : X _j

はX

_i

• h _ij > 0 : X _j

は

X _i

• h _ij < 0 : X _j

は

X _i

• h _ij = 0 : X _j

はX

_i

の分解過程に影響を与えない

• g _ij > 0 : X _j

は

X _i

• g _ij < 0 : X _j

は

X _i

• g _ij = 0 : X _j

はX

_i

• h _ij > 0 : X _j

は

X _i

• h _ij < 0 : X _j

は

X _i

• h _ij = 0 : X _j

はX

_i

の分解過程に影響を与えない物質Xiの分解過程

1 1

ij ij

n n

g h

i i j i j

j j

d X X X

dt a b

= =

= Õ - Õ

物質Xiの

合成過程物質Xiの分解過程

1 1

ij ij

n n

g h

i i j i j

j j

d X X X

dt a b

= =

= Õ - Õ

物質Xiの合成過程

X ₁

X ₂

+

X ₃

+ +

X ₁ X ₁

X ₂ X ₂

+

X ₃

+

X ₃ X ₃

+ +

遺伝子ネットワーク

X X _i _i

g _ij h _ij

X X _j _j

X X _i _i

g _ij h _ij

X X _j _j X X _i _i

g _ij h _ij

X X _j _j

S-system

１つにまとめられた物質Xiの分解過程

(5)

しかし、式(２)では、推定すべきパラメータの個数が非常に多いため(合計 2n(n+1)個)、限られた範囲の時間でモデルを完全に決定できるような実用的な探索アルゴリズムを開発することは容易ではない。我々はこれまで、遺伝的アルゴリズム(Genetic Algorithm、以下 GA と略す)のうち、実数値 GA を基本とする実数値パラメータの最適化法を独自に開発し、

上記の逆問題に適用してきた。その結果、960 個のパラメータを含む 30 遺伝子ネットワーク（相互作用のループ構造を含む）の相互作用を推定することが可能となった。また、生体内の真のネットワーク構造が予めわかっているわけではないので、実験データを再現するネットワーク構造（満足構造とよぶ）を１つではなく、できる限り多数、効率よく探索し、それらをユーザに示すことで、対話的に遺伝子ネットワーク推定をすすめるシステム (interactive inference system of genetic networks)を開発している。さらに、探索をより高速化するために、図３に示す探索用のクライアント・サーバシステムも開発している。

図３マスタ・ワーカモデルを用いた実数値

GA

の並列実装．

個体キューへは

S-system

の各パラメータの符号（＋、ゼロ、−）のみのセットが送られ、クライアント・

スレッドを介して各ワーカ上で、符号を考慮した探索範囲内（実数値）で探索を行い、最良解の評価値を該当する世代交代スレッド内の子個体に通知する。

細胞周期に関わるシグナル情報伝達系のシステム解析

Wet 実験の測定技術の進歩により、多数の遺伝子およびタンパク質などを同時に観測することが可能となり、生体内反応系の動的な振る舞いを効率よく解析する環境が整いつつあ

(6)

る。このようなデータを理論的かつ網羅的に解析するには、遺伝型から表現型にいたる生体内反応系を統合した数理モデルの適用が有効である。しかしながら、従来の生体内反応系のモデル化および動的制御の研究の多くは、遺伝子、代謝、細胞、組織といったある一定の階層を対象にしてきたため、複数の階層を考慮したシステムのシームレスな数理モデルの構築法および解析法の確立が必要となってきた。これら階層の時間および距離といったスケールファクターの値はそれぞれ大きく異なり、また、反応過程（現象）の支配方程式も様々に定義できる。さらに、生体内では階層間の相互作用により、生体内反応系特有の top‑down および bottom‑up 双方向の制御様式を実践している。このような生体内反応系の動特性を理解することは、有用物質の生産方策の設計（バイオ・アーキテクチャー）、疾患の治療方策の検討（オーダメード医療）などに不可欠である。したがって、Wet 実験でえられた膨大な知見を基に、遺伝型から表現型に至る大規模階層構造を如何にモデル化し、

階層内および階層を越えた制御方策を高速に精度よく見出すのかが、今後、生体内階層型統合システムの解析および制御の研究において最重要課題となる。

（１）階層系のシームレスな大規模数理モデルの構築

解析対象として、真核細胞内の細胞周期シグナル伝達機構を取り上げる。生体内マルチスケール・マルチフィジックス現象は、個体差をもつ細胞の集合体（組織）のダイナミクスであり、これをがん組織の発達で考える場合、各細胞内の細胞周期機構の破綻から忠実に議論する必要がある。単一細胞内では、細胞周期は、G1、S、G2、M および G0 期から構成され、細胞質内および核内で進行する反応過程を有する。核内の反応は、分子数が少ないため分子衝突の確率過程を考慮した支配方程式を設計することが望ましい。一方、核内に比して分子数が多い細胞質内の反応は、オルガネラなどの影響により分子が局在化する支配方程式を設計する必要がある。さらに、DNA 損傷、細胞間伝達物質（形質転換因子およびサイトカインなど）により惹起される遺伝子発現量の変化が、細胞周期機構の安定性に大きく影響する（図４参照）。細胞のガン化の数理モデルを構築するには、図４の中で、ガン化に関与するであろう化学種の反応過程および反応場を精査し、適切な支配方程式を設計しなくてはならない。腫瘍細胞は正常細胞に比して、揺らぎまたは時間遅れを創出する不確定性の高い反応系から構成される統合システムであり、細胞ごとに固有のダイナミクスを示す。

（２）階層系の高速・高精度数値計算法の構築

解析対象とする系は、階層構造を持った大規模系であり、階層内の状態方程式（連立微分方程式系）の数学的表記法が階層によって異なり、さらに、状態変数の変化の時間スケールが階層によって異なることを考慮しなければならない。階層によって時間スケールが異なる場合は、時間スケール変換関数を階層間ごとに設定することで、時間スケールの違いを克服する。たとえば、図５において、下位の層（X(t)）は、時間スケールが一番小さ

(7)

図４細胞周期機構に影響するシグナル伝達機構と DNA 修復機構 (D.Hanahan, R.A.Weinberg, Cell, 100, 57‑70, (2000))

朱書きの化学種（タンパク質、遺伝子）は細胞のガン化に関与する。外郭および内郭の膜は、それぞれ細胞膜および核膜を表し、反応場によって支配方程式の設計基準が異なる。また、細胞間を移動する物質（サイトカイン、形質転換因子など）により、単細胞がシステムとして統合され、多細胞系となる。

いため、数値計算の計算刻み幅（可変幅）はかなり小さくなる。しかし、上位の層（Y(t)）

では時間スケールがそれほど小さくないため、下位の層ほど計算刻み幅を小さくする必要が無く、下位のいくつかの計算ステップをまとめて１タイム刻みにしてよい。上位の層の適切な計算刻み幅を決定する指標として、下位の層の暫定区間の微分平均値と最終微分値を比較して、両者の差が許容誤差範囲内ならば、暫定値を確定値に設定し、暫定区間の確定値の積分平均値を算出し、それを上位の層の値算出に適用する。なおその場合、上位の層の計算刻みは暫定区間となる。許容誤差範囲内でなければ、範囲内になるまで、暫定区間を小さくする。このアルゴリズムを基本に、下位から上位への多層の時間スケールを数値的に調整する。高速高精度数値計算法については、時間刻み可変の Gear 法に加えて、新たに、計算条件の難度に応じて自動的に項数を変え、かつ桁落ち誤差を最小限に押さえる適切な刻み幅を迅速に決定できる Taylor 級数法を導入する。

（３）階層系の高速・高精度計算のためのハードウェア構築

支配方程式を時間スケールごとにクラスタリングし数値計算を実施する。master と

worker

から構成される数値シミュレータの計算機システムを図６に示す。数理モデルの常

(8)

図５階層系の数値計算のための時間スケール関数の適用

微分方程式群（o）、確率微分方程式群（s）、偏微分方程式群（p）を時間スケールごとに例えば５段階（クラスタ１〜５）に分割し、階層を定義する。クラスタ１に属する方程式群の反応速度が最も速く（数値計算の時間刻み幅が狭い）、クラスタ５に属する支配方程式群の反応速度が最も遅い（CPUの負荷が最も軽い）。クラスタ１〜５の計算結果は数値積分の繰返し演算が完了するたびに

master

の共有メモリに書き込まれ（Write）、必要に応じて他のクラスタの数値積分計算で化学種の濃度（値）が引用（Retrieve）される。共有メモリには、クラスタ５の数値積分の任意の繰返し演算回数分の数値計算結果を保存する。また、

master

は、クラスタ５の数値積分の繰返し演算が完了するたびに、ハードディスクへクラ

スタ１〜５の直近の計算結果を格納（Save）する。ここで、

CPU

から共有メモリへの

Write

および

Retrieve

の排他制御をいかに効率化するかが問題となる。

多細胞系の数値シミュレーションでは、臨床上

1cm

を超えたがん組織が治療対象となることから、100×100個の二次元配列された細胞シート（2cm相当）の計算条件を確保したい。多細胞系の数値シミュレーションの計算機システムを図７に示す。多細胞系の数値シミュレーションでは、細胞数

n

個に相当する単細胞系の計算機システムを

n

セット準備し、

さらに、経時的に細胞数が増加するため、単細胞系の計算機システムを自己増殖的に発生させなくてはならない。

最後に、図６、図７において、ハードウエア構成の見積り（希望も含めて）を行ってみる。

時間スケール関数の適用

下位の層の変化量を基に上位の層の適切な刻み幅を決定下位の層の変化量を基に上位の層の適切な刻み幅を決定

各階層の適切な時間刻み幅が必要各階層の適切な時間刻み幅が必要

X 0

下位の層：

X (t)

確定値

Y 0

上位の層：

Y(t)

tim e

X 0

下位の層：

X (t)

確定値

Y 0

上位の層：

Y(t)

tim e

時間スケールが異なる複数の階層の数値計算を同時に進行

tX1

暫定値暫定値

tX2 tX3 tX4 tX5 tX6 tX7 tX8 tX9

Y 1

tX10

If （

−

tX’10

）＞許容誤差

then

暫定区間を縮小して、

[tX 0,..,tX9]

で再試行

'

X

If （ t

−

tX’10

）＞許容誤差

then

[tX 0,..,tX9]

で再試行

'

X t

暫定区間の微分値の積分平均値

( ) t X '

と最終微分値

(tX’9 )

を算出暫定区間の微分値の積分平均値

( ) t X '

と最終微分値

(tX’9 )

を算出

If （

−

tX’9

）＞許容誤差

then

[tX 0,..,tX8]

で再試行

'

X

If （ t

−

tX’9

）＞許容誤差

then

[tX 0,..,tX8]

で再試行

'

X t

暫定区間暫定区間

暫定区間の微分値の積分平均値

( ) t X '

と最終微分値

(tX’8)

を算出暫定区間の微分値の積分平均値

( ) t X '

と最終微分値

(tX’8)

を算出

X tX t

X1 X2 X 3 X 4 X5 X6 X7 X8

X1 X2 X 3 X 4 X5 X6 X7 X8 tX9 tX9 tX9 tX1 0 tX 11 tX1 0 tX 11 tX1 0 tX 11

If （

−

tX’8

）＜許容誤差

then

暫定値を確定値に設定

暫定区間の確定値の積分平均値

( )を算出

を

Y1

の算出に適用する

Y1を用いて tX9

から数値計算を再開

' X t

X X t

t

If （

−

tX’8

）＜許容誤差

then

( )を算出

を

Y1

Y1を用いて tX9

' X

If （ t

−

tX’8

）＜許容誤差

then

( )を算出

を

Y1

Y1を用いて tX9

' X t

X X t

t

(9)

CPU p1 CPU s1 CPU o1

クラスタ１

Write Retrieve

ハードディスク

Save

共有メモリー

(Time= t-4 , τ , Τ ) (Time= t-3 ,τ,Τ) (Time= t-2 ,τ,Τ) (Time= t-1 ,τ,Τ) (Time= t -1 ^,τ,Τ) CPU pk

CPU sk CPU ok

クラスタｋ

CPU p5 CPU s5 CPU o5

クラスタ５

CPU X (master) CPU p1

CPU s1 CPU o1

クラスタ１

Write Retrieve

Save

共有メモリー

(Time= t-4 , τ , Τ ) (Time= t-3 ,τ,Τ) (Time= t-2 ,τ,Τ) (Time= t-1 ,τ,Τ) (Time= t -1 ^,τ,Τ)

共有メモリー

(Time= t-4 , τ , Τ ) (Time= t-3 ,τ,Τ) (Time= t-2 ,τ,Τ) (Time= t-1 ,τ,Τ) (Time= t -1 ^,τ,Τ) CPU pk

CPU sk CPU ok

クラスタｋ

CPU pk CPU sk CPU ok

クラスタｋ

CPU p5 CPU s5 CPU o5

クラスタ５

CPU p5 CPU s5 CPU o5

クラスタ５

CPU X (master)

図６単細胞系の数値シミュレータの計算機システム

数理モデルの常微分方程式群（o）、確率微分方程式群（s）、偏微分方程式群（p）を時間スケールごとに例えば５段階（クラスタ１〜５）に分割し、階層を定義する。Masterは、各クラスタの数値計算結果を共有メモリに保存し、クラスタ５の数値積分の繰返し演算が完了するたびにハードディスクへクラスタ１〜５の計算結果を格納する。なお、クラスタ１〜５に含まれる

CPU

は

worker

である。

C1 C2

C3 C#

C#

C# C#

C#

Cm

Cn

C#

CPU p1 CPU s1 CPU o1 クラスタ１

Write Retrieve

共有メモリー (Time= t-4 ,τ,Τ) (Time= t-3 ,τ,Τ) (Time= t-2 ,τ,Τ) (Time= t-1 ,τ,Τ) (Time= t-1,τ,Τ) CPU pk

CPU sk CPU ok

クラスタｋ CPU p5 CPU s5 CPU o5

クラスタ５

CPU X (master) CPU p1

CPU s1 CPU o1 クラスタ１

共有メモリー (Time= t-4 ,τ,Τ) (Time= t-3 ,τ,Τ) (Time= t-2 ,τ,Τ) (Time= t-1 ,τ,Τ) (Time= t-1,τ,Τ) 共有メモリー (Time= t-4 ,τ,Τ) (Time= t-3 ,τ,Τ) (Time= t-2 ,τ,Τ) (Time= t-1 ,τ,Τ) (Time= t-1,τ,Τ) CPU pk

クラスタｋ CPU pk CPU sk CPU ok

クラスタ５ CPU p5 CPU s5 CPU o5

クラスタ５

CPU X

(master)

C1

C2

C3 C#

C#

C# C#

C#

Cm

Cn

C#

CPU p1 CPU s1 CPU o1 クラスタ１

共有メモリー (Time= t-4 ,τ,Τ) (Time= t-3 ,τ,Τ) (Time= t-2 ,τ,Τ) (Time= t-1 ,τ,Τ) (Time= t-1,τ,Τ) CPU pk

クラスタ５

CPU X (master) CPU p1

CPU s1 CPU o1 クラスタ１

共有メモリー (Time= t-4 ,τ,Τ) (Time= t-3 ,τ,Τ) (Time= t-2 ,τ,Τ) (Time= t-1 ,τ,Τ) (Time= t-1,τ,Τ) 共有メモリー (Time= t-4 ,τ,Τ) (Time= t-3 ,τ,Τ) (Time= t-2 ,τ,Τ) (Time= t-1 ,τ,Τ) (Time= t-1,τ,Τ) CPU pk

クラスタｋ CPU pk CPU sk CPU ok

クラスタ５ CPU p5 CPU s5 CPU o5

クラスタ５

CPU X (master)

図７多細胞系の数値シミュレータの計算機システム

ｎ個の細胞からなる多細胞系の数値シミュレーションは、図８の単細胞系の計算機システムをｎセット準備する。共有のハードディスクに全細胞の計算結果を格納し、グラフィカル・ユーザ・インターフェースを用いて計算結果を可視化する。

1 単細胞系の数値シミュレータ（図６）の共有メモリの容量

図６の計算機システムに、化学種 200 個、反応速度定数 800 個から構成される数理モデルを適用し、クラスタ１〜５の各階層間に数値積分計算の繰返し回数が 10 倍

(10)

の差をもつと仮定する。化学種 200 個に倍精度浮動小数点型の変数（double）を適用した場合、メモリの最大必要容量は 160MB（=200 個×8byte×10

⁵

）である。そして、

反応速度が最も遅いクラスタの数値積分に着目し、10 回分（10 time step）の計算結果を保存すると、メモリの最大必要容量は 1.6GB（=160MB×10）となる。さらに、

反応速度定数 800 個に要するメモリ、計算機システムに費やすメモリなどを考慮すれば、化学種 200 個、反応速度定数 800 個から構成される単細胞系の数値シミュレーションに必要な共有メモリの容量は 2GB 弱となる。

2 多細胞系の数値シミュレータ（図７）の共有ハードディスクの容量

単細胞 10000 個の、図７の数値シミュレータの反応速度が最も遅いクラスタの数値積分計算 100 回分（100 time step）の計算結果を共有ハードディスクに格納する場合、200TB 弱（160MB×10000 個×100）の容量を必要とする。

おわりに

システム生物学の

4

つの研究戦略のうち、システム同定・推定とシステム解析について、

並列コンピューティングシステムがどのように活用されるのか概説した。

第一の大規模遺伝子ネットワークの相互作用推定については、大規模数値最適化を如何に効率よく行うかが最大のポイントである。これまで、ヒューリスティックな探索で、遺伝的アルゴリズムを代表とする進化アルゴリズムが数多く提唱されてきたが、最適解が予めわかっている種々のベンチマーク関数を設定し、最適解を如何に早く探索できるかでアルゴリズム評価を行っている。しかし、ここで紹介したように、生物の問題は、どれが最適解か特定できず、その評価は実験的（分子生物学的）実証にゆだねるしかない。したがって、情報科学の立場からは、如何に質のよい解候補をなるべく多く（もれなく）拾い上げるようなアルゴリズムを作るかが課題である。

第二の階層性をもつ生体内反応系のマルチスケール・マルチフィジックス現象の計算基盤の構築および数値シミュレーションは、国内外で未だ実践されていない。また、細胞周期の一部分は数理モデル化されつつあるが、細胞周期全体を表現する数理モデルは構築されていない。もちろん、腫瘍細胞のダイナミクスを微視的に表現する数値シミュレーションも実現していない。したがって、このような研究内容は、構築する計算モデルおよびそのアプリケーションのいずれにおいても非常に独創性が高く、様々な用途において将来性のあるテーマであることから早期着手すべき課題と考える。しかし、ハードウエアの見積りでも述べたが、かなり大規模の計算システムを設計しなければならず、予算的にも一研究室でできるようなものではなく、情報基盤センターのような計算センターがコアになって推進するテーマであろう。

九州大学学術情報リポジトリ