卒

題

業

論

文

目

モンテカルロ法を用いた化学反応モデルの FPGA によるシミュレーション
年度

平成 17 年度
所

属

慶應義塾大学
理工学部

情報工学科

指導教員

天野 英晴 教授氏

名60026142

吉見 真聡

卒業論文要旨
学科

情報工学

学籍
番号

60026142

フリガナ

ヨシミ マサト氏名

吉見 真聡(論 文 題 名 )モンテカルロ法を用いた化学反応モデルの FPGA によるシミュレーション

^
内
容
の
要
旨
は
約
25
行
程
度
で
記
入
の
こ
と
_

(内 容 の 要 旨)

近年の分子生物学の急速な進展により，種々の生物の遺伝子配列や，細胞内の代謝系などに関する定量的データが蓄積されつつある．これらデータを用いて細胞や個体を計算機上でシミュレートし，実験を行うことなく様々な条件下における生物の分子レベルでの挙動を解明できるようになることが期待されている．
生物学実験は、倫理的な問題が絡むこと，個体の成長に時間を要することなどか
ら，比較的古くから計算機を利用したシミュレーションについて研究が行われてきた．そのため近年では、そのいくつかが実用化され実際に利用されている．
しかし，計算機の処理能力は急速に向上しているにもかかわらず、それらシミュ
レーション実行に要する時間は，実時間の数倍から数万倍にも及び，計算時間に関わる問題は解決が困難である．細胞シミュレーションという問題は非常に不規則で複雑であるため，SMP やクラスタなどを用いた並列計算処理による高速化は適していない．
その一方で，分子動力学や流体力学といった分野では，専用計算機が性能や処理速度の面で大きな成果を挙げてきており，生物学の分野でも，ハードウェア処理による高速化への期待は大きい．しかし，生物のシミュレーションは様々なアルゴリ
ズムが組み合わせて使用されるため，専用ハードウェアの設計を行うことは難しい．
そこで，近年活発化している FPGA などを用いた可変構造型ハードウェアによる計算システム (リコンフィギャラブルシステム) を応用することで解決を図る．回路を再構成することでアルゴリズムの柔軟性を確保し，問題を専用計算機のように直接ハードウェア処理することで高速化を実現する．
本研究では，生物学計算指向リコンフィギャラブルプラットフォームである ReCSiP
をターゲットとして，Gillespie の開発した確率モデルシミュレーションアルゴリズ
ムを実装し，シミュレーションによる評価を行った．その結果，Lotka 反応モデルにおいて AMD 社 AthlonXP2800+を用いたソフトウェアの実行時の 53.48 倍の高速化を実現した．

慶應義塾大学理工学部

目次
第1章

1

緒論

第 2 章 背景
2.1 Gillespie のアルゴリズム . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.1 解析的モデルでの解法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.2 確率的モデルによる化学運動学の定式化のための物理的基礎 . . . . . . . . .
2.1.3 確率的シミュレーションアルゴリズム . . . . . . . . . . . . . . . . . . . . .
2.1.4 確率モデルの利点と制限 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.5 Lotka モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 主な細胞シミュレータ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 Virtial Cell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 STOCKS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.3 STOCHSIM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3
3
4
4
10
11
12
13
13
17
20

第 3 章 ReCSiP の構成
3.1 ReCSiP の設計方針 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 シミュレーション高速化のためのストラテジ . . . . . . . . . . . . . . . . . . . . . .
3.3 ReCSiP の構成 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 全体像 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.2 ReCSiP ボード . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.3 Verilog ライブラリ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.4 ソフトウェア . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23
23
23
24
24
25
28
31

第 4 章 設計と実装
4.1 変数の定義 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 シミュレーションの流れ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 実装 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.1 実装モジュール . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.2 浮動小数点実数の算術演算器 . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.3 一様乱数生成器 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.4 対数テーブル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.5 Lotka 反応モデルの実装 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.6 変数保存テーブル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.7 入出力モジュール . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34
34
34
36
36
37
38
39
40
45
47

i

第 5 章 評価
5.1 ソフトウェアでの計算時間 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Vertex-II への実装結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3 シミュレーション結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.1 ハードウェアとソフトウェアの動作比較 . . . . . . . . . . . . . . . . . . . .
5.4 性能評価 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

49
49
50
52
54
55

第 6 章 結論
6.1 確率モデルを FPGA でシミュレーションする有効性 . . . . . . . . . . . . . . . . . .
6.2 今後の課題 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2.1 除算器の実装による並列化 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2.2 精度向上 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2.3 モデル自動生成 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

57
57
57
57
57
59

参考文献

61

図目次
2.1
2.2
2.3
2.4
2.5
2.6

分子の衝突 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
衝突体積 ∆Vcoll . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
仮定する反応 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
確率モデルシミュレーションアルゴリズム . . . . . . . . . . . . . . . . . . . . . . . 11
finite-volume formalism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
Virtual Cell による神経細胞内のカルシウム濃度と信号の変化 . . . . . . . . . . . . . 17

3.1
3.2
3.3
3.4
3.5
3.6

ReCSiP の構成 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ReCSiP ボード . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ReCSiP ボードの構成 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ローカルバスの read/write 動作 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
BlockRAM コントローラ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ReCSiP の浮動小数点演算器 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25
25
26
28
30
31

4.1
4.2
4.3
4.4
4.5
4.6

組み込み乗算器を使用した浮動小数点乗算器 . . . . . . . . . . . . . . . . . . . . . .
32 ビット LFSR による疑似乱数生成モジュール . . . . . . . . . . . . . . . . . . . .
Lotka 反応モデルの確率モデルシミュレーションの計算過程 . . . . . . . . . . . . .
シミュレーション 1 つが 1 サイクル中でモジュールを使用する部分 . . . . . . . . .
37 段パイプラインの Lotka 反応モデルモジュールが使用する部分 . . . . . . . . . .
Lotka 反応モデルモジュールを 2 つ使用したシミュレーションモジュール . . . . . .

38
39
41
44
46
47

5.1
5.2
5.3
5.4
5.5
5.6
5.7
5.8

ソフトウェアで出力される分子数変化 . . . . . . . . . . . . . . . . . . . . . . . . . .
ソフトウェアの乱数種を変えたときの変化 . . . . . . . . . . . . . . . . . . . . . . .
Software result1 の末端 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Software result2 の末端 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
実装した Lotka 反応セットモジュールのシミュレーション結果 . . . . . . . . . . . .
並列動作する Lotka 反応シミュレーションの一部 . . . . . . . . . . . . . . . . . . .
ソフトウェアによる実行結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ハードウェアによる実行結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

50
51
51
51
53
54
55
55

6.1
6.2

ホストと協調動作する対数テーブル . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
より周期の長い乱数生成器の一例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

iii

表目次
3.1
3.2
3.3
3.4
3.5
3.6

ReCSiP ボードの構成部品 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ローカルバスの信号線 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
標準ローカルバスインターフェイスのレジスタ割当 . . . . . . . . . . . . . . . . . .
標準ローカルバスインターフェイスのメモリアドレス割当 . . . . . . . . . . . . . .
浮動小数点演算器の仕様 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
市販 IP コアの浮動小数点演算器の仕様 . . . . . . . . . . . . . . . . . . . . . . . . .

4.1
4.2

浮動小数点乗算器の比較 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
パイプライン化 Lotka 反応モジュールに必要な演算器 . . . . . . . . . . . . . . . . . 45

5.1
5.2
5.3

ソフトウェアでの処理能力 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
入出力付き Lotka 反応セットの合成結果 . . . . . . . . . . . . . . . . . . . . . . . . 52
各モジュールの合成結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

iv

26
27
29
30
31
31

第1章

緒論
1980 年代後半からの分子生物学の急速な進展により，種々の生物の遺伝子配列や，細胞内の代謝系などに関する定量的データが蓄積されつつある．これらのデータを用いて，細胞や個体を計算機上でシミュレートすることにより，実験を行なうことなしにさまざまな条件下における生物の分子レベルでの挙動を知ることができるようになることが期待されている．
実験と計算機シミュレーションの反復という研究手法の有効性は，既に他の多くの分野において有効性が実証されつつある．生物学の分野では，まだモデル化に必要なデータが決して充分とはいえず，シミュレーション自体はまだ研究途上の課題であるといえるが，生物学実験には常に倫理的な問題が絡むことや，個体の生長に時間を要することなどから，特に高等生物を対象にした実験は困難であり，シミュレーションへの期待はきわめて大きい．
したがって，計算機を利用してシミュレーションを行う，という研究は比較的古くから行われており，1983 年に開発された Kinsim[1] に汎用シミュレータのはじまりを見ることができる．1980
年代から 1990 年代にかけて Gepasi[2], DBsolve[3] などに代表される，代謝系のモデル化・シミュ
レーション・パラメータ最適化を行うソフトウェアが開発されるようになり，1990 年代後半になると計算機の性能向上や入手可能なデータの増加に伴って，細胞内の代謝反応すべてをシミュレー
ションするための E-Cell[4], The Virtual Cell[5] といったシミュレータの開発が行われはじめている．
また，代謝系を微分方程式にモデル化して解析的にシミュレーションを行う解析モデルに対して，
化学反応システムで発生する反応を確率を使って計算する確率モデルが 1977 年に Gillespie[6] によって提案された．1990 年代後半より，細胞の生長と分裂の挙動をシミュレーションする STOCKS[7]，分子の立体構造や共有結合の規則変化をシミュレーションする STOCHSIM[8] といった，確率モ
デルで代謝系へアプローチするシミュレータが開発されてきている．
しかし，これらのシミュレータを使用して細胞全体の代謝活動をシミュレーションするために要する時間は実時間の数倍から数万倍にも及んでおり，計算時間の問題は深刻である．その一方
で分子動力学や流体力学といった分野では，専用計算機が性能や成果の面で大きな成果を挙げてきており，ハードウェア処理による高速化への期待は大きい．
生物のシミュレーションにおいてはさまざまなアルゴリズムが組み合わせて使用され，新しいアルゴリズムも次々に導入されるなど，専用ハードウェアの設計は難しいと考えられる．しかし，可変構造型ハードウェアによる計算システムは，問題を専用計算機のように直接ハードウェア処理することによって高い性能を発揮することが可能できる上に，回路をソフトウェア的に書き換える柔軟性を兼ね備えており，生物学分野への応用が有望であると考えられる．
本研究では，FPGA を搭載した生物学計算向けハードウェア ReCSiP(Reconfigurable Cell Simulation Platform)[9] をターゲットとして，確率モデルによる化学反応モデルのシミュレーションア
ルゴリズムを実装し，その性能評価を行った．本論文の第 2 章では，化学反応システムに対する1

1. 緒論
確率モデルでのアプローチの方法と，近年開発が進められている細胞シミュレータについて述べ，
第 3 章では細胞シミュレーションプラットフォーム ReCSiP の構成と機能について述べる．第 4 章
では，ReCSiP をターゲットとした化学反応システムの確率モデルシミュレーションの一例として
Lotka 反応モデルのシミュレーションを高速化する手法と実装について検討し，第 5 章でその性能評価の結果について述べ，第 6 章で結論と今後の展開について述べる．

2

第2章

背景
近年，生物学の分野では遺伝子ネットワーク，代謝経路などを同定するためのデータを高速に得るための高スループットな実験方法の開発が盛んになっており，生物に関する数値的なデータを得ることが容易になりつつある．90 年代に始まった各種ゲノムプロジェクトの成果が出始めて
おり，2002 年にはヒト，イネ等の比較的遺伝資料の大きな生物の遺伝子の高精度解読の完了が宣言された．また，ゲノムと同様，様々な生物の代謝経路も研究されている．
さらに，これらの実験から得られた大量のデータを元に，計算機を用いて網羅的な解析や数理モデリングとシミュレーションを行う動きも本格化している．
この章では，細胞内代謝反応に代表される化学反応システムのシミュレーションの手法について，確率モデルによるアプローチである Gillespie のアルゴリズムについて述べる．次に，解析モ
デルと確率モデルについて，主なシミュレータの現状について述べる．

2.1 Gillespie のアルゴリズム
ある空間上の化学反応システムについて，その挙動を数学的に記述する 2 つの方法がある．
解析モデル (deterministic approach) は，システムの時間の変化を，微分方程式の組で決定されるものと見なす方法である．
また，確率モデル (stochastic approach) は，システムの時間の変化を，単一の微分差分方程式に支配されるランダムウォーク過程の一種と見なす方法である．確率モデルの計算はマスタ方程式と呼ばれる確率・統計的なアプローチから開発された手法であるが，マスタ方程式は数学的に解答することが困難である．しかし，マスタ方程式を直接的に扱うことなく，それと等価な数値計算を行う方法を用いることで，計算機を使った計算を行うことができる．
Gillespie の開発した “stochastic simulation algorithm” [6] (以下 Gillespie のアルゴリズム) は，この確率モデルを用いて化学反応システムの挙動を計算するアルゴリズムである．
Gillespie のアルゴリズムは次のような問題に解答することを目的として開発された．
¶³

N 種の化学物質が一様に分布する一定の体積 V の中で，相互反応する M 種類の化学反応経路があり，ある初期時間に各化学物質の分子数が与えられた場合，これらの分子数は時間の経過とともにどのように変化するか．
µ´

32. 背景2.1.1

2.1. Gillespie のアルゴリズム

解析的モデルでの解法

この問題を解くために，従来，問題を常微分方程式で表現して解析する解析モデルでのアプロー
チが使用されてきた．
ある時刻 t における V 内の i 番目の化学物質の数 (分子集合レベル．一般的には分子数など) が，
時刻 t を変数とする通常非線形の関数 Xi (t) (i = 1, . . . , N) の式で記述されている化学反応システムを仮定する．このシステムで M 種類の各化学反応が連続的に速度変化すると見なせるならば，以下のような一次の常微分方程式を構成することができる．

dXi
dt= fi (X1 , · · · , XN )(i = 1, · · · , N)

(2.1)

右辺の関数 fi は，反応経路 M の構造と反応速度定数によって決定される．これらの方程式は，反応速度方程式 (reaction-rate equation) と呼ばれる．X1 (t), . . . , XN (t) の解を導くことが，各化学物
質の時間変化になる．
このような反応速度方程式の解析的な解法は，通常単純なシステムで利用されることが多いた
め，計算機でこれらの方程式を数値的に計算することが必要になる．
この方法は，微分方程式で近似していることから，化学反応システムの時間変化が連続的かつ決定的であることを暗黙に仮定している．しかし，分子数のような分子集合レベルは離散的な整数値であるから，化学反応システムは時間的に非連続なプロセスである．さらに，発生する反応は決定的なものではない．もし，分子の挙動がなんらかの方程式に支配されるものと見なしたとしても，システム内のすべての分子の位置と速度を計算しない限りは，システムの挙動を計算す
ることは不可能である．
式 2.1 のような方程式を解く方法は，いくつかの化学反応システムの解析には非常に良い精度を示している．しかし，解析モデルによるシミュレーションは，システムの時間挙動を微分方程式の組で記述するため，初期値の小さな違いが，シミュレーションの進行とともに指数関数的に拡大してしまう．また，遺伝子発現に代表されるような，ごく短い時間で小数の分子が急激に反応する場合には，反応発生にランダムな要素が含まれ，適切なシミュレーションができない．
生態系や微小な生化学システム，非線形なシステムの挙動に注目が集まっている現在においては，分子集合レベルを十分高精度に計算できる方法が必要になる．
Gillespie のアルゴリズムは，化学反応システムが離散的かつ確率論的な挙動を示す，という前提から，化学反応システムの挙動を計算する計算機解法である．

2.1.2

確率的モデルによる化学運動学の定式化のための物理的基礎

化学反応は，適当な 2 種類以上の分子が衝突したときに発生する．確率的方法は，熱平衡状態の分子反応システムで起こる分子の衝突を計算するものである．

分子の衝突ある体積 V 内で，平衡状態にある 2 種類の気体分子 S 1 と S 2 の混合物からなるシステムがあるとする．単純化のため，S 1 分子と S 2 分子はそれぞれ半径 r1 ，r2 の球体であると考える．S 1 分子と S 2 分子の中心が重なったときに必ず反応が発生し，半径 r12 = r1 + r2 の分子になるとする (図
2.1)．

4

2. 背景2.1. Gillespie のアルゴリズム

r1

r1

S1

r2

S1

r1

r2

S2

S1

r2

S2

S2

Collision !

r1

r1

S1

r 1+r 2

r2

r2

S1
S2

S12

S2

図 2.1: 分子の衝突

V 内でこのような衝突が発生する頻度 (反応の速度) を計算する．
従来の反応速度を導出する方法を以下に述べる．まず，任意の分子の対から反応が起こる分子を選択する．そして，S 2 分子に関して S 1 分子の相対的な運動速度 v12 を求める．さらに，微小時
間 ∆t の間に，S 1 分子が S 2 分子に対して移動する分の体積 (Collision Volume : 衝突体積) を計算する．衝突体積 ∆Vcoll は以下の式 2.2 で求められる (図 2.2)．
2∆Vcoll = πr12
· v12 ∆t

(2.2)

S1

v12 ∆tS2

r1

r12 =r1 +r2

r22∆V =πr12 v12 ∆tcoll

図 2.2: 衝突体積 ∆Vcoll

もし，時間 t で S 2 分子の中心が ∆Vcoll 内に存在していたとすると，2 つの分子は微小時間 (t, t+∆t)で衝突する．

5

2. 背景2.1. Gillespie のアルゴリズム

この方法は，微小時間 ∆t を小さくしていくことで，反応速度の精度は上がることになる．また，
極限 ∆Vcoll → 0 なので，衝突体積中に存在する分子は 0 か 1 になってしまう．このように，微小時間で計算していく方法は，より厳密であるとは言い難い．
このことから，以下のように確率を用いた方法を考える．
システムは熱平衡状態であるから，分子は V 内で一様に分布している．そのため，時刻 t において任意の S 2 分子の中心が ∆Vcoll にある確率は，単純に ∆Vcoll /V で与えられる．よって，S 1 分
子と S 2 分子の反応速度分布を平均した場合，以下のような式 2.3 が得られる．これは，システム全体の体積に対する衝突体積の割合に等しい．

V

coll


=

V

2 v ∆tπr1212= S 1 -S 2 の分子が次の微小時間 ∆t で衝突する確率の平均
V

(2.3)

マクスウェルの速度分布より，相対速度の平均 v12 は (8kT/πm12 )1/2 と等しい．k はボルツマン
定数，T は絶対温度，m12 は換算質量 m1 m2 /(m1 + m2 ) である．上記の計算は S 1 分子，S 2 分子がV 内にそれぞれ 1 個ずつある場合であるが，S 1 分子の数として X1 ，S 2 分子の数として X2 が与えられたとすると，式 2.3 は，式 2.4 のように表される．

V

coll

V


=

2 v dtX1 X2 πr1212= 微小時間 (t, t + dt) で S 1 -S 2 の衝突が発生する確率
V

(2.4)

式 2.4 により，反応発生数を確定的には計算できないが，有限時間内に V 内で発生する衝突の確率を計算することができる．式 2.4 は決定的な反応速度過程の代わりに，確率的なマルコフ過程を構成しており，過去のシステムの挙動には依存せず，その時点での分子数と “単位時間あたりの衝突確率” によって，熱平衡状態の分子を特徴付けるものである．

反応確率速度定数 — “The Stochastic Reaction Constant”
化学反応の “反応速度” は，“適当な分子の単位時間あたりの 衝突確率 = 単位時間あたりの反応発生確率”と表現することができることを前節で述べた．
V 内で S 1 分子と S 2 分子が式 2.5 のように反応すると仮定する．

R1 :

S 1 + S 2 → 2S 2

s1

(2.5)

s2
Collision !

s2s2
図 2.3: 仮定する反応6

2. 背景2.1. Gillespie のアルゴリズム

式 2.3 より，2 つの分子の物理的性質とシステムの温度だけに依存する定数 c1 を設定する．

c1 dt = 微小時間 dt で特定の S 1 -S 2 分子対が作用 (反応 R1 が発生) する確率

(2.6)

もし，V 内の時間 t において，S 1 分子の分子数が X1 ，S 2 分子の分子数が X2 存在するならば，以下の式が成立する．

X1 X2 c1 dt = 反応 R1 が微小時間 (t, t + dt) で V 内のどこかで発生する確率(2.7)

式 2.7 を一般化して，熱平衡状態にある体積 V 内の N 種類の S i 分子が Xi の数だけ存在する混合物について， M 種類の化学反応 Rµ が相互作用すると仮定すると，分子の物理的性質とシステ
ム温度に依存する式 2.8 で表される M 種類の定数 cµ を仮定できる．

cµ dt = 反応 Rµ に関わる分子が無限微小時間 dt で衝突する確率の平均(2.8)

式 2.8 における “平均” とは，その反応 Rµ に関係する分子の組み合わせの総数に cµ を掛けた値
が，反応 Rµ が次の微小時間 (t, t + dt) で V 内のどこかで発生する確率になることを意味している．

反応速度定数と反応確率速度定数の関係
反応確率速度定数 cµ は反応速度定数 kµ と密接に関係している．式 2.5 のようなある反応 R1 について，以下のような関係がある．

k1 =Vc1 hX1 X2 i
hX1 ihX2 i(hxi は x の平均値)(2.9)

右辺について，hX1 X2 i = hX1 ihX2 i と見なせるので，

k1 = Vc1(2.10)

もし反応 R1 が 2 つではなく 3 つの反応分子を持つならば，V の代わりに V 2 を使用する．また，もし反応 R1 がただ一つの反応分子を持つ (isomerization : 異性化) ならば，V = 1 になる．
次に，反応 R1 の逆反応を以下の式 2.11 から考える．

R2 :2S 1 → S 1 + S 2(2.11)式 (2.5) から，定数 c2 を決定する (S 1 の分子の特定の組が微小時間で衝突し，R2 の反応を発生する確率の平均)．V 内の S 1 分子の組み合わせの数は X1 (X1 − 1)/2! になる．

Vc2 h X1 (X2!1 −1) i Vc2k2 ==hX1 ihX1 i2(2.12)

確率的モデルの時間計算

M 種類の相互反応をする N 種類の分子からなるシステムの時間経過を計算する．これは基礎的な仮説である式 2.6 から求めることになる．

7

2. 背景2.1. Gillespie のアルゴリズム

マスター方程式
確率モデルの化学反応システムの時間経過計算は，システムのマスタ方程式を設定し，それを解くことで求められる．マスタ方程式を数学的に解答することは困難で，シミュレーションでは後に述べるマスタ方程式を解くことと同義の処理を行ってシステムの時間経過を計算する．
マスタ方程式の最も重要な要素は主要確率関数 (grand probability function) である．化学反応シ
ステムにおいて主要確率関数は式 2.13 で表される．

P(X1 , X2 , · · · , XN ; t) ≡V 内の時刻 t で，S 1 分子が X1 存在し，
(2.13)かつ S 2 分子が X2 存在し · · · ，かつ S N 分子が XN 存在する確率式 2.13 は，時刻 t のシステムの状態確率 (stochastic state) を表している．たとえば，Xi に関する
k 次の P のモーメント (k 乗平均) は，式 2.14 で与えられる．

Xi(k) (t) ≡

∞
X

...

X1 =0

∞
X

Xik P(X1 , · · · , XN ; t)(2.14)

XN =0(i = 1, · · · , N; k = 0, 1, 2, · · · )Xi に関する k 次の P のモーメントは，“時刻 t で V にある分子数 Xi の平均の k 乗” である．ここでの “平均” とは，2.6 で定義された確率モデルでの平均反応発生確率を使い，同じ初期値 (Xν )で，時間 0 から t まで非常に多くのシミュレーションを実行を繰り返したときの Xi の平均のこと意味している．
毎回のシミュレーションで S i 分子の分子数 Xi は異なるが，これらの値の k 次の累乗の平均値(k)は，シミュレーションを非常に多く実行することで，Xi に収束していく．k = 1, k = 2 で表され(k)
る 1 次モーメント，2 次モーメントは特に有用である．k = 1 での Xi は時間 t で V 内に存在する
S i 分子の数の平均を示しており，k = 2 では式 2.15 を計算することで，母平均まわりの 2 次モー
メント (分散) を求めることができる．
Xi(k) (t) ≡

∞
X
X1 =0

...

∞
X

Xik P(X1 , · · · , XN ; t)(2.15)

XN =0(i = 1, · · · , N; k = 0, 1, 2, · · · )(1)

解析的な反応速度の式 2.1 で現れている Xi (t) は，1 次のモーメント Xi を近似しているが，完全な等式になることは少ない．
マスタ方程式は P(X1 , · · · , XN ; t) で表される時間経過の関数である．P(X1 , · · · , XN ; t) は，システ
ムが時刻 t + dt で状態 (X1 , · · · , XN ) になる 1 + M 通りの方法があり，それぞれの確率の合計である．この P は，式 2.6 から確率の加法と乗法を使用して求められる．


M
M

 X
X




P(X1 , · · · , XN ; t + dt) = P(X1 , · · · , XN ; t) 1 −aµ dt +Bµ dt(2.16)µ=1µ=1

ここで，ある数量 aµ を定義する．

aµ dt ≡ cµ dt × { 状態 (X1 , · · · , XX ) で反応 Rµ が起こる分子の組み合わせ数 }= 時刻 t で状態 (X1 , · · · , XN ) のシステムで，(t, t + dt) において Rµ 反応が発生する確率
8(2.17)

2. 背景2.1. Gillespie のアルゴリズム

式 2.16 の第 1 項は，システムが時刻 t において状態 (X1 , · · · , XN ) になり，時刻 t + dt においてその状態 (反応が起こらない) のままでいる確率である．第 2 項の Bµ dt は，時刻 t でシステムが状態(X1 , · · · , XN ) になり，時刻 t + dt で Rµ の反応を受ける確率である．
式 2.16 からマスタ方程式の導出をすることができる．
M h
X
i
∂
P(X1 , · · · , XN ; t) =Bµ − aµ (X1 , · · · , XN ; t)∂tµ=1(2.18)

特別な場合を除いて，マスタ方程式の記述は簡潔であるが，それを解くことは困難であることが多い．さらに，反応速度方程式とは異なり，独立変数の数と性質が原因で，計算機を使って数値的に解を導くことも容易ではない．全ての反応 Rµ が単純な単分子反応でなければ，モーメン
トを導出するための式は常に高次のモーメントを含んでいるため，特定の Xi や ∆i (t) のモーメン
トの時間変化を方程式で表して計算しようとすることは，無限大の時間が掛かってしまう．マスタ方程式は正確で簡潔に表現することはできるが，実際的な数値計算を行う上では有用ではない．
そのため，確率モデルでは，結果的にマスタ方程式と同じ解答を得られる方法を用いて計算機でのシミュレーションが可能な方法を用いる．

反応確率分布関数
化学反応システムについて，確率的方法を用いた場合に，時間経過を計算する方法について述
べる．
時刻 t で状態 (X1 , · · · , XN ) のシステムがあるとすると，このシステムの時間経過を計算するには，いつ，どの反応が起こるのかを明らかにする必要がある．
そのため，以下の式 2.19 で定義される関数 P(τ, µ) を導入する．

P(τ, µ) ≡ V 内の状態 (X1 , · · · , XN ) で，微小時間 (t + τ, t + τ + dτ) の間に Rµ の反応が起こる確率(2.19)
関数 P(τ, µ) は反応確率分布関数と呼ばれる．τ は次の反応が起こるまでの時間，µ は発生する反応の種類を示している．
この化学反応システムにおいて発生しうる反応 Rµ について，式 2.20 で定義される関数 h を計算する．
hµ ≡ 状態 (X1 , · · · , XN ) において，Rµ に関係する分子の組み合わせ数(µ = 1, · · · , M) (2.20)

もし反応 Rµ が S 1 + S 2 → S という形式であるならば，hµ = X1 X2 となり，Rµ が 2S 1 → S という形式ならば，hµ = 12 X1 (X1 − 1) になる．一般的に，hµ は変数 X1 , · · · , XN の複合関数である．
状態 (X1 , · · · , XN ) において，時刻 (t, t + dt) で反応が起こらない確率を P0 (t) とする．また，Rµが時刻 (t + τ, t + τ + dτ) で発生する確率を aµ とすると，aµ は式 2.21 で表される．

aµ dτ ≡ hµ cµ dτ(2.21)

よって，P(τ, µ) は以下のように定義される．

P(τ, µ)dτ = P0 (τ) · aµ dτ

9(2.22a)2. 背景2.1. Gillespie のアルゴリズム

P
右辺の P0 (τ) の導出について， 1 − ν aν dτ0 が状態 (X1 , · · · , XN ) の時間 dτ0 で反応が発生しな
い確率である事を利用し，以下の式から P0 (τ) を得る．


M
X


P0 (τ0 + dτ0 ) = P0 (τ0 ) · 1 −aν dτ0 (2.22b)
ν=1

 M

 X

0
P0 (τ) = exp −aν dτ (2.22c)

ν=1

式 2.22c を式 2.22a に代入することで，反応確率密度関数 P(τ, µ) を計算することができる．
aµ exp(−a0 τ) if 0 6 τ < ∞ and µ = 1, · · · , M(2.23)P(τ, µ) = 0otherwise

2.1.3

確率的シミュレーションアルゴリズム

本節では，化学反応システムにおいて，いつ，どの反応が発生するのかを計算する方法を述べる．
これは，“単位間隔の一様な乱数” から反応確率密度関数 P(τ, µ) の乱数を作り出して τ, µ を計算することで明らかになる．(0, 1) の一様乱数 r1 ，r2 から，以下の式を用いて τ と µ を生成する．

 
τ = a10 ln r11
µ−1
X
ν=1

aν < r2 a0 6µ
X

(2.24a)aν(µ : Integer)(2.24b)

ν=1

この 2 式より P2 (µ) = aµ /a0 に従う整数の乱数 µ と，確率密度関数 P1 (τ) = a0 exp(−a0 τ) に従う乱数 τ が生成され，P1 (τ) · P2 (µ) = P(τ, µ) を求めることができる．
シミュレーションは Step 1. から Step 3. の繰り返しで進行する (図 2.1.3)．

Step 0. 初期化
M 種類の反応定数 c1 , · · · , c M と，N 個の初期分子数 X1 , · · · XN に初期値を格納する．時刻 t と反応数カウンタ n を 0 に初期化する．乱数生成器を初期化する．
Step 1.
この時点での分子の数ごとに， M つの値 a1 = h1 c1 , · · · , a M = h M c M を計算し，a0 に aν の合計値を代入する．

Step 2.
(0, 1) の一様乱数 r1 と r2 を生成し，τ と µ を計算する．

10

2. 背景2.1. Gillespie のアルゴリズム
Step 0. Initialization
set X i c it =0set Random generator

Step 1. calculateax

set h i X
set a i = h ic i
set a 0= Σa i

Step 2. generate r1 & r2
calculate τ & µ
set r 1 = random1
set r 2 = random2
calculationτ
calculationµ

Step 3. Reactions
set τ = t + τoccur reaction µXi

図 2.4: 確率モデルシミュレーションアルゴリズム

Step 3.
τ によって t を増加する．反応 Rµ が影響を与える分子の数を調節する (例えば，Rµ が式 2.5 の反応だったならば，X2 が 1 増加し，X1 が 1 減少する)．反応数カウンタ n をインクリメントする．
これを Step 1. に戻って繰り返す．

Step 1. から Step 3. の繰り返しの中のどこか，あるいは任意の t または n の間隔において，(X1 , · · · , XN , t) の値を出力する．また，t あるいは n が何らかの値に到達したときや，a0 が 0 になったときには計算を終了する．

2.1.4

確率モデルの利点と制限

確率的モデルのシミュレーションアルゴリズムの利点と制限を述べる．
このシミュレーションアルゴリズムは，基礎的な仮定 (式 2.6) の形式で記述された化学反応シス
テムについて，正確なシミュレーションが可能である．また，確率モデルのシミュレーションア
ルゴリズムは，従来の微分方程式の数値解法のように，微小時間 ∆t ごとに計算を繰り返す方法をとらず，“次の反応までの時間” を計算するため，分子数が急激に変化するようなシステムにも有効である．このシミュレーションアルゴリズムを計算機上で実行した場合，メモリスペースの要
求量が少ないことも利点の一つとして挙げられる．N 種類の分子と M 種類の反応では，N 種類の11

2. 背景2.1. Gillespie のアルゴリズム

分子数と M 個の cν ，M + 1 個の aν の変数を格納するメモリスペースが確保できれば良い．このことは，FPGA 上に確率モデルシミュレーションアルゴリズムを実装する上でも大きな利点となる．
また，シミュレーションの過程や結果を読み取ることが容易であることも利点のひとつである．
シミュレーションサイクルごとに分子数が数値的に評価できるため，平均，分散，相関などを計算することが容易である．
確率モデルのシミュレーションアルゴリズムは，反応数に比例した計算時間が掛かる．そのため，シミュレーションできる分子の種類，数，反応数は適切な数を指定し，計算時間を適切に抑えなければならない．シミュレーションが可能なシステムは分子が一様分布している状態に限ら
れ，分子数に偏りがある場合は，システムの大きさを変更するなどして対応する必要がある．また，このシミュレーションアルゴリズムから信頼できる結果を得るためには，信頼性のある乱数生成器を使用しなければならない．この信頼性に関する標準を得ることは非常に困難である．さらに，良い統計精度を得るためには，複数回のシミュレーションの評価を取る必要があり，その実行回数は計算時間などのコストとのトレードオフになる．

2.1.5

Lotka モデル

化学反応の確率モデルシミュレーションを開発した Gillespie は，確率モデルシミュレーションを適用した例として，4 つのモデルについて言及している．それは，不可逆異性化反応，Lotka 反応，Blusselator，Oregonator の 4 つである．不可逆異性化反応は逆反応の起こらない不可逆変化する分子の反応であり，Gillespie はこの反応のシミュレーションを用いて確率モデルの具体的な検証を行っている．その他の 3 つのモデルはより複雑なモデルとして挙げ，検証している．不可逆
異性化反応は，後に述べる Lotka モデルの一部として使われる基本的なモデルである (式 2.25c のみで記述されるモデルが不可逆異性化反応である)．
今回の実装では，Gillespie の述べた “より複雑なモデル” から，Lotka 反応モデルを ReCSiP 上に実装して高速化を図ることを考える．
Lotka 反応モデルは，1920 年に Lotka が以下のような自触媒反応を観察したことから研究が始まった．
c1X¯1 + X2 −→ 2X2(2.25a)c2

X2 + X3 −→ 2X3
c3

X3 −→ Z

(2.25b)(2.25c)

このモデルは，Volterra が開発した捕食者-被食者関係に対応する生態系の数学的モデルから研究されているモデルである．Volterra は上記の反応を以下のような微分方程式で記述し，これに対応する反応速度方程式の使用法について研究した．

dX1= c1 X1 X2 − c2 X2 X3
dt(2.26a)dX3= c2 X2 X3 − c3 X3
dt(2.26b)12

2. 背景2.2. 主な細胞シミュレータ

式 2.25b は捕食種 X3 が被食種 X2 を食べて自分を再生産する様子を表現している．また，式 2.25aは X2 が食料 X1 を食べて自分を再生産する様子を記述している．ここでは X1 はただ消耗されるだけで，増加しないものと仮定している．異性化式 2.25c は X3 の死を表現している．
定常状態は dX2 /dt = dX3 /dt = 0 の条件を満たしているシステム状態であることは明らかである．
Lotka 反応モデルは化学反応の単純なモデルの一例であるが，これは以下のような単純なルールに支配される生態系を考えることでどのような反応であるかを容易に理解することができる．ここでの Lotka 反応モデルのシミュレーションに際しては，前述の式 2.25 に加え，X2 の死，およびX1 は反応によって減少しない，という条件を追加している．
¶

³

• 平原
– 一定の広さの平原があり，草が繁茂している場所が点在している– 平原の外に出ることはできず，また平原の中に何かが入ってくることもない• ウサギ– ウサギは，草を見つけるとそれを食べ，繁殖して 2 匹に増える– オオカミに食べられたり，一定の確率で自然死する• オオカミ– オオカミは，ウサギを見つけるとそれを食べ，繁殖して 2 匹に増える– オオカミは，一定の確率で自然死するµ´

この平原にウサギ 1000 匹，オオカミ 1000 匹を一様に分散させて投入したところ，時間の変化
とともにそれぞれの数はどのように変化するか，を計算するモデルが Lotka モデルの単純な形であると考えることができる．一般的に，草がウサギに食べられても減少しないという条件のもとで，適切な cµ を設定することで，ウサギとオオカミの数は震動することが知られている．

2.2

主な細胞シミュレータ

現在までに開発された，あるいは開発中の細胞シミュレータについての調査結果を述べる．まず，現在では解析モデルでのシミュレータが多く開発されているため，代表的な解析モデル細胞シミュレータである Virtual Cell について，および確率モデルのシミュレータである STOCKS，
STOCHSIM について述べる．

2.2.1

Virtial Cell

概要
Schaff らが開発した Virtual Cell[5] は単一の細胞の構造とその内部の反応をシミュレーションする汎用的なフレームワークである．Virtual Cell は細胞のメカニズムをテストする総合的なシステ
ムを提供するとともに，細胞内の運動と分布を総合的に扱うフレームワークを提供する．
Virutual Cell は細胞内で発生する個々の反応について，生化学的なデータ，あるいは電気信号な
どの電磁気学的なデータを微小な画素で表現し，細胞の時間変化を画像で表示する．Virtual Cell
13

2. 背景2.2. 主な細胞シミュレータ

はあらゆる分子のメカニズムを速度方程式や膜輸送の記述で定義し，個々の反応を計算機上に集約させる目的で開発が続けられている．

Virtual Cell が実行する計算
Virutual Cell は単一の細胞を微小な立方体 (element と呼ばれる) に分割し，微小時間ごとにそれぞれの element ごとに解析的な計算を行う，という方法を取っている．
また，Virtual Cell には以下の特徴がある．
• シミュレーションする問題は偏微分方程式で記述され，変数はシステムの状態を表し，次の微小時間のシステムへの入力にもなる．
1. 拡散/電位勾配による分子の流れ
2. 荷電した分子の分布による電束密度3. cable 理論によって求められる電流密度4. 系から生成される物質
5. 膜を介した物質の輸送
上記の項目 1〜4 は全ての element について，項目 5 は膜 (elemet を構成する物質がそれ以外と大きく異なっている部分) のみで定義される．

• 基本的な物理法則の計算式が実装され，シミュレーションごとに再定義する必要が無い• 変数は element の位置と状態 (分子濃度と電位) の関数で表される→−→−
項目 1 の分子の流れ Ji は細胞の挙動の解析に用いられる．Ji は NernstPlank 方程式によって，式2.27 記述される．

→−
Ji = −Di ∇Ci − zi µiCi ∇Vi
Ci
µi
R

濃度
移動度
気体定数

Di
V
T

µi =

拡散係数
ポテンシャル
絶対温度

xi
F

Di F
RT

(2.27)

イオン価数
ファラデー定数

式 2.27 の右辺第 1 項は拡散泳動を，第 2 項は電場の流れを示している．第 2 項は急速に反応が起こる場合や，高解像度の空間を解析する場合に使われる．
式 2.27 から，ある微小体積内における分子 i の濃度の時間変化は式 2.28 で表される．

δCi→−= −∇ · Ji + S i (C1 · C2 · · · C N )δt
S は分子の生成・消費に関わる複数の反応から総合的に決定される反応速度である．
単一の反応では，S i は速度方程式によって与えられる．
k1

A+BC
k−1

14

(2.28)(2.29a)2. 背景2.2. 主な細胞シミュレータ

S C = k1 [A][B] − k−1 [C](2.29b)分子は複数の反応で生成・消費されるので，S i の値は個々の反応速度の合計を与える．
微小時間ごとに微小体積上に定義される濃度は，分子が element 内で生成・消費される量と，element 間で流動する量との合計である．
項目 2 と 3 は，電子を持つ分子の分布と，移動度ベースの電位を記述する高解像度解析用の冗長的な要素である．項目 5 は膜を通過する分子を制御する法則である．2 つの隣接 element で性質が異なっている部分を膜と定義し，膜内外の分子について，双方向移動が記述されている．

シミュレーションの流れ
Virutual Cell のシミュレーションの目的は，なんらかの外乱を与えた後の恒常性維持機構の挙動を調べることにある．そのためには濃度や電位などについて適切な初期条件を与える必要がある．
定常状態を初期値とすることで，外乱に対する挙動を調べやすくすることができる．初期条件は実験的に概算されているものといないものがある．適当な値を定常状態になるまでシミュレーショ
ンすることで定常解を求めることができるが，初期微動を最小化する方法を使用すれば，非線形最適化問題として解くことができるようになる．
Virtual Cell が行うシミュレーションは以下のような流れで進行する．
まず，実験から得られた 3 次元データから，以下の項目の複合で細胞の形状を分析幾何学的に近似する．
• 球体，円柱のような平面プリミティブ• フラクタルソリッドのような数学的立体モデル• ポリゴン面 (CAD モデル)
シミュレーションは化学的な現象を分子，次元ごとの element 数，element のサイズ，シミュレー
ション時間などを定義して柔軟に描き出すことが可能である．分子の運動は，element 内の濃度の変化と，element の境界を横切る分子の流れの合計として描写される．分子の運動は，分子の拡散，イオン化傾向，電位，化学反応，分子の空間分布などの関係を表している．膜は，異なる性質を持つ element の境界で表現される．
シミュレーションの出力になるものは，各 element の時間ごとの分子の濃度，電位，膜上に存在している分子の状態や密度などである．
膜反応モデルと element 内のモデルはモデルの変数や各種パラメータの関数として記述された
反応速度のセットで構成されており，これらはデータベースに格納され，必要に応じてモデルオ
ブジェクトライブラリを介して呼び出される．この DB 化により，ユーザは毎回数値モデルに手を加える必要から解放される．

finite-volume formalism
Virtual Cell は finite-volume formalism と呼ばれる計算方法を利用している．
シミュレーションは計算対象空間を一定サイズの微小な立方体 (element) に区切って計算が行わ
れる．
15

2. 背景2.2. 主な細胞シミュレータ

式 2.30 で分子の流れを表現する (2.5)．
I
→
− →
−
−
Ji · d S = Areaboundary · [JX |X+ + JX |X− + JY |Y+ + JY |Y− + JZ |Z+ + JZ |Z− ]i(2.30)

S

隣接する element が膜だった場合やシミュレーションの端だったりした場合は，その方向の面で膜輸送の処理など特別な処理を行う．

Jx |x+
Jx | x-yzx

図 2.5: finite-volume formalism
式 2.31 のように，反応速度は微小時間で線形な部分と定数の合計であると近似して計算を行う．

S A ≈ − fA linear ([A], [B], · · · ) · [A] + fA cons tan t ([C], [D], · · · )(2.31)

未知の特性については，定数や適当なモデルで近似する．
膜反応は，膜上のチャネルが選択的に膜内外の分子を透過させるような処理を行う．チャネル
は複数の状態を持ち，膜内外の電位差や分子の濃度によって制御される．状態遷移を化学量論的に記述することによって，膜上の分子密度を記述するだけで膜反応シミュレーションを行えるようになる．

実行例 — Calcium wave in neural cell
Virtual Cell の実践テスト例として，カルシウムの生化学レベル，電気生理学レベルでの相互作用のメカニズムが提示されている．カルシウムは細胞内で信号を伝達するセカンドメッセンジャー
であり，その機嫌と規則と役割を解析することは大きな研究課題である．
実際の神経細胞におけるカルシウムを媒介とした信号伝達は以下の手順で進行する．
1. イノシトール三リン酸 (IP3 ) が神経細胞内で高濃度になる
2. 小胞体が IP3 に刺激され，Ca2+ を放出する
Virtual Cell で IP3 と IP3 R(IP3 受容体) の運動を視覚化する (図 2.6)．これは生物実験的に干渉することが難しい部分である．図 2.6 の上方の神経細胞の画像はカルシウムの濃度を，下方は IP3の濃度の分布を表している．このような図は任意のシミュレーション時間間隔で出力され，動画として確認することも可能である．このシミュレーションで実験で見られるものと比較して，非常に類似している Ca2+ の動きのパターンと速度が得られた．

Virtual Cell の計算時間
図 2.6 に示された神経細胞のシミュレーションは，反応のモデルと分子の数を制限して適度な2 次元セルモデルに近似して計算が行われている．計算は 27000 の element と 3 つの非線形偏微
分方程式，7 つの非線形常微分方程式から構成される．評価は Silicon Graphics 社の R8000MIPS
CPU(300Mflops) のワークステーションを使用した．2 秒間のシミュレーションに 2 日を要している．

16

2. 背景2.2. 主な細胞シミュレータ

図 2.6: Virtual Cell による神経細胞内のカルシウム濃度と信号の変化

Virtual Cell の展開
Virtual Cell は現在バージョン 3.1 であり，より多様なモデルに対応し続けている．また，シミュレーションしたデータからリポジトリを構築し，World Wide Web で共有されている．
Virtual Cell は Japa applet で記述され，プロジェクトの Web ページから Virtual Cell を使用することが可能になっている(注 1) ．

2.2.2

STOCKS

Kierzek は原核生物の遺伝子発現を解析するソフトウェア STOCKS(STOChastic Kinetic Simulation) を開発した．STOCKS は標準 C++で記述されたソフトウェアプログラムである(注 2) STOCKS
は化学反応系を Gillespie のアルゴリズムを利用して解析する．

STOCKS の特徴
STOCKS は Gillespie が提案していた計算方法に加え，細胞内代謝シミュレーションを可能とするために主に以下の 2 つの機能が実装されている．
• 体積増大および細胞分裂に対応
STOCKS では，細胞を化学反応システムと見なしてシミュレーションする場合に，シミュレーション中に細胞が生長する場合や，細胞分裂が起こる場合に対応した実装がなされている．
この機能は以下の項目が前提として仮定されている．

(注 1)
(注 2)

http://www.nrcam.uchc.edu/
2004 年現在バージョン 1.02．GNU GPL に基づいて公開 http://poczta.ibb.waw.pl/stocks/17

2. 背景2.2. 主な細胞シミュレータ– 細胞は一世代で体積が 2 倍になる
– 細胞分裂が発生すると体積は初期値に戻る

これは，Gillespie のアルゴリズムにおいてシステムの体積 V の代わりにシミュレーショ
ン時間 t を変数とする関数 V(t) を使うことで実現している．

V(t) は式 2.32 で表される式である．T は細胞の一世代の時間である．式 2.32 により，細胞の体積が細胞の一世代の時間で線形に増大していく反応を表現することができる．

V(t) = 1 +tT(2.32)

具体的なシミュレーション方法は，毎サイクルの時刻 t から V(t) を再計算し，cµ の値をV(t) の値を使って更新していくことで実現する．t = T すなわち V(t) = 2 になったとき，細胞の一世代が終わったことになり，細胞分裂が発生する．
細胞分裂が発生すると，体積および cµ の値が初期値に戻り，そのときに存在する各分子の数が 2 で除算される．DNA エレメントをモデル化する分子数は除算の前に 2 倍しておき，分裂後も同じ数を保持するようになっている．
このような線形の体積変化で細胞生長，細胞分裂を表現する方法は，原核生物のシステムを対象としたシミュレーションでは妥当な結果が得られている．

• ランダムプール
任意の数の分子が存在する “ランダムプール” を用意しておき，他から小さな影響を大量に受けて数が変動する分子がある場合に，分子の数を変動させる機能が追加されている．
STOCKS では，RNA ポリメラーゼのモデリングにこのランダムプールの機能を使用している．
RNA ポリメラーゼのランダムプールは以下のような方法が実装されている．
– ランダムプール内の分子は，分子数がサイクルごとに毎回変更される– 分子数は hµ ，aµ を計算する前に，そのときの分子数の平均と標準偏差を用いて正規分布する乱数を得，その値が設定される– 体積の増加とともにプール内の分子数の平均値は増大するが，濃度は一定のままになる
ランダムプールは RNA ポリメラーゼ以外の分子にも適用でき，また反応分子数の明示的なモデル化が難しい大規模プロセスのバッファにもなると Kierzek は予測している．
ただし，ランダムプール機能は，多くの場合に妥当な結果が得られているとは言え，数学的に正しさが証明されているわけではないので，ランダムプール機能を使用する場合には注意が必要になる．

計算コストへの対策
STOCKS は計算コストが大きく，実行時間は数時間から数日に及ぶため，UNIX オペレーティ
ングシステムのバックグラウンドジョブとして実行することに向いた設計が必要になる．
そのための一つの方法として，入力ファイル，制御ファイル，出力・ログファイルの 3 種類のテキストファイルでシミュレーションを定義している．
18

2. 背景2.2. 主な細胞シミュレータ• 入力ファイル
入力ファイルにはシステムの仕様が格納される．反応物質の構成と初期値，および確率モ
デルの速度定数が記述されている．同時に，分裂の際のルールも記述されている．
• 制御ファイル
制御ファイルにはどの反応物質の数をモニタするのかが記述される．プログラムは制御
ファイルに記述されたシミュレーション時間間隔で反応物質の数を記録する．モンテカルロ
反復の回数，乱数の種も制御ファイルで設定する．
• 出力・ログファイル
出力・ログファイルには出力とログが出力されるディレクトリパスを記述する– 出力ディレクトリに書き出されるファイルは軌跡 (Trajectry) ファイルと呼ばれる．軌
跡ファイルは 2 カラムのテキストファイルで，制御ファイルに記述された時間間隔で時間の関数として分子の数を記録している．Gnuplot のような作図ソフトウェアで使用しやすい形になっている– ジョブが中断された場合に復帰できるように，システム内の全ての反応分子数を格納
している

乱数生成器

STOCKS では Marsaglia と Zaman の 2144 サイクルの乱数発生器を使っている．

ユーティリティプログラム

STOCKS にはデータの解析を支援する 4 つのユーティリティプログラムが提供されている．
• 軌道平均算出プログラム
モンテカルロ反復により，同じ分子をモニタしているシミュレーションの軌跡ファイルが複数作成される．これらについて，軌道平均算出プログラムは指定した時間間隔で分子数の平均・標準偏差を計算する．
• 分子数加減算プログラム
2 つの軌道ファイルから，ある時間における分子数を加減算するプログラムである．いく
つかの分子数の合計を知りたい場合に使用する．
• 分子数平均算出プログラム
軌道ファイルの所望の部分について，分子数の平均を計算するプログラムである．分子数が定常レベルに達しているときに使用する．
• 軌道-1 次関数比較プログラム
軌道ファイルの分子数が描く形について，適当な部分に一次関数を合わせることによって，軌道の挙動を確認するプログラムである．分子数が一定速度で増加または減少している
ときに使用する．

19

2. 背景2.2. 主な細胞シミュレータ

STOCKS の計算時間
Kierzek は STOCKS の提案に際し，2 つの反応について例示している．Kierzek が言及している
ように，計算の目的は，反応物質の詳細なモデル構築ではなく，ソフトウェアの性能評価にある．
• 原核生物の転写と翻訳の開始頻度と遺伝子発現の確率変動提案した大腸菌 (E.coli) のタンパク質合成速度レベルと，LacZ 遺伝子の mRNA レベルに関する実験データから Kierzek らが提案した運動モデルをシミュレーションするものである．
シミュレーションは 100 個の独立した軌道について，細胞 10 世代の期間 (2100 秒間) の計算を行い，100 回のモンテカルロ反復の結果を蓄積し，タンパク質の時間変動をプロットした．PentiumIII 800MHz プロセッサでこのシミュレーションに約 22 時間掛かっている．
• E.coli の LacZ，LacY の遺伝子表現およびそのタンパク質の酵素/転送活性STOCKS が酵素反応と少数の反応分子を含むシステムに適用できることをテストするものである．細胞 1 世代 (2100 秒) の 1 つの軌道の計算に，PentiumIII 800MHz プロセッサを用いて約 2.5 時間掛かる．また，単一の軌道について細胞 10 世代のシミュレーションに約
90 時間掛かる．

2.2.3

STOCHSIM

STOCHSIM[10] は Gillesie のアルゴリズムの弱点を克服するものとして Shimizu が開発している確率モデルシミュレータである．
Gillespie のアルゴリズムの弱点克服
従来の解析モデルと同様，Gillespie のアルゴリズムはシステム内の分子の種類と数のみを扱っており，分子の偏在については考慮していない．Gillespie のアルゴリズムが適用できるためには，分子が一様分布していることが前提になっている．これは，異化活性化や親和性による結合など，
多重変形するタンパク質がシステム内に存在する場合に問題になる．分子のそれぞれの状態は異なる特性を持っている場合があり，モデルを詳細化していくと，反応数が数百万に上ることがあるかもしれない．解析モデルや Gillespie のアルゴリズムは反応数に比例してシミュレーション時間が増加するため，これらの機能の実装が困難になっている．
STOCHSIM では，同じ分子でも 3 次元構造や性質が異なるものが存在するシステムのシミュ
レーションへの対応がなされている．

STOCHSIM の特徴• 多状態分子の定義
STOCHSIM は基本的な単分子反応・2 分子反応でシステムの挙動をシミュレートする．シ
ステム内のそれぞれの分子はオブジェクトで定義されている．このオ ブジェクトは “多状態分子 (Multistate molecules)” と呼ばれ，内部に状態を保存しておく分子表現が可能になっている．オブジェクト内部にある状態を示す変数は，状態によって反応発生確率が異なる場合に対応し，反応がもたらす，共有結合修飾や立体構造変化のような分子相互の反応を表現するために使われる．
20

2. 背景2.2. 主な細胞シミュレータ• 反応確率と計算時間の定義
STOCHSIM ではシミュレーション時間を確率モデルで計算せず，ある一定の微小時間間隔ごとに計算を繰り返す，という方法を取っている．シミュレーションの初期化の段階で，全ての反応についての反応確率がユーザ定義による速度定数から計算され，ルックアップ
テーブルに格納される．シミュレーションの 1 サイクルで進行する時間は，そのサイクルで起こりうる反応について，反応速度を保証できる長さに調節される．
この基本設定が行われた後は，単分子・2 分子反応がランダムに選択された分子で発生する，という単純な計算の繰り返しでシミュレーションが進展していく．シミュレーションされるシステムには疑似分子と呼ばれる擬似的な分子が定義されている．疑似分子は STOCHSIM
でのみ定義される仮想的な分子で，実際に存在するわけではない．STOCHSIM は単分子反応と 2 分子反応の反応を簡潔に記述するために，選択した 2 つの反応分子のうち，片方が疑似分子だった場合を単分子反応であると定義している．疑似分子の数は単分子反応の発生頻度を示すものである．

• 実装
STOCHSIM は標準 C++で記述され，Perl/Tk で記述された GUI が提供されている．オリジナルの STOCHSIM シミュレータは Microsoft Windows 用のシミュレータとして開発され，
MS/MFC 製品をベースとした GUI が提供されている．UNIX,Linux，MacOS のような他のOS で STOCHSIM を実行するために，Perl で記述された Tk インタフェースが提供されている．
出力はユーザ定義の時間間隔での分子数や分子の状態であるが，空間の分子の状態や分布を表すグラフィカルスナップショットなども出力することが可能である．

計算アルゴリズム
STOCHSIM のシミュレーションは以下のような流れで実行される．時間は非連続の微小時間に量子化される．それぞれの微小時間ごとに，システム内に存在する分子のうち 1 種類のオブジェ
クトがランダムに選択される (ここで選択する分子は疑似分子ではない)．次に，疑似分子を含めたシステム内に存在する分子の中から 1 つのオブジェクトを選択する．もし，2 つの分子が選択されたならば，対応する分子対で発生する 2 分子反応が発生する．もし 1 つの分子と 1 つの疑似分子が選択されたならば，最初の分子の単分子反応が発生する．
発生した反応に対応する反応確率がルックアップテーブルから読み出される．それから乱数を生成し，反応確率と比較することで反応の発生を判定する．
この結果，反応が発生したならば，分子数の増減，分子の状態の変化など，発生した反応に対応するシステムの更新が行われる．
分子が選択された後に反応が発生する確率は以下のように計算される．
• もし最初に分子 A が選ばれ，次に疑似分子が選ばれた場合，A が単分子反応する確率は，P=

k1 × n × (n + n0 ) × ∆tn0(2.33)• もし 2 つの分子が選択された場合，それらが反応する確率はP=

k2 × n × (n + n0 ) × ∆T
2 × NA × V
21(2.34)

2. 背景2.2. 主な細胞シミュレータ

記号の意味は以下の通り．

• n システム内の分子の数• n0 システム内の疑似分子の数• k1 単分子反応速度定数 (s−1 )• k2 二分子反応速度定数 (M −1 s−1 )• ∆t タイムスケール (分割した微小時間)(s)• NA アボガドロ定数• V システム体積 (l)疑似分子の数は，単分子反応のうち最も高速な反応の反応確率を，2 分子反応のうち最も高速な反応の反応確率に可能な限り近づけた値が設定される (式 2.35)．式 2.35 の [x] は， x に最も近い自然数に近似することを示している．

"

k1,maxn0 = 2 × NA × V ×k2,max#(2.35)

多状態分子 — multiple state molecule
STOCHSIM でそれぞれの分子オブジェクトのインスタンスに内部状態を保持している．酵素の活性やタンパク質の信号伝達は，共有結合修飾やリガンドやサブユニットの結合，タンパク質の構造変化のような多数の要素で制御される．これらは分子オブジェクトのインスタンスの内部状態でモデル化されている．内部状態 1 ビットのフラグ列で管理されており，分子の状態を表現す
る．たとえば，外部のリガンドが膜輸送受容体へ結合するかどうかを決定するフラグなどがある．

開発状況
STOCHSIM は 2004 年段階のバーションは 1.4 である．微小時間での計算，多状態分子など，
Gillespie のアルゴリズムに比べて現実的なモデル，たとえば熱力学的により現実的なモデルを扱うことができるが，計算時間やメモリ要求はより大きくなっている．開発当初のバージョン 1.0 では均一に混合された溶液に関する反応系を取り扱うことができたが，バージョン 1.2 で空間を 2 次
元格子に分割し隣接した場所にある分子との相互反応をシミュレーションすることができるようになっている．これは，細胞内の膜上に密集した受容体などの研究に使用できる．空間の表現方
法は現在も拡張が続けられている．

22

第3章

ReCSiP の構成
本章では、細胞シミュレーションプラットフォーム ReCSiP(Reconfigurable Cell Simulation Platform) の設計方針について述べ、続いてその支援ツール群を含めたシステム構成について述べる。

3.1 ReCSiP の設計方針
ReCSiP は，FPGA を搭載した 1 枚のボードを一般的な PC/WS 環境に付加することで，現在はクラスタなどによって実現されている強力な計算能力を実現することを目標としており，全体の設計においては，特に以下の 4 点に重点がおかれている．
柔軟性 高い性能とともに，さまざまなアルゴリズムが実装できる柔軟性を備えること互換性 ホストとなるシステムをなるべく選ばないようなアーテキクチャであること
容易性 ハードウェアや低レベルプログラミングに関する知識のないユーザでも簡単に利用できる
ようなプログラミングインターフェイスを備えること
可搬性 アクセラレータのハードウェアがない環境でも，ソフトウェアによって代行処理が可能で
あること

3.2

シミュレーション高速化のためのストラテジ

現在，バイオインフォマティクスの分野では BLAST[11] などを用いた遺伝子配列の解析が盛んであり，さまざまな研究機関で大規模な PC クラスタが稼働している．これは，遺伝子配列の解析という問題は主に sequence matching であり，容易に並列化し，かつ高い並列度で実行することが可能であるためである．
Sequence matching を行うための計算機としては，FPGA を用いて Sun のワークステーションに
attach する形のアクセラレータである Splash-2[12] が過去に大きな成果を挙げており，分子動力学の分野では，多体問題を解くために開発された専用計算機 GRAPE[13] が，蛋白質の構造解析などで成果を発揮しつつある．
代謝系シミュレーションにおいては，個々の化学反応は互いに独立な現象であるため，問題そのものに並列性は充分にあるものの，• 複数の反応に関与する物質が多く，その依存関係が極めて複雑であること
23

3. ReCSiP の構成

3.3. ReCSiP の構成• 物質の移動するモデルを扱う場合，計算の途中でシミュレーションする系全体に関するデー
タが必要となることがあるといった問題があり，クラスタなどを用いた場合には通信がボトルネックになる可能性を含んでいる．この場合に，FPGA などに複数の演算器を載せて並列処理を行えば，通信オーバーヘッドを発生させることなく効率のよい処理を実現することができると考えられる．
しかしながら，FPGA を用いたシステムをシミュレータから利用する場合，ハードウェアに密接したプログラムを書く必要があり，低レベルプログラミングに不慣れなユーザにとっては困難である．また，FPGA は並列処理による高速な演算処理を可能とするが，複雑な処理を行うことは困難であるため，• ユーザプログラムの開発容易性• ホストのプログラムと FPGA の処理の適切な切り分けの実現などを考慮した，ホスト–アクセラレータ協調処理環境を構築する必要がある．ReCSiP は，図 3.1に示すように，ソフトウェア・ハードウェアの両面でシミュレーションシステム開発に必要なリソースを提供しており，これらを利用することで FPGA を利用したアプリケーションを容易に開発したり，既存のアプリケーションを移植することで高い性能を得ることが可能になる．
このように ReCSiP は細胞シミュレーション (代謝シミュレーション) を主眼において設計されているが，実際にはシミュレータなどの特定のアプリケーションにターゲットを絞ったものではなく，シミュレーションの他にも顕微鏡などの測定機から送られてくる大量のデータをリアルタイムに処理 [14] することができる予定であったり，あるいは一度ディスク上に保存されたデータを高速に処理するシステムを ReCSiP の上に構築することも可能である．

3.3 ReCSiP の構成
3.3.1

全体像

ReCSiP のシステム構成を図 3.1 に示す．
ReCSiP が基本コンポーネントとして提供するのは，• ハードウェア資源– ReCSiP Board: FPGA を搭載したボード本体– ReCSiP Verilog Library: FPGA に演算用のモジュールを実装する際に用いる基本的なモ
ジュール類• ソフトウェア資源– ReCSiP Driver: ReCSiP ボードのデバイスドライバ– ReCSiP API: ReCSiP ボードを制御するための基本的なインタフェース– ReCSiP 拡張 API 生成ツール: FPGA に搭載したモジュールを制御するコードを書くための補助ツール

24

3. ReCSiP の構成

3.3. ReCSiP の構成

Controller
SolverPathway
Map

SolverPathway
Map

SSRAM

SolverPathway
Map

SSRAM

SolverPathway
Map

SSRAM

SolverPathway
Map

User Program(Simulator)

Interface Code

ReCSiP Driver
Operating System

SDRAM

SSRAM

FPGA

Host Program

ReCSiP Board

図 3.1: ReCSiP の構成
SSRAM

Virtex-II
XC2V6000

SSRAM SDRAM

図 3.2: ReCSiP ボード
といったものである．
ユーザは，FPGA 上に目的とする演算モジュール (solver) を実装，あるいは既存の solver を利用
し，ReCSiP の提供するソフトウェア資源を利用することで，対象アプリケーションの高速化を実現する．
ReCSiP の提供する個々の基本コンポーネントについては以下で述べる．

3.3.2

ReCSiP ボード

ReCSiP ボード (図 3.2) は，FPGA による高速な演算を実現するために，FPGA に PCI インターフェイス，メモリ等を接続したボードである．ボードの構成を図 3.3 に，ボードに搭載されている主要なコンポーネントを表 3.1 に示す．

25

3. ReCSiP の構成

3.3. ReCSiP の構成

Micron 256Mb SDRAM
MT48LC16M16A2TG-7E x2Expansion
Connector

XilinxMicron 18Mb SSRAM

Virtex-II

MT58L1MY18DT-7.5 x2

XC2V6000-4BF957C

Local Clock

Local Clock
Generator

64bit Local Bus

USB Port

MultiLINX
Cable

Configuration via
USB/MultiLINX

Configuration
Control

1.5V
DC Power
SupplyQuickLogic
Configuration via
PCI

QuickPCI

DC 5V

3.3V

QL5064-PB456C-66B

ReCSiP Board

64bit/66MHz PCI Bus

図 3.3: ReCSiP ボードの構成
表 3.1: ReCSiP ボードの構成部品

Core FPGAPCI Interface
Memory (SRAM)
Memory (DRAM)

メーカー

シリーズ

型番

XilinxQuickLogic
MicronMicronVirtex-II
QuickPCI
SyncBurst SRAM
Synchronous DRAM

(XC2V6000-4BF957C)(QL5064-PB456C-66B)(MT58L1MY18DT-7.5)(MT48LC16M16A2TG-7E)

Core FPGA
コアとなる FPGA には，Xilinx 社の Virtex-II シリーズの製品である，XC2V6000-4BF957C を採用している．ボード上のほぼすべてのコンポーネントは Virtex-II に接続されており，Virtex-II は計算処理を行うとともに，ボード上のコンポーネントの制御を行う．
ボードに搭載されている Virtex-II は，6M ゲート相当で 8,448 の CLB (Configurable Logic Blocks)と 2,592kbits の内部 RAM，684 のユーザ I/O ピンを持っており，大規模な演算ロジックやコント
ローラを構成することが可能である．

PCI インターフェイス
PCI インターフェイスは，QuickLogic 社の FPGA である QL5064-PB456C-66B が用いられている．この LSI は，• PCI インターフェイスロジックと，ユーザがプログラム可能で高速な動作が可能なアンチヒューズ型 FPGA を 1 チップに統合• PCI インターフェイスはマスタ DMA 転送が可能で，32bit/64bit, 33MHz/66MHz, 3.3V/5V の各規格に対応• PCI インターフェイスと FPGA 部分のクロックは分離可能

26

3. ReCSiP の構成

3.3. ReCSiP の構成

表 3.2: ローカルバスの信号線

SYSCLK

inなし

ローカルバスクロック．基板上のオシレータから QL5064と Virtex-II のクロックピンに供給される

L RSTXout

負

システムリセット信号．PCI バスのリセット信号に連動して動作

L DRQXL DAKXL ADSX

out
in
i/o

負

ターゲット転送要求信号

負

L WEX

i/o

なし/負

L BSTMXL RDYOX
L RDYIX
L AD[63:0]

i/o
out
in
i/o

負

L BEX[7:0]

i/o

なし/負

L INTX
L S64X
L DBSYX

in
in
out

負

L DRQX に対する acknowledge
転送サイクルの開始とアドレスの転送を示すアドレススト
ローブ．ターゲットアクセス時は出力 (QL5064→Virtex-II)
で，マスタアクセス時は入力
Write 動作を示す信号で，ターゲットアクセス時は出力，マ
スタアクセス時は入力
バースト転送を示す
QL5064 がデータ転送可能であることを示す
Virtex-II がデータ転送可能であることを示す
アドレス/データ信号線．ターゲットアクセス時は最初に
L ADSX と同時にアドレスが転送され，続いてデータが転送される．マスタアクセス時は L ADSX と同時にアドレスが，続いて転送バイト数が転送され，続いてデータが転送される
バイトイネーブル信号．ターゲットアクセス時には出力，マスタアクセス時には入力
割り込み信号
コントロールレジスタ選択信号
PCI bus の状態を示す．Low の場合には PCI ライトマスタ
アクセスは行えない

負

負
負
正

負
負

といった特徴を備えており，ReCSiP ボードのカードエッジにある PCI バスコネクタに接続されている．このコネクタは 64bit/66MHz，3.3/5V 対応になっており，66MHz のバスに接続した場合にはボード上のジャンパで 33MHz 動作を強制することもできる．

ローカルバス

QL5064 の FPGA 部分にはローカルバスを制御するためのロジックが書き込まれ，PCI バスとローカルバスの間のブリッジとして動作する．ローカルバスは表 3.2 に示す信号線で QL5064 とVirtex-II の間を接続しており，PCI バスとはクロックが切り離されている．ローカルバスのクロッ
クはボード上の水晶発振器から供給されており，任意の周波数のものに交換することが可能である．
ローカルバスの target read/write 動作におけるタイミングを図 3.4 に示す．なお，ローカルバスのアドレス空間は 27bit である．

27

3. ReCSiP の構成

3.3. ReCSiP の構成

SYSCLK
L_DRQX
L_DRQX
L_S64X
L_ADSX
L_AD[63:0]

A

D0

D1

D2

D3

A

D0

D1

D2

D3

L_WEX
L_BSTMX
L_RDYOX
L_RDYIX

Target Read

Target Write

図 3.4: ローカルバスの read/write 動作

RAM
メモリは 18Mb の同期 SRAM である，Micron 社の MT58L1MY18DT-7.5 を 2 つずつ 4 セット(注 1) と，256Mb の同期 DRAM である，Micron 社の MT48LC16M16A2TG-7E を 2 つ実装している．
各メモリのデータ幅は 16bit であり，2 つを 1 セットで並列に Virtex-II に接続することで 32bitのデータ幅にしている．SRAM の各セットの容量は 4MB，計 16MB であり，DRAM は 64MB である．
各メモリへ供給するクロックは，Virtex-II 上の DLL(Delay Locked Loop)) により，水晶発振器から供給されるクロックを分周または逓倍して供給することができる．クロックスキューを低減するため，ボード上のクロック配線から DLL にフィードバックを行うことも可能である．
SRAM は制御が簡単で柔軟なアクセスパターンに対応できるため，シミュレーションの中間デー
タなどを格納しておき，SDRAM は計算結果などを保存するために用いられる．

電源

PCI バスから供給される電源は最大 25W と定められている [15] が，SRAM の最大消費電力は 1
チップあたり 1.7W で，8 チップでは 13.5W となる．これに FPGA の消費電力が加わると電圧低下などの危険があるため，ボード上にレギュレータを搭載し，電源を供給している．
FPGA が用いる 1.5V の電源を供給するのは Linear Technology 社の LT1585A-1.5 で，最大出力は 5A である．メモリ等の 3.3V の電源は LT1584-3.3 が供給し，最大出力は 7A である．それぞれのレギュレータへの電力は，PC に内蔵される一般的な ATA 電源コネクタから供給される．

3.3.3

Verilog ライブラリ

ReCSiP では，ReCSiP ボード上の Virtex-II に回路を実装する際に必要となる，いくつかのモジュールを標準で提供する．ここでは，これらの実装について述べる．
(注 1)

最初の 1 枚の基板に限り，8Mb の Micron MT58L512L18DS-7.5 であり，容量は半分となる．

28

3. ReCSiP の構成

3.3. ReCSiP の構成

表 3.3: 標準ローカルバスインターフェイスのレジスタ割当
アドレス

方向

内容

0x0000 0000
0x0000 0008
0x0000 0010
0x0000 0018
0x0000 0020
0x0000 0028

W
W
W
W
R
W

マスタ動作時のボード側物理アドレス
マスタ動作時のホスト側物理アドレス
write 転送長さ (bytes)
read 転送長さ (bytes)
割り込み制御レジスタ
ユーザ回路リセット

標準ローカルバスインターフェイス

QL5064 とのハンドシェイクを行い，ホストとのデータ転送を行うのがこのインターフェイス
モジュールである．標準で実装されている制御レジスタを表 3.3 に示す．
マスタ動作時のボード側／ホスト側物理アドレスを設定し，write 転送長さ，あるいは read 転
送長さのレジスタに転送長を設定すると，転送長設定レジスタへの書き込みがトリガとなって，
QL5064 経由でマスタモード転送を開始する．
割り込み制御レジスタは，ReCSiP ボードが PCI バスの INTA#をアサートして割り込みを発生する場合に用いる．Virtex-II 側では，INTA#をアサートする前にこのレジスタを non-zero な値にしておき，ホストからレジスタに対するリードアクセスが発生するまで保持する．ReCSiP ボードがINTA#をアサートしていない間は，このレジスタの値は 0 に保持される．
ユーザ回路リセットレジスタへの書き込みアクセスにより，ユーザ回路リセット信号が発生するようになっており，Virtex-II 内のすべてのモジュールは，このユーザ回路リセット信号を参照
しなければならない．
．リセット動作時はユーザロジックとローカルバスインターフェイスからク
ロック制御系の DLL にリセット信号が出たのちに，DLL による位相制御がロックしたことを確認してから，ユーザロジックへのリセットが解除される．これにより，DLL を利用するユーザは，
DLL のステータスを気にせずにロジックを設計することが可能になる．
また，空いているアドレス空間にユーザが任意のレジスタを実装することも容易である．空き
アドレス空間は，これらの標準レジスタの後から，表 3.4 に示すメモリ類の前までであるが，システムの拡張を考慮して 0x50 より下のアドレス空間は予約されている．
BlockRAM コントローラ
多くのユーザは 32bit 幅のメモリを必要とするが，ローカルバスの幅は 64bit である．したがって，メモリに書き込む際には 64bit から 32bit への変換を要する．Virtex-II の内部 RAM (BlockRAM) を用いる場合に，この変換を行うモジュールが BlockRAM コントローラである．
このコントローラには，1 つの 64bit インターフェイスと，2 つの 32bit インターフェイスがあり，8bit 幅の dual port な BlockRAM が 4 つ接続される．BlockRAM 自体には byte write enable の機構がないため，8bit 幅の RAM を複数用いることで byte write enable に対応している．64bit イ
ンターフェイスは，ローカルバスインターフェイスに接続され，32bit インターフェイスはユーザ
ロジックに接続される．

29

3. ReCSiP の構成

3.3. ReCSiP の構成

A

D

Q

W_ BE

64bit Interface

ENA_

A

Control Signal

D

Q

W_

32bit Interface A

A

D

Q

W_

32bit Interface B

BlockRAM Controller

CLK
CLK2

Port A Interface (8bit x 4=32bit)

A D Q
Port AA D Q
Port APort B Interface (8bit x 4=32bit)

A D Q
A D Q
Port APort B

A D Q

A D Q
Port B

Port AA D Q
Port B

BlockRAM
BlockRAM(8bit x 1024 word)
BlockRAM(8bit x 1024 word)
BlockRAM(8bit x 1024 word)

A D Q
Port B(8bit x 1024 word)

図 3.5: BlockRAM コントローラ
表 3.4: 標準ローカルバスインターフェイスのメモリアドレス割当
アドレス

内容

0x0200 0000
0x0280 0000
0x0300 0000
0x0380 0000
0x0400 0000

SRAM0
SRAM1
SRAM2
SRAM3
DRAM

ふたつの 32bit インターフェイスは，BlockRAM と同じタイミングで使用することができる．
32bit インターフェイスの片方 (A) は 64bit インターフェイスと排他利用であり，もう一方 (B) は常にアクセスが可能である．BlockRAM 同様に，32bit インターフェイス A と 32 ビットインター
フェイス B，または 64bit インターフェイスと 32bit インターフェイス B は dual port RAM として動作し，同時アクセスが可能である．
64bit インターフェイスからの read/write を 32bit に変換するため，64bit インターフェイスがenable されている間，BlockRAM のポート A にはシステムクロックの 2 倍の周波数のクロックが供給され，BlockRAM コントローラがデータ幅変換を行う．

SSRAM/SDRAM コントローラ
SSRAM/SDRAM のローカルバス上のアドレスを表 3.4 に示す．これらのコントローラは現在デバッグ中である．

浮動小数点演算器

ReCSiP では，Verilog-HDL で記述された浮動小数点演算器として，単精度の加減算器と乗算器を用意している．表 3.5 に示すように比較的コンパクトかつ高速に作られており，Virtex-II 用の30

3. ReCSiP の構成

3.3. ReCSiP の構成

表 3.5: 浮動小数点演算器の仕様
加減算器
パイプライン
動作周波数
使用スライス数
フォーマット

乗算器

5段
8段
83MHz
84MHz
384 (1.13%)1125 (3.33%)
IEEE-754 準拠 単精度浮動小数点

表 3.6: 市販 IP コアの浮動小数点演算器の仕様
加減算器
パイプライン
動作周波数
使用スライス数
フォーマット

乗算器

4段
7段
58MHz
74MHz
605 (1.79%)
742 (2.20%)
IEEE-754 準拠 単精度浮動小数点

IP コアとして Digital Core Design 社が提供しているもの (表 3.6) と比較しても遜色ないものであることがわかる．

3.3.4

ソフトウェアReCSiP が提供するソフトウェアは，
1. 基本的なボードとの通信機能を提供するドライバ
2. ユーザプログラムからドライバを利用するための API
3. FPGA に実装した回路とのインターフェイスを作成するための拡張 APIから構成される．
ドライバは ReCSiP ボードに割り当てられた 128MB のアドレス空間に対して読み書きを行うためのインターフェイスを提供しており，次に示す API を利用してアクセスすることが可能である．

a[31:0]addsubfloataddsub32b[31:0] clk

a[31:0]b[31:0]floatmult32

out[31:0]

out[31:0](a) 加減算器(b) 乗算器

図 3.6: ReCSiP の浮動小数点演算器

31

clk

3. ReCSiP の構成

3.3. ReCSiP の構成

ReCSiP API• デバイスの使用宣言／開放– int

rc_open();// device open

– void rc_close(); // device close
ReCSiP ボードを使用する場合には rc open() によって使用を宣言する．複数のプロセスによるボードの利用は許されず，既に使用しているプロセスが存在している場合にはエラーが返される．ReCSiP ボードは rc close() によって解放される．
• PIO による 32bit 転送– long rc_single_read(unsigned long device_addr);– void rc_single_write(unsigned long device_addr, long data);
ReCSiP ボードのローカルアドレスを指定して 32bit の read/write 操作を行う．
• PIO によるブロック転送– void rc_block_read(void *user_addr, unsigned long device_addr,
size_t size);– void rc_block_write(unsigned long device_addr, void *user_addr,
size_t size);
ホスト側のポインタと ReCSiP ボードのローカルアドレスを指定してその間でブロック転送を行う．write 操作時には PCI ブリッジによって自動的にバーストサイクルが発生することがある．

• バスマスタ転送– void rc_master_read(void *user_addr, unsigned long device_addr,
size_t size);– void rc_master_write(unsigned long device_addr, void *user_addr,
size_t size);
使い方は PIO によるブロック転送と同じであるが，これらの関数を用いた場合には ReCSiP
ボード上の QL5064 の DMA によってデータが転送される．

• メモリ空間へのマップ– void *rc_mmap(unsigned long device_addr, size_t size, int prot);– int

rc_munmap(void *ptr, size_t size);rc_mmap は，ReCSiP ボードのローカルアドレスとサイズを指定し，ユーザプログラムのメ
モリ空間にマップし，rc_munmap は，これを解除する．
• その他
32

3. ReCSiP の構成

3.3. ReCSiP の構成– void rc_wait();– void rc_load(char *file);– void rc_reset();これらは ReCSiP ボードからの割り込み待ち，Virtex-II のコンフィギュレーション，およびVirtex-II 上のユーザ回路のリセット操作を行うための関数である．

拡張 API

API を直接使用する場合，ユーザは Virtex-II 上に存在する資源が，ReCSiP ボード上のローカル
アドレス空間にどのようにマッピングされているかを常に意識しなければならない．
拡張 API は，Virtex-II 上に構成されたレジスタやメモリのアドレスマップから，それらをポインタとして名前ベースの操作で扱えるようなヘッダファイルを生成するプログラムによって提供され，ユーザはこのヘッダファイルを呼び出すことで mmap() などの操作を意識せずに，ボード上のリソースにアクセスすることが可能である．

33

第4章

設計と実装
4.1

変数の定義

Gillespie のシミュレーションアルゴリズムを設計・実装するにあたり，以下のように変数を定
義する．
Xν
hµ
cν
aµ
a0
τ
µ

4.2

:
:
:
:
:
:
:

システム内に存在する S ν 分子の数
反応 Rµ にかかわる分子の組み合わせ数
システム温度と分子の物理的性質に依存する確率反応速度定数
各反応の発生確率を示す指標
発生反応と時間経過計算の基礎的な値
次の反応が発生するまでの時間
次に発生する反応番号

シミュレーションの流れ

Gillespie のアルゴリズムを用いて，Lotka 反応モデルのシミュレーションを行うモジュールを
ReCSiP 上に実装した．
Lotka 反応モデルは，初期値 Xν (ν = 1, · · · , 4) と cν (ν = 0, · · · , 3) を設定し，(0, 1) の一様乱数
r1 ，r2 を毎回生成した上で 1 反応ごとに以下のステップの計算を行うことでシミュレーションが進行する．
Step 1. hν の計算
シミュレーションのサイクルの実行時に，その時点の分子の数をもとに，hν を計算する．

h1 = X1 · X2

(4.1a)

h2 = X1 · X3

(4.1b)

h3 = X3(4.1c)

h4 = X4(4.1d)

34

4. 設計と実装

4.2. シミュレーションの流れ

Step 2. aν および a0 の計算
Step 1 で求められた hν と初期設定されている cν を掛け，aν を計算する．さらに，aν を全て加算し，a0 を算出する．
a1 = h1 · c1

(4.2a)

a2 = h2 · c2

(4.2b)

a3 = h3 · c3

(4.2c)

a4 = h4 · c4
4
X
a0 =
ai

(4.2d)(4.2e)

i=1

Step 3. τ の計算
a0 と (0, 1) の一様乱数 r1 から，次の反応が起こる時間 τ を計算する．
1 1
τ=
ln
a0
r1

(4.3)

Step 4. µ の計算
aν と a0 を比較し，“次に発生する反応” を示す µ を求める．
µ−1Xaν < r2 a0 6

µ=1

µ
X

aν

(4.4a)

ν=1

Step 5. 反応結果のフィードバック
求められた µ の値から，対応する反応を実行し，分子数を増減する．Lotka モデルの場合は，式4.5 で表される式が反応で起こる分子数の変化を示している．
if µ = 0
 X2 + 1(4.5a)
(X − 1)1


 X2 − 1 X3 + 1

if µ = 1

 X3 − 1 X4 + 1

if µ = 2 X2 − 1 X4 + 1

if µ = 335(4.5b)(4.5c)(4.5d)

4. 設計と実装

4.3. 実装

Step 1 から 5 の一連の処理がシミュレーションの 1 サイクルになる．実装における Lotka モジュールにおいて，この処理を実行する．

4.3

実装

4.3.1

実装モジュール

実装した Gillespie のシミュレーションアルゴリズムの実装にあたり，基本的な演算器以外に実装したモジュールの名称と機能を定義する．

Lotka 反応セットモジュール:

Lotka 反応モデルモジュール:

反応実行モジュール:対数テーブル:

変数保存テーブル:

今回の実装の最上位モジュール
Lotka 反応モデルモジュールを 2 つ搭載し，メモリへの入出力を持つ
Lotka 反応モデルシミュレータ
Gillespie のアルゴリズムを用いてシミュレーションを実行する
各種演算器と反応実行モジュールを内部に持つLotka 反応の分子数変更モジュール式 4.5 に従い，入力された反応により分子数の増減を行う
32 ビット幅 15 ビットエントリのデュアルポートテーブル(0, 1) の乱数 r1 より計算された ln (1/r1 ) の値が格納される
32 ビット幅 37 エントリのデュアルポートテーブル
シミュレーション中の分子数を格納する

Gillespie のアルゴリズムを用いて確率モデルシミュレーションを実装するために，必要となる
基本的な演算モジュールは以下のものが挙げられる．
• 32 ビット 整数 加減算器
• 32 ビット 整数 乗算器• 単精度 浮動小数点加減算器• 単精度 浮動小数点乗算器• 整数-浮動小数点実数変換器• 浮動小数点比較器• 一様乱数生成器• 対数格納用テーブル• 変数用テーブル36

4. 設計と実装

4.3.2

4.3. 実装

浮動小数点実数の算術演算器

Gillespie のアルゴリズムでは，(0, 1) の乱数が使用され，対数計算を含めた有効数字の大きな実数が使用できることが望まれる．そのため，浮動小数点実数を扱えることが必要になる．実装す
る浮動小数点実数の形式は，IEEE754 に準拠した単精度浮動小数である．
シミュレーションアルゴリズムに使用する演算器は，ReCSiP が標準的に提供している浮動小数
点加減算器を使用している．
現行の ReCSiP 標準の単精度浮動小数点乗算器は，整数乗算部分が通常の RTL で記述されており，Vertex-II の特有の機能である “18 × 18 ビットの乗算ブロック (embedded multiplier)” を使用していない．そのため，使用面積が大きくなってしまっている．これを解決するために，Vertex-IIの乗算ブロックを使って整数演算器および浮動小数点乗算器の整数乗算器部分を実装した．乗算
ブロックは 64 ビットまでの任意ビットの整数乗算器に使用することが可能である．整数乗算器はXilinx 社提供のツール ISE CORE Generator 6.1.03i を使用して自動的に生成した．
実装した浮動小数点乗算器の性能を表 4.1 に示す．表 4.1 のスライス数，乗算ブロック数は，ReCSiP に搭載されている FPGA である Vertex-II XC2V6000 の全スライスに対する割合を示している．ReCSiP で標準的に提供されていた浮動小数点乗算器と比較して，約 1/7 の面積の面積と高速化を実現できていることがわかる．
表 4.1: 浮動小数点乗算器の比較

ReCSiP 標準
パイプライン
動作周波数
使用スライス数
使用乗算ブロック
フォーマット

乗算ブロック使用

8段
7段
84 MHz110 MHz1125 (3.33%)
159 (0.471%)
0 (0.00%)
4 (2.78%)
IEEE-754 準拠 単精度浮動小数点

組み込み乗算器を使用した浮動小数点乗算器は，図 4.3.2 のような過程を経て処理を行っている．
計算の過程は以下の通り．

• phase 1.
– 乗数 a，被乗数 b の符号ビット (最上位ビット) の排他的論理和を取る– a，b の指数ビット (31〜24 ビット目の 8 ビット) を整数加算する (出力は桁上げ含め 9
ビット)– a，b の仮数ビット (23〜1 ビット目の 23 ビット) を組み込み乗算器を使用して整数乗算
する (4 段パイプライン)• phase 2〜4.
整数乗算の計算を待つ間，符号部はビット数 ×6 段，指数部はビット数 ×5 段の長さのシフ
トレジスタに格納し，次の入力に対応する• phase 5.
37

4. 設計と実装

4.3. 実装asign

exponent
8 bit

matissa
23 bit

1 bit
reg

1 bit
xor

8

integeradder

23

in 8bitout 9bit

1 bit
reg

1 bit
reg

1 bit
reg

1 bit
reg

q
9

9 bit
reg

9

9 bit
reg

9

9 bit
reg

9

9 bit
reg

sign

9
integersub

8exponent
8 bit

in 9bit
8bit
out 9bitbmatissa
23 bitsign

88exponent
8 bit

integer multiplier(4 stages pipelined)matissa
23 bit23

48

in 48bitout 23bit
8bit

in 24 bitout 48 bit

1

2

normalize

3

4

23

5

23 bit
reg

6

23

7

図 4.1: 組み込み乗算器を使用した浮動小数点乗算器– 仮数部の整数乗算の結果を正規化し，1.x × 2y の形にし，x の上位 23 ビットを出力する– y の値が 1 であったら 127 を，0 ならば 126 を次の指数部の整数減算器の入力とする• phase 6.
指数部の値から，正規化モジュール (normalize) から入力された値を減算し，乗算結果の指数レジスタ (8 ビット) に格納する• phase 7.
符号ビット，指数ビット，仮数ビットを乗算結果として出力する
phase 5. において，127 または 126 を指数ビットから減算する理由は，浮動小数点の指数部はイ
クセス表現を取っており，+127 を 0 とみなしているため，phase 1. での指数ビットで整数加算された余分なイクセス表現部を取り除く必要があるからである．

4.3.3

一様乱数生成器

一様乱数生成器として，32 ビットの線形フィードバックシフトレジスタ (LFSR:Linear Feedback
Shift Register) を用いた M 系列乱数発生器を実装した．
LFSR は排他的論理和で帰還をかけたシフトレジスタによるカウンタである．構成するビット数を n とすると，2n − 1 までの疑似乱数を生成することができる．確率モデルシミュレータでは，
15 ビットと 23 ビットの一様乱数が必要になるため，32 ビットの LFSR を 2 つ生成し，それぞれ

38

4. 設計と実装

4.3. 実装

上位の 15 ビットを取るモジュール rgen1 と，23 ビットを取るモジュール rgen2 を実装した．この動作を図 4.2 に示す．

1 bit shift / clock
16

31

0

図 4.2: 32 ビット LFSR による疑似乱数生成モジュール
このモジュールは 1 クロックごとに以下のような動作を行う．

• 32 ビットのレジスタのうち，31 ビット目と 16 ビット目を取りだし，排他的論理和を取る• レジスタ全体を 1 ビット右シフトし，0 ビット目に排他的論理和を取った結果のビットを格納する• 最上位ビットから必要な分のビット列を出力する(0, 1) の浮動小数型実数の乱数を得るためには，23 ビットの疑似乱数 M を生成し，その上位に001111111 という 9 ビットを付加し，32 ビットの浮動小数を生成する．付加した 9 ビットのうち，最上位は符号ビット (0 は正)，続く 8 ビットは浮動小数の指数部 (01111111(2) = 127(10) はイクセス表現で 127 が加算されているため指数部は 0 となる) を意味している．この過程で生成された乱数値 r2 は，23 ビットの疑似乱数を M とすると，
r2 = 20 × 1.M(2) = 1.M(2)(4.6)

である．これは，小数点以下 23 ビットの (1, 2) の一様乱数と見なすことができる．そのため，浮動小数点減算器を使い，r2 から 1.0 を減ずることで，(0, 1) の実数型一様乱数を得る．

4.3.4

対数テーブル対数を使った計算に用いるテーブルは Vertex-II 上に搭載されている BlockRAM(XC2V6000 の使用可能サイズは 288KB) に実装した．このテーブルには (0, 1) の一様乱数 r1 を使用した 1/ ln(1/r1 )の値が格納されているが，便宜上これを “対数テーブル” と呼ぶことにする．BlockRAM はデュアルポートで使用し，確率モデルシミュレータを 2 つ FPGA 上に生成し，1 つのテーブルから同時に 2 つ値を読み出せる構成とした．また，現時点では対数テーブルのデータをシミュレータ側から書き換える必要が無いため，この対数テーブルは読み出し専用で構成している．
テーブルは 32 ビットのデータを 15 ビットエントリ (215 = 32768 個) 格納できるサイズで構成
し，テーブルの中に書き込む値は，0 から 1 の間をエントリの数 (15 ビットエントリなので 0 から1 の間を 215 = 32768 個に分割) に等分し，それぞれの値を r1 とし，その r1 から 1/ ln(1/r1 ) を計算
し，32 ビット浮動小数点表現で格納している．シミュレータでは，LFSR で 15 ビットの乱数を生39

4. 設計と実装

4.3. 実装

成し，その値をアドレスとしてテーブルから 32 ビットのデータを引いてくることで，1/ ln(1/r1 )の浮動小数点型実数を得ることになる．
対数テーブルは，Vertex-II の BlockRAM 上に 32 ビット幅 13 ビットエントリのテーブルを 4 つ構成し，15 ビット乱数の上位 2 ビットから値を取り出すテーブルを選択するモジュールを実装した．テーブルを 4 つに分割した理由は，ISE で BlockRAM を構成する場合に，巨大なテーブルにデータを書き込む場合に非常に長時間かかることと，将来の対数テーブル書き換えに対応するた
めである．

4.3.5

Lotka 反応モデルの実装

Gillespie のアルゴリズムを用いて Lotka 反応モデルのシミュレーションを実行する Lotka 反応
モデルモジュールを実装した．このモジュール 1 サイクルのデータフローは図 4.3 のように進行する．
Lotka 反応モデルモジュールの 1 サイクルは，以下に示す 7 つのフェーズで実行される (phase 0.
は初期化の段階)．各フェーズの (Step X) は，4.2 節での Step に対応している．
phase 0.
シミュレーションの実行前に各種初期化を行う．

• 反応に関わる分子の物理的性質とシステムの温度のみに依存する cµ の初期値を設定する• 乱数種 2 つを乱数生成器用に設定する• 各分子の初期分子数を設定する．パイプライン化した場合，パイプラインの数だけ入力を繰
り返す
phase 1.
Step 1. : シミュレーションサイクル第一段階．サイクル開始時の分子数から，反応ごとに関わる分子の組み合わせを示す値 hν を計算する．
h1 = X1 · X2
h2 = X1 · X3
h3 = X3
h4 = X4• 変数保存テーブルから読み出された分子数から，分子の衝突で反応が発生する可能性がある
分子，または単一の分子で反応が起こる可能性がある分子を選択する• 異種の分子の場合，反応に関わる分子の組み合わせの数を整数乗算により計算し，整数を浮動小数点実数に変換する• 衝突無しで反応する可能性のある分子 (X2 ，X3 ) も分子数を浮動小数実数に変換する

40

4. 設計と実装

4.3. 実装x1
32
32 bit integer
multiplierconstant c
C1
C2
C3
C4

x2

x3

32

32

64

64

integer-FP
transfer

integer-FP
transfer

h1
32

32
32
32
32

phase 132 bit integer
multiplier

32 bit FP
multiplierinteger-FP
transfer

integer-FP
transfer

h2

h3

h4

32

32

32

32 bit FP
multiplier

32 bit FP
multiplierphase 232 bit FP
multiplier

32 bit FP
adder

phase 332 bit FPadder
32
32 bit FP
adder
a0

32

a1

a2

a3

phase 423 bit random
generator

32

constant
1

23

32 bit width
15 bit entry
lntabler2

FP compareter

32 bit FP
multiplierphase 5

3

mu

r1
32

32

32 bit FP
multiplier32

15 bit random
generator
15

32 bit FP
substracter
32

32phase 6

Lotka
phase 7
X1

X2

32

X3
32

32

X4

32

32

BlockRAM
X1

X2

32

X1

X3

32

X2

X4

32

X3

32

X4

1/tau

図 4.3: Lotka 反応モデルの確率モデルシミュレーションの計算過程

41

4. 設計と実装

4.3. 実装

phase 2.
Step 2. 前半 : aν を計算する．
a1 = h1 · c1
a2 = h2 · c2
a3 = h3 · c3
a4 = h4 · c4
• システム温度と分子の物理的性質で決定される固定値 cµ と，phase 1. で計算した分子の組
み合わせ数 hν を乗算し，aν を求める
phase 3.
Step 2. 後半 : phase 2. で求められた aν を全て加算することで，次の反応までの時間 τ，次に起こる反応 µ を計算するための要素となる a0 を計算する．
a0 =

4
X

ai

i=1

• phase 2. の結果を順番に加算する (最後の加算の結果が a0 になる)．a1 から a3 まで順番に加算した結果を phase 6. で使用する．そのため，反応が 1 種類増えるごとに，シミュレーショ
ンの使用する加算器が 1 つずつ増え，この加算の部分の長さが増加することになる．
phase 4.
Step 3. : (0, 1) の単精度浮動小数型一様乱数 r2 を計算する．
• 一様乱数生成モジュールが出力する 23 ビットの乱数の上位に (00111111)(2) を付加し，(1, 2)の浮動小数点実数の一様乱数を得る• 前項の一様乱数から 1.0(浮動小数点型 3F800000(16) ) を減算し，(0, 1) の浮動小数点実数の一様乱数を得るphase 5.
Step 3. : 15 ビットの乱数から対数テーブルを引いて 1/ ln(1/r1 ) の値を得 (r1 は (0, 1) の単精度浮動小数点型実数)，次の反応が起こるまでの時間の逆数 1/τ を計算する．
1a0=
τ ln(1/r1 )• 一様乱数生成モジュールが出力する 15 ビットの乱数から，その値をアドレスと見なし，対数格納テーブルから値を読み出す• a0 の値と対数テーブルから得られた値を乗算し，1/τ を得る
42

4. 設計と実装

4.3. 実装

phase 6.
Step 4. : phase 2. で求めた aν と phase 3 で求めた a0 を比較し，次に発生する反応 µ を求める．
µ−1Xaν < r2 a0 6

µ=1

µ
X

aν

ν=1

• a0 の値と phase 4. で生成した (0, 1) の実数を乗算する• 浮動小数点比較器でその値が 0, a1 , a2 , a3 のどの間にあるのかを求め，(0, a1 ] にあるならばµ = 0，(a1 , a2 ] にあるならば µ = 1 というように，µ を求める
phase 7.
Step 5. : 発生した反応 µ に応じて，phase 1. で使用した分子数 Xν を増減する．ここで，phase 5.
の 1/τ と分子数 X2 , X3 を外部メモリに出力する．
X2 + 1(if µ = 0)

X2 − 1andX3 + 1(if µ = 1)

X3 − 1andX4 + 1(if µ = 2)

X2 − 1andX4 + 1(if µ = 3)

• 変数保存テーブルから分子数を読み出し，phase 6. で得られた µ をもとに，分子数を加減算する• 結果の分子数を次サイクル用に phase 1. へ出力する• 結果の 1/τ と分子数を外部に出力する実装したモジュールごとの処理の流れを示すと，図 4.4 のようになる．

使用している整数・浮動小数演算器はパイプライン化されているため，1 クロックごとに連続的に値を投入することが可能である．しかし，浮動小数点加減算は，出力の結果を次の入力で使うため，単体の演算器に入力値を選択するレジスタが必要になり，結果が出るまで 6 クロック (入力時にセレクタが追加され 1 クロック増える) 要している．
各演算モジュールにより，シミュレーション 1 サイクル (反応が発生する間隔) に必要なクロッ
ク数は，図 4.4 より，40 クロックとなる．
各演算モジュールを 1 つずつ使用してこの Lotka 反応モジュールを構成すると，使用率が最も高いのは浮動小数点乗算器で，1 つのシミュレーション計算に 1 サイクル 40 クロックのうち，6クロック分使用している．この構成で単純に並列化すると，40/6 = 6.7 より最大で 6 つの並列化が可能となる．しかし，この方法はモジュールに投入する値を制御する回路が複雑になってしまう．また，浮動小数点乗算器以外のモジュールの使用率が低く，シミュレータをさらに複数用意して共有することは，回路がさらに複雑になり，信号線を引き回すために起こるオーバヘッドの増加や動作周波数の下落，設計の複雑化などのデメリットが考えられる．

43

44

Lotka

memory I/O

integer
multiplierinteger tofloating
floating
multiplierfloating
add & subfloating
comparator
random 23
generator
random 15
generator
lntable

0

1

2

3

5

phase 1

4

6

7

8

9

11

phase 2

10

12

13

14

15

16

17

18

19

20

22

23

phase 3

21

24

25

26

27

29

30

phase 4

28

31

32

34

36

37

phase 635phase 5

33

38

39

40

4. 設計と実装
4.3. 実装

図 4.4: シミュレーション 1 つが 1 サイクル中でモジュールを使用する部分

4. 設計と実装

4.3. 実装

図 4.3 から，phase 1. と phase 2. における整数乗算および浮動小数点乗算，整数-浮動小数変換器はそれぞれ並列動作させることができることがわかる．また，実装したすべてのモジュールは
パイプライン動作し，かつ入力動作が 1 クロックで完了することから，演算数分の演算器を使用
し，1 サイクルにつき 1 クロックだけモジュールを使用するようにそれぞれのモジュールを接続
し，Lotka 反応モデルモジュールをパイプライン化する．表 4.2 に，パイプライン Lotka モデルモ
ジュールに必要なモジュールのうち，複数使用するモジュールとその数を示す．
各モジュールは面積が比較的小さく，Vertex-II に十分収めることができ，また設計と信号線の制御の容易性の面からもこの設計は妥当であると言える．
表 4.2: パイプライン化 Lotka 反応モジュールに必要な演算器

32 bit 整数乗算器
32 bit 整数-浮動小数変換器
単精度 浮動小数点加減算器
単精度 浮動小数点乗算器

2
4
4
6

各乱数生成器，浮動小数点比較器，反応実行モジュールはパイプライン化が可能である．このパイプライン化により，あるモジュールの出力が次のモジュールの入力に直接信号を接続することができるため，浮動小数点加減算器の入力時にセレクタが不必要になり，図 4.4 よりも 1 サイ
クルに掛かるクロック数が少なくなっている．37 クロック 1 サイクルのパイプライン Lotka 反応
モデルモジュールの 1 つのシミュレーション計算に使用する演算モジュールとその処理の流れを図 4.5 に示す．
パイプライン化により，それぞれの演算モジュールは 1 サイクルにつき 1 クロック使用されるため，データを 1 クロックずつずらして入力することで，1 つの Lotka 反応モデルモジュールは37 個のシミュレーションを同時に計算することが可能である．

4.3.6

変数保存テーブル

Lotka 反応モデルモジュールでは，それぞれのシミュレーションサイクルの最初で現時点の分子数が入力され，サイクルの最後で次の反応までの時間 τ と，反応の結果の分子数が出力される．
そのため，実装したパイプライン化 Lotka 反応モジュールでは，反応前の分子数を保存しておく必要がある．実装した Lotka 反応モデルは 4 つの分子からなるモデルであり，37 段パイプラインのため 32 ビット幅で深さ 37 の変数格納用テーブルが 4 つ必要になる．このテーブルは出力された分子数を対応する場所に書き込み，その時点の入力に対応する分子数を Lotka 反応モジュールに読み出す，という動作をするため，Vertex-II に搭載されているデュアルポートの BlockRAM を使用して構成した．変数保存テーブルを構成するために必要な BlockRAM のサイズは，
32[bit] × 37[pipeline] × 4[molecules] = 4736[bit] = 592[Byte](4.7)となる．この変数保存テーブルの 2 つのポートの一方をポート A，他方をポート B とすると，ポート A は読み出し専用，ポート B は書き込み専用に使用し，それぞれアドレスを 0 から 36 までの値をインクリメントしながらそれぞれのポート A で読み出し，ポート B で書き込みを行う．浮動小数点比較器が µ の値を出力してから，メモリから分子数をポート A から読み出し，反応実行モ

45

46

Lotka

memory I/O

lntableinteger
multiplier 1
integer
multiplier 2
integer tofloating 1
integer tofloating 2
integer tofloating 3
integer tofloating 4
floating
multiplier 1
floating
multiplier 2
floating
multiplier 3
floating
multiplier 4
floating
multiplier 5
floating
multiplier 6
floating
adder 1
floating
adder 2
floating
adder 3floating
substractor
floating
comparator
random 23 bit
generator
random 15 bit
generator

0

1

3

4

phase 1

2

5

6

7

9

10

phase 2

8

11

12

13

14

15

16

17

18

20

phase 3

19

21

22

23

25

phase 4

24

26

27

28

30

31

phase 529

33

phase 6

32

34

35

36

37

4. 設計と実装
4.3. 実装

図 4.5: 37 段パイプラインの Lotka 反応モデルモジュールが使用する部分

4. 設計と実装

4.3. 実装

ジュールを通して反応による分子数の増減を行った上で，読み出したアドレスをポート B で書き換えるため，2 つのポートのアドレスが一致することはなく，RAW ハザード，WAR ハザードは発生しない．

4.3.7

入出力モジュール

Lotka 反応を ReCSiP 上で動作させる上で，現実的な実装方法を述べる．
対数テーブルはデュアルポートで構成されており，1 クロックで 2 つの値を読み出すことができる．
Lotka 反応モデルモジュールの面積は Vertex-II 全体に比べて小さいことが考えられる．BlockRAMの総量に比べて対数テーブルの占める割合は大きいため，Lotka 反応モジュールを 2 つ 1 組で構成
し，BlockRAM を効率的に使用することにする．
initialize
x1

x3
x2

32

c1
x4

randomseed 1c3
c2

randomseed 2

c4

32

32

32

15 bitrandomgenerator
15

15 bitrandomgenerator
15

lntable (dual port)
32

32

Lotka
alpha
Lotka
bravo
Lotka
Simulator
32 32 32 32x1 x2 x3 x4

32

1/tau

32 32 32 32x1 x2 x3 x4

32

1/tau

図 4.6: Lotka 反応モデルモジュールを 2 つ使用したシミュレーションモジュール以降は 2 つの Lotka 反応モジュールを 1 組として，Lotka 反応セットモジュールとして扱う．Lotka
反応セットモジュールの入力は，リセット信号が無効になった後，初期化値として 4 つの分子の数の値 (x1 〜 x4 ) の値と 2 つの乱数種を 32bit の値で与える必要がある．入力値は 1 組を 1 クロック47

4. 設計と実装

4.3. 実装

ずつ 37 クロックに渡って入力する．38 クロック目から 1 クロックごとに出力として反応後の分子の値の組 (x1 〜 x4 ) が Lotka 反応モデルモジュール alpha と bravo より合計 2 組ずつ出力される．
Lotka 反応モデルモジュールの 1 つの反応は 37 クロックごとに結果を出力する．37 段パイプライ
ンで全ての演算器を 100% 使用するため，1 クロックごとに分子数と 1/τ が出力される．毎回の出力の組とともに，前回の反応からその反応が発生するまでの時間の逆数 1/τ が単精度浮動小数点として出力される．
現行の設計では浮動小数点の除算器が実装されておらず，τ の値はホスト側で計算しなければならないため，複数サイクルの反応をまとめてメモリに出力する，という方法が使用できない．そのため，1 サイクルごとの結果をメモリに記録する必要がある．また，結果の出力において，結果を書き出す SRAM が 4 つであることから，分子数 x2 , x3 を 16 ビットで併せて 32 ビット，32 ビットの 1/τ を 1 つの SRAM に書き出し，シミュレーション結果を検討する．

48

第5章

評価
5.1

ソフトウェアでの計算時間

Lotka 反応セットモジュールを Vertex-II 上で実行する上で，現在あるソフトウェアシミュレー
ションと比較した際の性能向上について述べる．ソフトウェアで Gillespie の確率モデルシミュレー
ションアルゴリズムの Lotka 反応モデルについてシミュレーションを行った場合の処理能力について評価を行った．
ソフトウェアによって Lotka 反応モデルの確率的シミュレーションを行うプログラムを C 言語で記述し，50 万反応 (50 万サイクル) のシミュレーションを 1 回実行させた場合に要する時間を表 5.1に示す．また，1 秒当たりに実行できるサイクル数 (反応数) を表に含めた．処理時間はプログラム
を 10000 回実行し，平均を取ることで算出した．1 秒当たりに実行できるサイクル数 S [Mcycle/sec]
は，以下の式から算出した．
0.5[Mcycle]
(5.1)
S =
処理時間 [µsec]

表 5.1: ソフトウェアでの処理能力
プロセッサ

メモリ

環境

処理時間 [µsec]

cycle/sec

Xeon 2.8GHz Dual(HT off)
Pentium4 2.53GHz
PentiumIII 800MHz

4GB
1GB
512MB

Linux2.4.21+gcc2.95.3
Linux2.4.21+gcc2.95.3
Linux2.4.21+gcc2.95.3

214219244607512386

2.33M
2.04M
0.98M

AthlonXP 2800+ 2.08GHz
Pentium4 2.6GHz
Celeron 450MHz

2GB
1GB
512MB

FreeBSD 4.8+gcc2.95.3
FreeBSD 4.8+gcc2.95.3
FreeBSD 4.8+gcc2.95.3

172226
192496740459

2.90M
2.59M
0.68M

UltraSPARCIIIcu 1.2GHz
UltraSPARCII 300MHz

4GB
2GB

Solaris8+gcc2.95.3
Solaris8+gcc2.95.3

5559071847605

0.90M
0.27M

ソフトウェアで反応を 50 万サイクルシミュレーションした結果は，以下のようになる．比較対象としたソフトウェアプログラムでは，C 標準ライブラリを用いて乱数を生成している．図 5.1に，x 軸方向に時間，y 軸方向に分子数を取り，X2 (被食者) と X3 (捕食者) の数の変化をプロットしている．x 軸の単位は cν の単位によって変わるため，明示していない．一般的な化学反応の場合，49

5. 評価5.2. Vertex-II への実装結果

単位は [sec] である．図 5.1 のように，Lotka 反応モデルでは 2 つの分子数が時間の進行とともに振動する．
Software result1
3000

X2
X3

2500

population level

2000

1500

1000

500

0

0

2

4

6

8

10

time

図 5.1: ソフトウェアで出力される分子数変化
このソフトウェアプログラムについて，乱数の種を変えたときにどうなるかを述べる．ソフト
ウェアプログラムを乱数の種を変えて実行した時のグラフを図 5.2 に示す．
乱数種を変えるとシミュレーション結果にも変化が現れることがわかる．Gillespie のアルゴリ
ズムは確率モデルを用いて化学反応モデルのシミュレーションをするものであるため，サイクルごとに異なる乱数が与えられるならば，異なる結果が現れる．
シミュレーション結果は起こる反応に変化が現れるだけでなく，シミュレーションの時間も変化する．図 5.1 と図 5.2 は同じ反応を 50 万サイクルの間シミュレーションを実行した結果である
が，50 万回反応が発生するまでの時間は異なっている．
図 5.3 と図 5.4 で 50 万サイクルにかかる時間に違いが発生している．最初のシミュレーション
をシミュレーション A，次のシミュレーションをシミュレーション B とすると，50 万サイクル経
過の時間について，A は 9.936，B は 10.143 である．このように，次に反応が発生するまでの時
間 τ も乱数を用いて計算するため，シミュレーションを特定サイクル実行しても，シミュレーショ
ン時間が異なることがわかる．

5.2 Vertex-II への実装結果
現在，Gillespie のアルゴリズムを用いて基本的な Lotka 反応モデルのシミュレートを行うために，4.3 章で述べた設計・実装法を用いて Lotka 反応セットモジュールの論理合成と Vertex-II への50

5. 評価5.2. Vertex-II への実装結果

Software result2
3000

X2
X3

2500

population level

2000

1500

1000

500

0

0

2

4

6

8

10

time

図 5.2: ソフトウェアの乱数種を変えたときの変化

Software result1

Software result2

2400
2200
2000

2000

1800

1800

1600

1600

1400
1200
1000

1400
1200
1000

800

800

600

600

400

400

200

9.6

9.8

10
time

10.2

X2
X3

2200

population level

population level

2400

X2
X3

200

10.4

図 5.3: Software result1 の末端

9.6

9.8

10
time

10.2

図 5.4: Software result2 の末端

51

10.4

5. 評価5.3. シミュレーション結果

配置配線の結果を述べる．
HDL の記述は Verilog-HDL を使用し，論理合成/配置配線は Xilinx 社の ISE6.1.03i を用いて行った．
Lotka 反応セットモジュールに入出力モジュールを加えた状態で論理合成/配置配線を行った場
合のリソース使用量および動作周波数を表 (5.2) に示す．
表 5.2: 入出力付き Lotka 反応セットの合成結果
スライス数8376 (24.79%)

組み込み乗算器

ブロック RAM

動作周波数60 (41.67%)

66 (45.83%)

78 [MHz]この Lotka 反応セットモジュールは 37 クロックを 1 サイクルとして，1 クロックごとに分子数
X2 ，X3 および反応発生までの時間の逆数 1/τ を 2 組ずつ出力する．このことから，Lotka 反応モ
デルについてのシミュレーションを 37 × 2 = 74 個並列に実行することが可能である．
Lotka 反応セットに含まれている各モジュールの合成結果を以下の表 (5.3) に示す．
表 5.3: 各モジュールの合成結果

モジュール

スライス数

組込乗算器

BlockRAM

動作周波数

パイプライン

単精度浮動小数点加減算器
整数乗算器
Lotka 反応処理モジュール
一様乱数生成器 (15 ビット)一様乱数生成器 (23 ビット)
対数テーブル
単一の Lotka 反応モジュール

384 (1.13%)
159 (0.471%)
106 (0.313%)
97 (0.287%)
16 (0.0473%)
16 (0.0473%)
348 (1.03%)
3974 (11.76%)

0 (0%)
4 (2.78%)
3 (2.08%)
0 (0%)
0 (0%)
0 (0%)
0 (0%)
30 (20.83%)

0 (0%)
0 (0%)
0 (0%)
0 (0%)
0 (0%)
0 (0%)
60 (41.67%)
4 (2.78%)

83 [MHz]
110 [MHz]78 [MHz]
97 [MHz]
94 [MHz]
87 [MHz]
79 [MHz]80 [MHz]

5段
7段
4段
1段
1段
1段
1段
37 段

入出力付き Lotka 反応セット

8376 (24.79%)

60 (41.67%)

66 (45.83%)

78[MHz]37 段

単精度浮動小数点乗算器

5.3

シミュレーション結果

Lotka 反応セットモジュールを Verilog-HDL で記述し，シミュレーションを行った結果を以下に述べる．シミュレーションは synopsys 社 vcs を用いて行い，シミュレーションログファイルを出力した．シミュレーションログを perl5.6.1 で 74 のファイルに分割し，Gnuplot を用いてグラフ化した．シミュレーションログで出力されるシミュレーション結果は X2 , X3 の分子数，およびその反応までの時間の逆数 1/τ であるため，ログファイルをシミュレーションごとのファイルに分割する際に逆数を取ることで τ を求め．サイクルごとに τ を加算し，その分子数はシミュレーショ
ン時間でいつの状態なのかを計算している．初期の分子数 Xν と，システム温度と分子の物理的性質のみで決定する確率反応速度定数 cν は以下の値を使用している．

52

5. 評価5.3. シミュレーション結果

c1 = 0.0002
X1 = 100000

c2 = 0.01
X2 = 1000

c3 = 10
X3 = 1000

c4 = 10
X4 = 0Hardware result 0
2000

X2
X3

1800

population level

1600

1400

1200

1000

800

600

400

0

2

4

6

8

time

図 5.5: 実装した Lotka 反応セットモジュールのシミュレーション結果

53

10

5. 評価5.3. シミュレーション結果

図 5.5 は Lotka 反応モジュール A のうちの 1 つについて，50 万サイクルを実行したときのグラフである．実装した Lotka 反応セットモジュールは，これと同時に 74 個のシミュレーションを行っている．それらの結果の一部を図 5.6 に示す．
Hardware result 1

Hardware result 2

2200

2000

X2
X3

3500

1800

3000

1600

2500

1400

X2
X3

1800

population level

1200
1000
800

population level

1600

1400

1400
1200
1000
800

600

population level

2000

1600
population level

Hardware result 4

4000

X2
X3

2200

1800

2000

1500

1000

1200

1000

800

600

400
200

Hardware result 3

2400

X2
X3

2000

500

400
0

2

4

6

8

200

10

0

2

4

6

8

0

10

600

0

2

4

6

time

time

time

Hardware result 5

Hardware result 6

Hardware result 7

3000

3500

X2
X3

2500

8

3000

1600

2500

1400

0

2

4

6

8

10

time
Hardware result 8

1800

X2
X3

400

10

4000

X2
X3

X2
X3

3500

3000

2000

1500

population level

1500

population level

population level

population level

2000
1200

1000

2500

2000

1500

1000

500

0

0

2

4

6

8

1000

800

500

600

0

10

0

2

4

6

8

400

10

1000

500

0

2

4

6

time

time

time

Hardware result 9

Hardware result 10

Hardware result 11

2000

1600

X2
X3

0

10

1100
1000

1600

1500

1000

800

2

4

6

8

600

10

2

4

6

8

0

10

0

2

4

time

Hardware result 13

6

8

200

10

0

2

4

time

Hardware result 148

10Hardware result 16

3500

X2
X3

6
time

Hardware result 15

3500

X2
X3

1000

400
0

time

2500

1200

600

500

700

0

1400

800

800

400

10

X2
X3

2000

population level

1200

900

600

8

2200

2000
population level

population level

population level

1000

6

1800

1300

1200

4
time

1600

1400

2

2400

X2
X3

2500

1400

0

Hardware result 12

3000

X2
X3

1500

1800

8

3500

X2
X3

3000

3000

3000

2500

2500

2500

X2
X3

2000

1500

population level1000

population level

1500

population level

population level

2000

2000

1500

1000

1000

500

500

2000

1500

1000

500

0

0

2

4

6
time

8

10

0

0

2

4

6

8

0

10

500

0

time

2

4

6

8

10

0

0

time

2

4

6

8

10

time

図 5.6: 並列動作する Lotka 反応シミュレーションの一部

5.3.1

ハードウェアとソフトウェアの動作比較

Lotka 反応セットモジュールが正常に動作し，それとソフトウェアによるシミュレーション結果が一致することを確認する．
図 5.1 や図 5.5 より，ソフトウェア，ハードウェアともに，X2 ，X3 の数が振動していることがわかる．しかし，ソフトウェアでは標準ライブラリ，ハードウェアでは M 系列と，乱数生成規則が異なり，また，ハードウェアの浮動小数点演算器には丸め機能が無い．それらが出力される結果にどのような影響を及ぼすかを述べる．
ハードウェアとソフトウェアで乱数生成規則が異なるため，初期値や乱数種に同じ値を設定してシミュレーションを実行しても，少数の試行の結果だけでは，同じ形状の波が形成されることはない．乱数生成規則の有効性の有無については，今後検証を行っていくこととし，ここではハードウェアで実装したシミュレータが正しい動作をしていることを述べる．
このことを確認するため，ソフトウェア側にハードウェアと同等の 32 ビット M 系列乱数生成関数を実装した．シミュレーションにおける τ は，a0 と乗算した後に出力され，以降のシミュレー
ションで再び使用されることはない．そのため，τ による誤差の蓄積は無い．
54

5. 評価5.4. 性能評価

しかし，(0, 1) の実数から求められる µ は，発生する Rµ の決定に関わり，反応によって分子数が変更され，それが次の a0 の決定に使用されるため，ここで誤差が発生した場合，反応の進展に合わせて誤差が伝搬してしまう．丸め誤差の有無による誤差伝搬はここに現れる．
10 万サイクルの Lotka 反応モデルについて，乱数生成規則を M 系列とし，乱数種に同じ値を与え，使用する乱数列を等しくした上で，ソフトウェアとハードウェアの両実装を比較した．
Lotka 反応モデルの 10 万サイクルの反応について，ソフトウェアによるシミュレーション結果
を図 5.7 に，ハードウェアによるシミュレーション結果を図 5.8 に示す．
ハードウェアのシミュレータには浮動小数点除算器が実装されていないため，Perl を使用し，シミュレーションログから 1/τ の値を取り出し，逆数を取ることで τ の値を計算している．
SoftwareHardware

1600
1500
1400

1400

1300

1300

1200

1200

1100
1000

1100
1000

900

900

800

800

700

700

600

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

X2
X3

1500

population level

population level

1600

X2
X3

600

1.8

time

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

time

図 5.7: ソフトウェアによる実行結果

図 5.8: ハードウェアによる実行結果

ソフトウェア，ハードウェアのシミュレーションログから，発生する反応を示す µ の値はシミュ
レーションを通して同一であり，与えられる乱数列が同一の場合，丸め誤差による発生反応の変化，およびその誤差伝搬は生じていない．反応発生時間 τ は，対数テーブル内のデータは C で作成されたことから，対数テーブルからデータを読み出す乱数生成器側が同じであるので，ハード
ウェア，ソフトウェアともに同じ値が出力されている．τ は a0 との乗算およびログ出力後の除算により，出力される値に小さな違いがある．しかし，τ は誤差伝搬しないため，シミュレーション
を通して結果が大きく異なることはない．
このことから，ソフトウェアで実行するシミュレーションとハードウェアで実行するシミュレー
ションに大きな誤差は発生していないことが確認できた．
シミュレーションの結果が正しいことを判断するためには，これらのシミュレーションを多数繰り返し，平均や分散を取って調査する必要がある．

5.4

性能評価

実装した Lotka 反応セットは，37 クロックで 1 つの反応を処理することが可能である．Lotka 反応セットモジュールを 78[MHz] で動作させたときのシミュレーション 1 つあたりの処理時間は，
78[MHz]/37 = 2.11M[cycle/sec] となり，Pentiun4 2.53GHz とほぼ同等の処理速度を実現できる．
また，Lotka 反応セットモジュールはパイプライン化されており，37 個のシミュレーションを同時に行うことが可能である．さらに，内部には Lotka 反応モデルモジュールが 2 つ構成されて55

5. 評価5.4. 性能評価

おり，同時に動作させることができる．つまり，Lotka 反応セットモジュールは 37 × 2 = 74 個のシミュレーションを同時に処理することができる．そのため，Lotka 反応セットモジュールの全てを使って処理を行った場合，74 × 2.11 = 156.14M[cycle/sec] の処理速度を実現できる．これは，ソフトウェアで実行したときにもっとも高速だった AthlonXP 2800+の 53.84 倍である．

56

第6章

結論
6.1

確率モデルを FPGA でシミュレーションする有効性

本研究では，FPGA を用いた生物学的計算向けのプラットフォームである ReCSiP をターゲット
として，Gillespie のアルゴリズムを用いた確率モデルによる代謝シミュレーションコアの開発を行った．Gillespie が例示している確率モデルから Lotka 反応モデルをシミュレーションする回路を実装し，その評価を行った．その結果，Lotka 反応モデルについて AthlonXP2800+の 53.84 倍の速度でシミュレーションが可能であることがわかった．生物学計算をはじめとしたモンテカルロ
シミュレーションは大きな並列性を持つため，ReCSiP ボードのような FPGA による並列化，高速化の有効性を示すことができた．

6.2

今後の課題

6.2.1

除算器の実装による並列化

面積，BlockRAM，組み込み乗算器に余裕があるため，Lotka 反応セットモジュールを ReCSiP
上を 2 つ構成することも可能である．現行の実装では Lotka モジュール 1 つにつき 2 つの SRAM
へ結果を出力しているため，Lotka 反応セットモジュールが 1 つに限られてしまう．この解決法として，浮動小数点除算器を実装することで出力を間引くことが考えられる．結果格納メモリへの入出力を取り除くことによって，組み込み乗算器の数の上限まで，Lotka 反応セットを 2 つ ReCSiP上に実装することが可能である．このため，最大で 37 × 2 × 2 = 148 個のシミュレータを同時に実行することが可能になる．
また，浮動小数点除算器により，Gillespie の例示している Lotka 反応以外の “複雑な反応モデ
ル” も実装でき，より幅広いモデルのシミュレータを構築することが可能になる．

6.2.2

精度向上

対数テーブル
現行のシミュレータでは，対数テーブルが 32 ビット幅 15 ビットエントリであり，格納データ数が 215 = 32768 個しかなく，次の反応までの時間 τ の計算が雑なものになっていることは明らかである．最大で 215 の周期の実数の対数は乱数として用いるには不十分である．この問題に対して，ホスト PC との協調処理を行って解決を図る方法が考えられる．

57

6. 結論6.2. 今後の課題

すなわち，τ = ln(1/r1 )/a0 のうち，1/ ln(1/r1 ) を計算する役割をホスト PC が担当する処理を行う．このことによって，ハードウェアによる乱数生成に比べ，ホスト PC 側のソフトウェアで乱数生成アルゴリズムを使用し，信頼性の高い乱数を得ることが可能になる．実装方法は，以下のよ
うに対数テーブルを複数に分割し，読み出しを行うテーブル，値の更新を行うテーブルを構成することが考えられる (図 6.1)．読み出し，更新のテーブルを一定クロックごとに切り替えるセレク
タを付加することで，擬似的に (0, 1) の 32 ビット浮動小数点乱数を得ることが可能になる．この場合，ホスト PC から書き込まれる値もランダムになるため，15 ビットの乱数生成器は不要になり，また浮動小数点減算器も不要になる．
1/tau

lntableFrom Host PC

with control unit
32

32 bit width
15 bit entry

lntable 4lntable 2
32 bit width
13 bit entry

32 bit width
13 bit entry

lntable 1

lntable 3

32 bit width
13 bit entry

32 bit width
13 bit entry

32

32

32

32

alpha

bravo

charliedelta

Lotka_set_A

Lotka_set_B

図 6.1: ホストと協調動作する対数テーブル

乱数生成器の改良
乱数生成器は現在 32 ビットの LFSR による M 系列疑似乱数生成器である．この乱数生成期は最も長いものでも周期が 232 − 1 であり，複雑なモデルの数百万反応を超えるシミュレーションを視野に入れると，現行の乱数生成器では対応できないことが容易に予想できる．開発中の ReCSiP2
では自然乱数生成器が搭載される予定だが，自然乱数生成器は動作周波数が低く，全ての乱数の生成には利用できない．そのため，自然乱数生成器を使って乱数の種を生成し，そこから疑似乱数を作っていくことが考えられる．ハードウェアでの乱数生成は M 系列を用いたものが主である
が，この場合も，LFSR を 90 ビット程度の長いものを使用し，その LFSR を複数並べ，それらか58

6. 結論6.2. 今後の課題

ら 8 ビットずつ読み出し，それらを組み合わせて 32 ビットの長さの乱数を生成する，といった方
法で，乱数の周期性を伸ばす手法を取り，十分な検証を行う必要がある．図 6.2 にその一例を示す．これは，90 ビットの LFSR を 4 つ使用し，毎クロック最上位 8 ビットを取り出し，それらを組み合わせることで 32 ビットのビット列を生成するモジュールである．
90bit LFSR

90bit LFSR

8

8

90bit LFSR

90bit LFSR

8

8

8bit / clock
shift register

8

8

8

8

32bit random number

図 6.2: より周期の長い乱数生成器の一例

6.2.3

モデル自動生成

Lotka 反応モデルを実装したことで，化学反応モデルのモンテカルロシミュレーションは ReCSiPのような FPGA を搭載したデバイスで大きな効果が期待できることがわかった．Gillespie は Lotka
反応モデル以外にも複雑な化学反応モデルを例示しているし，STOCKS のようなシミュレータは入力ファイルを記述することで細胞のシミュレーションに対応している．
Gillespie のアルゴリズムのシミュレーションコアは，モデルごとの反応処理モジュールに加え，面積的に多くを占有する算術演算器が以下の数だけ必要になることがわかっている．反応分子の組み合わせ数を x とすると，hν を決定するための整数乗算器を x，hν と確率速度定数 cν を乗じる
ための浮動小数点乗算器を ν，それらを累積加算するための浮動小数点加算器が x − 1 必要になる．
ReCSiP の Vertex-II XC2V6000 に，反応処理モジュールおよびこの数の算術演算器が搭載できる範囲内であれば，Lotka 反応セットモジュールで示したように各サイクルをパイプライン化し，動作中は全てのモジュールの使用率が 100% でモジュールを有効に使用することが可能であり，さらに大きな並列性を得ることが期待できる．搭載できないほど大きな処理の場合は，演算器や使用モジュール数を制限した上で，それらのスケジューリングを行い，並列化を図る予定である．
モデルに合わせたシミュレーションモジュールを自動生成することから始め，将来的には，SBML(Systems
Biology Markup Language) で記述されたモデルから確率モデルシミュレータの Verilog コードを生成できるように展開する予定である．

59

謝辞
本研究の機会を与えてくださり，また常にご指導，ご助言をいただきました
天野英晴教授に，深く感謝いたします．毎回のアドバイスが構想の発展に繋がり，研究を続けるための大きなモチベーションになりました．
配属後の研究初期段階から多くの面でお世話になり，また論文チェックも行ってくださった
長名 保範氏に感謝致します．実装の楽しさを学びました．
実装の良否について多くの助言をいただき，また論文の査読を行ってくださった福島 知紀氏福田 静人氏に、感謝致します．たくさんの画期的なアイデアを発想できました．
研究を進めるにあたり，生物学分野，ソフトウェア工学分野からの多くの有益な知識を頂きました，現情報工学科斉藤博昭研究室，生命情報学科榊原研究室の松井 洋氏に感謝致します．
日々共に作業し，励ましてくださった
住吉 正人氏黒瀧 俊輔氏に感謝します．深夜に寂しくて泣き出さなかったのはお二人がいてくれたおかげです．
日々の研究に際し，多くの助言と，生活のアドバイスをいただいた
天野研究室の全てのみなさまに感謝します．この１年間で知識，サバイバビリティともに大きく成長できました．
家族から，祖母の川上 榮子氏に感謝します．彼女の有形無形の支援が無ければ，私はここまで来られませんでした．ありがとう．
今まで私と関わり，ここまで引っ張ってくれた全ての皆様に限りない感謝を．
ありがとうございました．

60

参考文献
[1] Bruce A. Barshop, et al. Analysis of numerical methods for computer simulation of kinetic processes: Development of kinsim — a flexible, portable system. Analytical Biochemistry, Vol. 130,
pp. 134–145, 1983.
[2] Mendes P. Gepasi: a software package for modelling the dynamics, steady states and controlof biochemical and other systems. Computer Applications in the Biosciences, Vol. 9, No. 5, pp.
563–571, Oct. 1993.
[3] I Goryanin, et al. Mathematical simulation and analysis of cellular metabolism and regulation.
Bioinformatics, Vol. 15, No. 9, pp. 749–758, Sep. 1999.
[4] Masaru Tomita, et al. E-cell: software environment for whole-cell simulation. Bioinformatics,
Vol. 15, No. 1, pp. 72–84, Jan. 1999.
[5] J. Schaff, et al. A general computational framework for modeling cellular struc ture and function.
Biophysical Journal, Vol. 73, pp. 1135–1146, Sep. 1997.
[6] Daniel T. Gillespie, et al. Exact stochastic simulation of coupled chemical reactions. The Journalof Physical Chemistry, Vol. 81, No. 25, pp. 2340–2461, Dec. 1977.
[7] Andrzej M. Kierzek. Stocks: Stochastic kinetic simulations of biochemical systems with gillespie
algorithm. Bioinformatics, Vol. 18, No. 3, pp. 470–481, Mar. 2002.
[8] Nicolas Le Novère, et al. Stochsim: modelling of stochastic biomolecular processes. Bioinformatics, Vol. 17, No. 6, pp. 575–576, Jun. 2001.
[9] 長名保範ほか. Reconfigurable hardware による細胞シミュレータの高速化手法. ARC2002
148-8, 情報処理学会, May. 2002.
[10] Carl Firth, Nicolas Le Novère, and Tom Shimizu.
ftp://ftp.cds.caltech.edu/pub/dbray/, 2001.

STOCHSIM,the stochastic simulator.

[11] Stephen F, et al. Basic local alignment search tool. Journal of Molecular Biology, Vol. 215, pp.
403–410, Oct. 1990.
[12] Jeffrey M. Arnold, Duncan A. Buell, and Elaine G. Davis. Splash 2. In Proceedings of the fourthannual ACM symposium on Parallel algorithms and architectures, pp. 316–322. ACM Press, 1992.
[13] Ebisuzaki T, et al. Grape: special purpose computer for simulations of many-body systems. In
Proceedings of the twenty-sixth Hawaii International Conference on System Science, Vol. 1, pp.
134–143. IEEE, Jan. 1993.
61

参考文献

[14] 田村友紀ほか. 細胞系譜構築システムの reconfigurable system による高速化. VLD2001 122,
信学報, Jan. 2002.
[15] PCI Special Interest Group. PCI Local Bus Specification Revision 2.3, Oct. 2001.

62