Quantitative Analysis for Evolution of Mutational Robustness and Bistability: Use of a Rare-Event Sampling Method

Abstract

生命は進化の過程で機能に加えて変異に対する頑健性を獲得した．この過程を数理モデルで理解するために，進化シミュレーションと比較する対照群をマルチカノニカル法で構築する手法を提案する．遺伝子制御ネットワークで計算を行い，進化が頑健性を増強させることと進化では双安定性の出現が遅れることを定量的に示した．

1.  進化と変異に対する頑健性

生命は進化によって機能を獲得してきた．この進化の過程は遺伝子型に生じる変異とその結果現れた表現型を持つ個体に対する淘汰の繰り返しであるが，遺伝子の変異によって容易に機能を失ってしまうような遺伝子型の系列は進化の過程で淘汰されてしまう．そのため，遺伝子の変異に対する頑健性を備えた遺伝子型が子孫を残しやすい．これは機能とは無関係に起きる淘汰なので，second-order selectionと呼ばれることがあり，進化という過程に特有の現象である¹⁾．実際，変異に対する頑健性は遺伝子の網羅的ノックアウトなどで実験的に確認できる．しかし，頑健性がどのように進化するかを実験的に調べるのは難しい．そこで我々は数理モデルを用いてこの問題に取り組んでいる．

数理モデルでは便宜上なんらかの適応度関数を設定して，それが大きくなるように進化させることが多い．要するに最適化なのだが，そこには進化ならではの現象があり，その代表が上述の変異に対する頑健性である．しかし，進化シミュレーションを行うだけでは，進化の特殊性を知ることはできない．我々は進化という過程そのものが持つ特殊性を数値計算で理解したい．そのためには進化シミュレーションと比較する対象がほしい．以下では遺伝子制御ネットワーク（GRN）のモデルを考えるが，その場合に適切な比較相手はランダムに作ったGRNの集合だろう．といっても，単純なランダムサンプリングは役に立たない．というのは，殆ど全てのランダムなGRNは機能を持たず，適応度が高いGRNは稀だからである．この問題を解決するために，我々はマルチカノニカル・モンテカルロ（McMC）法を用いた研究手法を提案している．

本稿の主眼は手法の紹介なので，まずMcMC法について解説する．生物学的なネットワークのモデルに対して初めてMcMC法を用いた研究は斉藤と筆者による²⁾．本稿と直接関係するのは永田と筆者³⁾および金子と筆者⁴⁾の研究であるが，本稿の後半では文献4の結果を紹介する．

2.  対照群を作る：マルチカノニカル法

McMC法はマルコフ連鎖モンテカルロ法（メトロポリス法）の一種である．ここではカノニカル分布を作り出すためのメトロポリス法は既知として話を進める．実際に我々が用いたのはMcMC法のバリエーションであるエントロピック・サンプリングだが，本稿ではまとめてMcMC法と呼ぶ．

McMC法はもともと統計物理学の分野で広いエネルギー範囲の微視的状態を均等にサンプリングするために提唱された⁵⁾．一次相転移を考えてみよう．シミュレーションで扱う有限系の場合には，転移点付近で高温相と低温相が自由エネルギー障壁で隔てられて共存するはずである．しかし，自由エネルギー障壁にあたる部分は状態数が少ないため，メトロポリス法では越えるのが難しく，共存するはずの二相の一方しか実現しないことが多い．そこで，なんとかして状態数が少ないエネルギーを重点的にサンプリングできれば，両相を頻繁に行き来させられるだろう．

そのために，McMC法ではエネルギーEを持つ微視的状態の重みをボルツマン重率ではなくexp(–S(E))とする．S(E)はエントロピーで，状態数W(E)とS(E) = log W(E)の関係にある．つまり，状態数の逆数に比例する重みにする．これなら，全てのエネルギーは等確率で出現するはずである．よく誤解されるが，ランダムサンプリングではないことに注意してほしい．

といっても，状態数があらかじめ分かっているなら，分配関数が計算できてしまうのだからシミュレーションの必要はないだろう．McMCの最大の発明は，重みが大雑把にexp(–S(E))に近ければ実用上は十分だと見抜いた点にある．少しくらい真のS(E)からずれていても，全てのエネルギーが十分にたくさんサンプリングできていさえすればいい．そこでシミュレーションを行う前にS(E)を「学習」によって事前に推定しておく．これをσ(E)としよう．シミュレーションに必要なのは異なるエネルギー間でのエントロピー差だけなので，σ(E)にも付加定数だけの不定性がある．重みをexp(–σ(E))としてメトロポリス法を行い，得られたエネルギーヒストグラムをH(E)としよう．すると，

  

$$\mathit{H}\left(E\right)\propto W\left(E\right)\text{exp}\left(-\mathrm{\sigma }\left(E\right)\right)$$

の関係があるから，逆に解いて状態数は

  

$$\mathit{W}\left(E\right)\propto H\left(E\right)\text{exp}\left(\mathrm{\sigma }\left(E\right)\right)$$

と推定できる．あるいは

  

$$\mathit{S}\left(E\right)=\mathrm{\sigma }\left(E\right)+\text{log}H\left(E\right)$$

によって，エントロピーの推定値が（付加定数を除いて）得られる．状態数が得られてしまえば，そこから有限温度のカノニカル分布を構成するのはやさしいが（reweightingと呼ばれる），今回の研究では必要ない．

この手法では，サンプリング間隔を十分に空けさえすれば，エネルギーの値ごとにそのエネルギーを持つ状態をランダムサンプリングできる．一般にはサンプル間に相関が生じるが，サンプルを十分に多くとれば問題にならない．なお，ランダムな初期状態から出発できる場合には初期緩和が必要ないことを指摘しておく（高温極限でのメトロポリス法に初期緩和が必要ないのと同様である）．

McMC法が提唱されると，適切な「エネルギー関数」さえ設定すれば非物理系にも適用できることが認識されるようになった⁶⁾．どんな対象であれ，稀な状態を効率よくサンプリングする手法として使えるわけである．特に，（1）原理的にランダムサンプリングできる（2）全状態数が分かっている問題では特定の性質を持つ稀な状態の数が計算できる，の二点が重要である．例えば，北島と筆者はMcMC法を用いて大きな魔方陣の数を推定した⁷⁾．

McMC法を進化の数理モデルに用いるには，遺伝子型を微視的状態とし，エネルギー関数として適応度を採用すればいい．これにより，広い適応度範囲に渡って各適応度を持つ遺伝子型をランダムにサンプリングできる．こうして得られた遺伝子型の集合を進化シミュレーションと比較するための対照群とする．

実際の計算ではエネルギーを有限幅の区間に分割し，各区間内にある状態の出現確率が全ての区間でほぼ等しくなるようにσ(E)を学習させる．特に各区間の中ではσ(E)を一定値とするのがエントロピック・サンプリングである⁸⁾．学習の手法としてはいろいろなものが考えられる．極端なことを言えば，勘で決めても構わないのだが，エントロピック・サンプリングのための学習法としてはWang-Landau（WL）法が簡便で効率もよい．WL法はもともとMcMC法とは独立にW(E)を推定する手法として提案された⁹⁾．しかし，詳細釣り合いを満たさない手法のため，単独で用いるのはあまり筋がよくなく，McMC法の学習手法としてのみ使うべきである．WL法の詳細については原論文に譲る．また，WL法だけで重みを決めるのではなく，WL法によって得られたσ(E)を使って一度McMC計算を行い，そのヒストグラムからσ(E)を更新してそれを本番のMcMC計算に用いることを勧める．このあたりの細かいことはMcMC法を使っている研究者なら知っているはずだが，あらわに書いた文献は意外にない．ちなみに，WL法を含め，学習の部分を並列化する手法は確立していないが，McMC法は独立な計算を並列に走らせることができるのでサンプル数を稼ぎやすい．

まとめると計算は以下の三段階からなる．

1．WL法によって最初のσ(E)を作る．

2．適当な長さのMcMC計算を行いσ(E)を改良する．

3．得られたσ(E)を使って本番のMcMC計算を行う．

3.  遺伝子制御ネットワークのモデル

我々はGRNのモデルとしてニューラルネットワーク型のモデルを用いる．遺伝子発現の詳細は忘れ，その制御関係だけを有向グラフで表したモデルである．7ノードからなる小さなグラフの例を図1に示す．ノードが遺伝子でエッジが制御関係を表している．外界からの入力を受け取るノードと目的のタンパク質を発現するノードをひとつずつ持つ1入力1出力のモデルで話を進めよう．

図1

7ノードのGRNの例．Iは入力ノード，Oは出力ノードを表す．また，赤線（矢印）は活性化，青線（T印）は抑制を表す．

ここでは与えられたネットワークに対して以下の離散ダイナミクスを仮定する．

  

$${\displaystyle x_i\left(t+1\right)=R\left(I{\mathrm{\delta }}_{i0}+\sum _{j\ne i}{J}_{\mathit{\text{ij}}}{x}_j\left(t\right)\right)}$$

  

$${\displaystyle R\left(x\right)=& \frac{1}{1+\text{exp}\left(-\mathrm{\beta }\left(x-\mathrm{\mu }\right)\right)}}$$

x_i(t)はi番目の遺伝子のtステップ目での発現量，J_ijはj番目の遺伝子からi番目の遺伝子への制御を表し，今回は+1が活性化，–1が抑制，0が制御関係なしのみっつの値を取れるものとする．入力遺伝子（0番）にだけ入力 $I\in \left[\mathrm{0,1}\right]$ が加えられる．Rはシグモイド型の応答関数で，以下でμ = 0，β = 2とする．

このGRNに対し，入力がI = 0とI = 1の場合をなるべく鋭敏に見分けることを要請し，適応度関数fを以下のように定義する．

  

$$f=\left|\stackrel{\_}{x_{\mathrm{\text{out}}}}\left(1\right)-\stackrel{\_}{x_{\mathrm{\text{out}}}}\left(0\right)\right|$$

$\stackrel{\_}{x_{\mathrm{\text{out}}}}\left(I\right)$ は入力Iのもとで定常状態に達した後の出力ノードの発現量を表す．McMC法の毎ステップで，適応度を求めるために定常状態に達するまでダイナミクスを二度走らせなくてはならないことに注意しよう．

以下の節ではノード数を32，エッジ数を80に固定している．進化シミュレーションでは点変異だけを考える．つまり，ランダムに選んだエッジをひとつ取り除き，ランダムに選んだノードペア間に新たなエッジを付け足す．個体数は1000として，適応度上位の500個体を残し，そのコピーに点変異をほどこす．文献4では900個体残す場合も計算しており，結果の大筋は変わらない．

4.  遺伝子型エントロピー

図2はMcMCで得られた適応度のヒストグラムと計算に使用した重みから遺伝子型のエントロピーを求めたものである． $\left[\mathrm{0,1}\right]$ の間を100区間に分割し，各区間で平均50000サンプルのGRNが得られるように設定した．適応度最大の区間の実現確率は適応度最小の区間と17桁ほど違っており，確かに稀なGRNが十分にたくさんサンプルされている．進化シミュレーションでの適応度の変化も同じ図に重ねてある．これを見ると，エントロピーが急激に下がりはじめる適応度で進化の速度も急激に遅くなっている．これは進化速度が大まかにはエントロピーで決まることを示唆している．

図2

（青破線）McMCで得られたヒストグラムから求めたエントロピー．縦軸（左）は状態数の常用対数を表している．適応度を100区間に分割して計算した．（オレンジ実線）進化シミュレーションによる適応度の変化．縦軸（右）は下に向かって世代が進む．1000個体中で150世代目に最も適応度が高かった個体の系統を辿り，それを10万系統について平均した．縦の点線はエントロピーが急減する適応度の目安である（文献4より改変）．

5.  変異に対する頑健性と双安定の進化

変異に対する頑健性の指標として，ここではGRNからエッジをひとつ取り除いた時の適応度を全ての可能な取り除き方について平均したものを用いよう．この量は適応度とともに増えるので，異なる適応度間での比較には意味がないが，McMC法と進化シミュレーションを同じ適応度で比較する目的には十分である．

図3Aは頑健性の指標を適応度に対して描いたものである．適応度が低いうちは両者が一致しているのに対し，適応度が高くなると進化シミュレーションの頑健性がMcMC法で得られたものより有意に大きくなっている．これより，進化が変異に対する頑健性を増強させることが定量的に確かめられた．そこで， $f\in \left[\mathrm{0.99,1}\right]$ の区間で不可欠なエッジ，すなわちそれを切ると適応度が0に近くなるようなエッジの数の分布を調べてみると，McMC法と進化シミュレーションでは分布が大きく異なり，進化では不可欠なエッジが少ない（図3B）．これより，不可欠なエッジを多く持つ遺伝子型が進化の過程で淘汰されたことが分かる．

図3

（A） 頑健性の指標の適応度依存性．進化シミュレーションで得られたGRNもMcMC法と同様に適応度で100区間に分類した．（B） $f\in \left[\mathrm{0.99,1}\right]$ での不可欠なエッジの数の分布．不可欠なエッジとは，削除すると適応度がほぼ0になるエッジを指す（文献4より改変）．

ここで用いたパラメーター値では非常に高い適応度でほぼ全てのGRNが双安定性を示す．つまり，入力を0から1まで少しずつ変えた時と逆過程とで出力ノードの発現量変化にヒステリシスが見られる．これはMcMC法で確かめられるので，進化経路に依存しない性質である．図4は双安定GRNの割合を適応度に対して描いたものである．進化でも最終的には全てのGRNが双安定になるが，双安定GRNの増加のしかたはMcMC法と進化で異なっており，同じ適応度でも進化のほうが双安定GRNの比率が低い．言い換えると，進化では双安定GRNが出現しづらい．双安定性を新たな表現型とみなすなら，新たな表現型の出現が進化では遅れるわけである．

図4

各適応度における全GRN中で双安定性を示すGRNの比率．McMCと進化シミュレーションのいずれもf→1で双安定GRNの比率が1に漸近する（文献4より改変）．

この理由を調べるためにGRNを単安定と双安定に分け，それぞれについて適応度と頑健性の分布をMcMCで求めたのが図5である．両者の分布は大きく異なっており，単安定GRNでは頑健性の指標が適応度とともに上がっていくのに対し，双安定GRNでは高い適応度でも頑健性が低いものが多く見られる．進化の過程では変異に対する頑健性が低い遺伝子型は淘汰されやすいので，結果的に双安定GRNは単安定GRNに比べて出現しづらくなると考えられる．つまりこの結果は，同じ適応度であっても変異に対する頑健性の違いによって特定の表現型が選択されづらくなる場合があることを示唆している．

図5

McMCで得られたGRNを双安定と単安定に分類し，それぞれについて適応度と頑健性の指標の分布を描いたもの．適応度の各区間でサンプル数の総和は約5万であり，カラーバーの単位は万である（文献4より改変）．

6.  おわりに

進化という過程の特性を調べるためにMcMC法で対照群を構築して進化シミュレーションと比較する手法を紹介した．GRNでの計算結果は，我々の手法が進化について非自明な結論を導くのに有効であることを示している．この手法はGRN以外の系の進化にも適用できる．ただ，進化シミュレーションに比べるとかなり重い計算なので，計算量は問題ごとに検討する必要がある．また，進化に限らず，機械学習の学習過程なども同じ手法で議論できることを指摘しておきたい．ぜひ様々な問題で試みていただきたい．

本稿では適応度関数が一次元的な場合だけを示したが，現実の進化では適応度関数をこれほど単純には定義できない．試みに我々はふたつの適応度関数を持つ系についても計算を行っている．この場合，二次元の適応度空間に対してMcMCを行うことになる．このような拡張では計算量が次元とともに急激に増えるため，現状では二次元止まり，あるいは非常に大規模な計算機を用いるとしても三次元程度が限度だろう．

謝辞

本稿で紹介した研究は大学院学生だった永田新太郎氏と金子忠宗氏との共同研究に基づいている．またMcMC法については斉藤稔氏と伊庭幸人氏との議論から得たものが多い．各氏に感謝する．

文献

• 1)  Wagner, A. (2005) Robustness and Evolvability in Living Systems, Princeton Univ Press.
• 2)   Saito,  N.,  Kikuchi,  M. (2013) New J. Phys. 15, 053037. DOI: 10.1088/1367-2630/15/5/053037.
• 3)   Nagata,  S.,  Kikuchi,  M. (2020) PLoS Comput. Biol. 16, e1007969. DOI: 10.1371/journal.pcbi.1007969.
• 4)   Kaneko,  T.,  Kikuchi,  M. (2022) PLoS Comput. Biol. 18, e1009796. DOI: 10.1371/journal.pcbi.1009796.
• 5)   Berg,  B. A.,  Neuhaus,  T. (1991) Phys. Lett. B 267, 249-253. DOI: 10.1016/0370-2693(91)91256-U.
• 6)   Iba,  Y. et al. (2014) Ann. Inst. Stat. Math. 66, 611-645. DOI: 10.1007/s10463-014-0460-2.
• 7)   Kitajima,  A.,  Kikuchi,  M. (2015) PLoS One 10, e0125062. DOI: 10.1371/journal.pone.0125062.
• 8)   Lee,  J. (1993) Phys. Rev. Lett. 71, 211-214. DOI: 10.1103/PhysRevLett.71.211.
• 9)   Wang,  F.,  Landau,  D. P. (2001) Phys. Rev. Lett. 86, 2050-2053. DOI: 10.1103/PhysRevLett.86.2050.

Biographies

菊池　誠（きくち　まこと）

大阪大学サイバーメディアセンター教授