E-CELL project
E-CELLSystemによるシアノバクテリア
概日リズムの表現法と経路予測システムの研究

政策・メディア研究科 三由文彦


INDEX

abstract
1 目的
2 シアノバクテリアの概日リズム
3 方法
3.1 スモルコフスキー式
3.2 シグナル伝達と遺伝子発現系のモデル表現法
3.3 E-CELL内部の計算方法
3.4 サンプルデータを用いた決定論と確率論モデルの比較
4 パラメータ推定
5 結果
6 考察
7 展望
7.1 計算速度の向上
7.2 モデルルール生成環境の構築
謝辞
参考文献





概日リズムは比較的独立したシステムで明示的な入力情報は光だけであり、少 数のコンポーネントから成っていることからコンピュータによる経路予測がし やすいと思われる。この経路予測を可能とするために、反応の活性化エネルギー を変数として含むスモルコフスキー式をベースにした表現法を開発した。この 表現法では、決定論と確率論シミュレーションを同一のパラメータで実行可能 であり、モデルの挙動にあわせて自動的にシミュレーション手法変換すること も可能になる。今回、概日リズムの遺伝子制御ネットワークやタンパク質相互 作用の予測を行なう目的で、この表現法を用いたシアノバクテリアのモデルを E-CELL System上で構築した。これらを用いて発現の時系列データやプロモー ターなど塩基配列情報から推測される遺伝子制御ネットワークのグループをも とに予測を行なう。


1 目的


ポストゲノム時代の幕開けに伴って、遺伝子制御ネットワークを構成するmRNAや 蛋白質などの時系列データが網羅的・大量に観測されることが予測される。遺伝 子ネットワークを予測する手法として、ブーリアンネットワークによるもの [17,18]や、S-systemによるもの[19,20]などが挙げられ一定の成果を得ているが、 詳細な遺伝子間相互作用を予測するレベルにまでは至っていない。

そこで今回、これらの時系列データを用いて私は遺伝子制御ネットワークにおい て経路予測する手法の確立を目指すこととした。その際、比較的独立したシステ ムを持ち、外部からの明確な入力情報が光だけであり、さらに少数の構成要素か らなる概日リズムが、その経路予測システム構築基盤として適当であると考えた。

概日リズムの数理モデルに関する研究は、すでにキイロショウジョウバエなどに おいて盛んに行われている[12,13,14]。これらは従来あった代謝の数理モデルと 同様に遺伝子の発現を決定的に捉えて反応速度式によって、数理モデル化を行っ ている。しかし、これらの決定論的な方法に関しては、ターゲットとする生物の 未解明部分や諸説を仮説に基づいたモデルが立てられている段階に留まっており、 実際にこれらを利用した網羅的なパラメータ推定によるネットワーク予測などの 試みはほとんど行われていない。

一方、決定論モデルでは、反応が同時に、多数、平均的に起こっている状態を表 すことに向いているため、この反対つまり、ランダム、少数、局地的に反応が起 こっている遺伝子制御ネットワークやシグナル伝達系に関して、その挙動を詳細 に再現することは難しい。また細胞内では、数種類のタンパク質が複合体を形成 することによって反応を起こす例がある。この時、数十のタンパク質が起こす反 応を再現しようとすると、その組合せは大変な数となり、その結果、時間幅毎に 反応速度式の計算をしなければならない決定論モデルでは、計算量自体も大変な 数になってしまう。このため決定論モデルでは、組合せをある程度省略した形で 表現せざるを得なくなってしまう。

以上の点から、決定論モデルでは少数のタンパク質同士による相互作用や遺伝子 制御ネットワークを、全ての反応を詳細に表現することは極めて困難であること がわかる。
そこで決定論モデルが抱える問題を解決する方法として、確率論モデルによる表 現方法がある。この表現方法を用いると先のランダム、少数、局地的な反応に関 しては、分子同士の確率反応として捉えてやることができる。また全通りの組合 せが表現しきれない問題に関しては、時間あたりに反応する分子数を定めて、分 子が持つ反応確率に従い反応する分子を決定する方法で解決できる。

しかし、確率論モデルにも問題がある。それは反応が確率に依存し、毎回反応の 挙動が違うために系としての平均的挙動を知りたい場合には何度もシミュレーショ ンを行なわねば成らない点である。概日リズムにおいては、系の平均的な挙動が まだ不明確な点が多いためこの点は問題といえる。

そこで以上から、私が考えた手法は、概日リズム機構の系の平均的な挙動に関し てを決定論モデルに基づいた網羅的なパラメータ推定を行って表現した後、その 決定論で予測したパラメータをもとに、反応を確率論モデルで再現し、その後、 より細かな反応まで表現できるように再度パラメータ推定を行ない、詳細なタン パク質相互作用や遺伝子制御ネットワークを決定していくというものである。こ の手順に沿ってモデルを作ることで決定論モデルと確率論モデル双方が抱える問 題点を補うことができる。この決定論と確率論の両方を組み合わせたモデル表現 方法は今までにも例がなく、既存の研究にはなかった新しい研究といえる。




2 シアノバクテリア概日リズム


シアノバクテリアの概日リズムの研究は、GrobbelaarらとMitsuiらのグループに よって、シアノバクテリアの窒素固定が概日時計に制御される反応であることが 示された[10,11]ことから始まった。

そしてシアノバクテリアの概日リズムの研究は、近藤らが Synechococcus sp. PCC7942における生物発光遺伝子の移入による発光リズムの 測定実験系を作りだしたことにより[7,8,9]、急速に進歩した。

現在までにシアノバクテリアの振動に直接関わっている遺伝子として、 kaiA,kaiB,kaiC[6]の三つが同定され、また周期の安定性に関わる sasA[3],cikA[2]と概日機構の出力系に影響していると思われるcmpA[4]いった遺 伝子も同定された。しかしながら、KaiAとKaiBの結合にはKaiCが必要である [5]ということ、KaiAがpkaiBCのポジティブフィードバック因子で、KaiCがネガ ティブフィードバック因子である以外は、これらのタンパク質がどのように相互 作用を行ない、実際に概日リズムに影響を及ぼしているかということに関しては、 ほとんどわかっていないというのが現状[1]で一応、次のようなモデル図が考え られるくらいである。

シアノバクテリア概日リズムのモデル図




3 方法


遺伝子制御ネットワーク予測を以下の流れで行なう。まず最初に、遺伝子制御ネッ トワークをシミュレーションモデルとしてE-CELL System上に構築する。この際、 反応経路が不明、または反応物質が不明なところに関しては、現在わかり得る限 りの分子生物学の知見を基に何十通り、何百通りの反応モデルを作る。次に、実 験から得られる遺伝子発現のデータや生物の表現様式に基づき理論的に可能な遺 伝子ネットワークの候補を反応モデルから特定する。そして、その特定されたモ デルをシステムダイナミックス的解析を行なうことで、そのネットワーク構成を 予測するということを行なう。

そこで今回、E-CELL Systemで経路予測することを考えて、私は詳細な分子間相 互作用を表現できるスモルコフスキー方程式をベースとした、新たなシグナル伝 達/遺伝子発現制御系のモデル表現方法を開発した。この表現方法の大きな特徴は、

  • 経路最適化への応用を考慮して極めて少数のパラメータで分子間相互作用が記述できるようになっていること

  • 決定論と確率論が同一の式とパラメータで表現できるため、経路最適化を決定論で行なって式の反応パラメータを決定した後、確率論に置き換え、系に刺激を与えてやることで、野性型以外の変異型の振舞いも表現し、その結果、未知の物質や経路を探索することも可能であること

  • 分子拡散・分子の大きさなどに関して拡張できるつくりになっていること

である。決定論モデルは計算過程が確率などの不確定要素に左右されないために、 物質や反応速度の初期値が同じであれば何度シミュレーションを行なっても同じ 結果であり、その系の平均的な振舞いを表すのに向いている。一方、確率論モデ ルにおいては、反応が確率によって依存するため、系における分子同士の詳細な 振舞いを探索するのに向いていると言える。 また決定論と確率論が同一のパラメータで表現できるようにした理由としては、 まず一つ目として、概日リズム機構ではまだ系の平均的な挙動が不明な点がある ため、まず平均的な挙動を再現することが必要であること。二つ目として、平均 的な挙動を再現する決定論モデルをもとにすることで、未知経路や未知物質など を含む遺伝子制御ネットワークやシグナル伝達経路の詳細な予測を確率論で行な う際、概日リズム発生に必要不可欠な因子の影響が最小限になるため、リズムの 安定性や位相の制御を行なっている因子など、通常では予測しにくい要素の探索 に集中できることがあげられる。以上二点から決定論と確率論が同一の式とパラ メータで表現できることが重要であることがわかる。 それでは以下、その式の導出過程及び、E-CELL Systemで作成、使用している式 の説明を行なう。


3.1 スモルコフスキー式


スモルコフスキー式から一分子の反応速度が導かれるまでの説明に関しては[22]を参考にして頂きたい。

より、一分子の反応速度は、

スモルコフスキー式とは、活性化エネルギーを含め分子レベルの相互作用を確率 過程で表現した反応速度式である。自由エネルギーを用いて確率過程に基づく反 応速度式を採用しているものとしてStochSim[15,16]などがあるのに対し、活性 化エネルギーを用いているものは極めて少ない。ここで自由エネルギーについて 考える。自由エネルギーが分子同士の反応において決定するのは、その分子の反 応方向のみである。それでは、分子の反応速度はどのようにして決まるのかとい うと、反応速度とは、活性化エネルギーの山とその場が持っているエネルギー、 分子運動、そして電子的性質など複雑なパラメータにより決定されるものである。 このことから自由エネルギーではその情報量が少なく、分子の状態を正確に表し ているとは必ずしもいえない。よって活性化エネルギーを用いていることとした。




3.2 シグナル伝達と遺伝子発現系のモデル表現法


シグナル伝達系

先に述べたスモルコフスキー式を基にして、実際にE-CELL Systemで表現できる 形にしたものを以下で説明する。

という反応の場合、決定論モデルでは、

とかける。

確率論モデルでは、先に導き出した以下のスモルコフスキーの反応速度式からを導き出す。

今、Xの個数を[X]と置くと 確率

とかける。そこで反応速度 rの項を、Wと置くと、確率 は、

とかける。ここでは定数だから、これをZと置くと、

よって、は下図の活性化エネルギーであるから、

とかける。逆反応も同様にして、

(EA:活性化エネルギー、ΔG:標準自由エネルギー変化)



左図より、X→Yの場合、反応が起こるために超えなければなら ないエネルギーの山がEAなのに対して、Y→Xの場合、 EA+ΔGでなければならない。

遺伝子発現系

決定論モデルと確率論モデルは同様の式で、

を以下のように表す。




3.3 E-CELL System内部の計算方法


ここでは、具体的に確率論に基づいた場合に分子同士の反応がどのように E-CELL System内で計算されるかというアルゴリズムについて段階毎に具体的な 数値を交えて説明する。

物質が1000個未満だった場合

今、X→Y で [X]が400個、が0.2、反応をE-CELL System で計算するため積分刻み幅であるステップサイズが 1/1000秒であるとすると、

一ステップに一度反応が起こる可能性がある確率は、

[X] / stepsize = 400 /1000 = 0.4

となる。

次にこの0.4の確率で実際に反応が起こるかどうかを以下のように表現する。

if( 0.4 > ランダム関数 ){ 一回、実際に反応が起こる }

ここでいうランダム関数とは、0〜1の間で乱数を発生してくれる計算関数のこ とをいう。確率論モデルでは、このランダム関数の精度が極めて重要である。

そして反応が起こることが確定した時に、実際に分子ができるかどうかを以下のように表現する。

if( 0.2 > ランダム関数 ){ 分子が一個できる }

つまり、ランダム関数が0.4だった場合には分子はできず、0.1といった場合に分子が生成するということである。

物質が1000個以上だった場合

今、X→Yで [X]が5400個、が0.3、反応をE-CELL System で計算するため積分刻み幅であるステップサイズが 1/1000秒であるとすると、

一ステップに反応が起こる可能性がある確率は、

[X] / stepsize = 5400 /1000 = 5.4

となる。

次にこの5.4の確率で実際に反応が起こるかどうかを以下のように表現する。 この場合、5回反応が起こることは決まっているので、残りの0.4の確率が反応が 起こるどうかを見てやればよいから、

if( 0.4 > ランダム関数 ){ 6回、実際に反応が起こる }

0〜1の間の乱数を均等に発生するランダム関数が0.4より小さい時、例えば、0.2の時には6回実際に反応が起こり、

if( 0.4 =< ランダム関数 ){ 5回、実際に反応が起こる }

0〜1の間の乱数を均等に発生するランダム関数が0.4以上の時、例えば、0.6の時には 5回実際に反応が起こるというように計算される。

そして反応が起こることが確定した時、実際に分子ができるかどうかを反応の起こる回数分だけ行なう。具体的にいうと、

if( 0.3 > ランダム関数 ){ 分子が一個できる }

例えば、2において6回実際に反応が起こると決まった場合には、上式の「分子が一個できるのか」という同じ試行を6回行なってやることで、一分子一分子に対して分子が生成するかどうかの判定を確率的論に基づき行なっている。

以上の過程を経て、確率論に基づく分子生成が起こる。




3.4 サンプルデータを用いた決定論と確率論モデルの比較


これまで説明してきたスモルコフスキー式をベースとした手法で決定論と確率論 が、E-CELL Systemにおいてどのように表現されるのかをToy Pathwayを用いて以 下に示した。Toy Pathwayは、下図にもあるようにAからB、Cを経てDへと流れる反応 である。



決定論モデルでToy Pathwayを100秒走らせたシミュレーション結果全体図




確率論モデルでToy Pathwayを100秒走らせたシミュレーション結果拡大図




決定論モデルでToy Pathwayを100秒走らせたシミュレーション結果拡大図




決定論モデルでToy Pathwayを100秒走らせたシミュレーション結果全体図





4 パラメータ推定


パラメータ推定手法に遺伝的アルゴリズムを用いた理由は、まず一つ目として 遺伝子制御ネットワークの予測を行なう際、未確定のパラメータがかなり多く、 パラメータの探索域も非常に広く、またどこに目的値があるかの検討すらつかな い。よって網羅的なパラメータ予測を行なう必要がある。しかしその際、他の手 法では局所解に陥ってしまう危険があり、初期の優秀個体のために目的値に収束 しない可能性がでてくる。この点、遺伝的アルゴリズムはユーザーが目的に応じ て設定を変えることができる。

次に、最適化するための基となるデータは、実験において反応活性データを測定 することが非常に難しいことから反応活性データではなく、測定可能なタンパク 質分子の数である。しかし我々が予測したいのはE-CELL Systemにおけるリアク ター内で定義される反応活性を表すパラメータであり、このパラメータの値が系 の挙動に与える影響が大きいか小さいかはわからない。つまり、反応自体を直接 チューニングするのではなく、反応した物質を間接的にチューニングするという ことになる。このため、E-CELL Systemにおいて最適解を見つけることは非常に 難しい。

以上のに二点より、今回、パラメータ推定の手法として遺伝的アルゴリズムを用 いることにした。遺伝的アルゴリズムに関しての詳細は[21]を参考にして頂きた い。




5 結果


現在のところ目的の解を得るまでには至っていない。以下に9回に及ぶチューニング結果を記す。





Case1においてGA終了後、得られた最適解を4800秒シミュレーションした結果



Case3においてGA終了後、得られた最適解を4800秒シミュレーションした結果



Case4においてGA終了後、得られた最適解を4800秒シミュレーションした結果



Case5においてGA終了後、得られた最適解を4800秒シミュレーションした結果



Case6においてGA終了後、得られた最適解を4800秒シミュレーションした結果



Case1218においてGA終了後、得られた最適解を4800秒シミュレーションした結果




6 考察


目的解が得られなかった数ある原因のうちの一つ評価関数に関する問題について 取り上げたい。現在、E-CELL Systemにおけるパラメータ推定機構における評価 関数は以下のようにしている。

:E-CELL での計算値, :実験から得られる時系列データ、n:測定できる変数の数, N:サンプリングポイント数

上記の式にあるように、計算で得られた変数と実験で得られた変数の二乗差を最小にするという方法、最小二乗法を用いてい る。 実際にパラメータ推定を行なってみると、KaiAとKaiCの値が一定値に収束して定 常状態になる場合を最適解としてしまう場合がほとんどであった。 初期値の段階で、目的解の周辺をパラメータ推定機構が探索していれば現在の評 価関数でも解は導ける。しかし、探索範囲が非常に広い今回のような場合では、 いきなり理想解の周辺をパラメータ推定機構が探索するという可能性は極めて少 なく、その結果、一定値に収束して定常状態になったものを最適解としてしまっ たと考えられる。よって今回用いた評価関数では、振動を検出する方法としての効力 が弱く、理想解へつながる解が例えあったとしても、最適解として認識できなかっ た可能性が高い。


7 展望





7.1 計算速度の向上


現在のE-CELL System(E-CELL Version 1.0)において未知変数が20程度の場合約 500世代のチューニングにかかある時間は約2週間程度であり、この計算速度では 仮想実験空間としての機能を果たさない。しかしながら、同プロジェクトの高橋 らを中心とする次世代E-CELL System(E-CELL Version 3.0)の開発が進み、本格 的なシミュレータとしての機能を果たすまでの期間が一年以上ある点を考えると、 現在のE-CELL Systemを大規模な改良を伴わないで高速化することは必須である といえる。 そこで考えられている方法がいくつかあるが、まず一つ目としてコンパイラを現 在使っているgccから、cxx(Compaq C++ Compiler)に変えることで約20\%程度の 高速化が見込まれる。ただし現在STL関係でソース自体に変更が必要なため、同 プロジェクトのE-CELL System設計・開発者である高橋恒一氏らと協力しながら 実現したい。またパラメータチューニング手法の改良として、gprofというツー ルを使ってプログラム中でどこがボトルネックになっているかを調べることも行 ないたい。これにより、ボトルネックとなっている部分を集中的にチューニング することでかなり計算速度が向上すると期待できる。


7.2 モデルルール生成環境の構築


現在、たかだが数種類程度による物質による遺伝子発現と相互作用をモデルとし て表現しているが、今後、反応系に関与している物質が増えていく。それに伴っ て記述しなければならないルールファイルもますます増えていくため、網羅的に ルールファイルを作成する環境を整える必要がある。具体的な方法としては、 E-CELL Systemの周辺モジュールの一つとして、必要な物質を与えてやることで 網羅的にルールファイルを記述してくれるツールの構築を考えている。


謝辞


数々のアドバイスを頂き、親身になって指導して頂いた慶応義塾大学環境情報 学部助手の中山洋一氏、パラメータ推定に関する有益なアドバイスをして頂い た齋藤裕介氏、専門的立場から意見を頂いた理化学研究所の齋藤輪太郎氏、奈 良先端技術大学院大学の森浩禎氏に心から感謝したい。そして総合的に支援・ 指導して頂いた慶応義塾大学環境情報学部教授の冨田勝氏にこの場を借りて心 から感謝の意を表したい。


参考文献


  • Iwasaki H, Kondo T. (2000) Plant Cell Physiol, 41(9), 1013-20.
  • Schmitz O, Katayama M, Williams SB, Kondo T, Golden SS. (2000) Science, 289(5480), 765-8.
  • Iwasaki H, Williams SB, Kitayama Y, Ishiura M, Golden SS, Kondo T. (2000) Cell, 101(2), 223-33.
  • Katayama M, Tsinoremas NF, Kondo T, Golden SS. (1999) J Bacteriol, 181(11), 3516-24.
  • Iwasaki H, Taniguchi Y, Ishiura M, Kondo T. (1999) EMBO J, 18(5), 1137-45.
  • Ishiura M et al. (1998) Science, 281(5382), 1519-23.
  • Kondo T et al. (1994) Science, 266(5188), 1233-6.
  • Kondo T, Ishiura M. (1994) J Bacteriol, 176(7), 1881-5.
  • Kondo T et al. (1993) Proc Natl Acad Sci U S A, 90(12), 5672-6.
  • Grobbelaar N, Huang T-C, Lin H-Y, Chow T-J. (1986) FEMS Microbiol.Let, 37, 173-177
  • Mitsui A, Kumazawa S, Takahashi A, Ikemoto H, Cao S, Arai T. (1986) Nature, 323, 720-722
  • Goldbeter A. (1995) Proc R Soc Lond B Biol Sci, 261(1362), 319-24.
  • Leloup JC, Goldbeter A. (1998) J Biol Rhythms, 13(1), 70-87.
  • Leloup JC, Goldbeter A. (1999) J Theor Biol, 198(3), 445-59.
  • Morton-Firth CJ, Shimizu TS, Bray D. (1999) J Mol Biol, 286(4), 1059-74.
  • Morton-Firth CJ, Bray D. (1998) J Theor Biol, 192(1), 117-28.
  • Akutsu T, Miyano S, Kuhara S. (1999) Pac Symp Biocomput, 17-28.
  • Ideker TE, Thorsson V, Karp RM. (2000) Pac Symp Biocomput, 305-16.
  • Savageau MA. (1976) Biochemical System Analysis : A Study of Function and Design in Molecular Biology, Addison-Wesley, Massachusetts.
  • Tominaga D et al. (1999) Computer Science and Biology (Proc of the German Conf on Bioinformatics(GCB'99)), 127.
  • 齋藤裕介. (1999) 生命と情報, 6, 329-339.
  • 北原和夫. (1994) 非平衡系の科学II 緩和過程の統計力学, 109-114.