- 元論文: Survival of Self-Replicating Molecules under Transient Compartmentalization with Natural Selection
自己複製分子 (replicator) が、寄生体（parasite）を伴った環境において一時的な区画化（transient compartment）の下で自然選択的に生存するか、生存するならば一時的な区画化は生存に対してどのような影響を与えるのかを調べた研究。
実験方法
次の 4 ステップを 1 round として、$R (\in \mathbb{N})$ round 実行する。
-
プールから区画に分子を割り当てる（Inoculate）
- 区画の総分子数 $n = m + y \in \mathbb{N}$ は平均 $\lambda$ のポアソン分布に従うものとみなしてサンプリングする。ここで $m$ は replicator の数で $y$ は parasite の数を表す。
- replicator の数 $m$ はパラメータ $x$ を持つ二項分布に従うものとみなしてサンプリングする。
-
各区画で replicator と parasite の時間発展を行う（Maturate）
- Inoculate ステップ で得た $m$ と $y=n-m$ を初期値として、次のダイナミクスを実行する。
$$ \begin{aligned} \dot{m}(t) &= \alpha m(t)^2\\ \dot{y}(t) &= \gamma m(t) y(t) \end{aligned} $$
- 上のダイナミクスを見れば分かるように $m$ も $y$ もひたすら増え続けるため、あるタイミング（これを停止時間 $T$ としている）で止める必要がある。元論文では、この停止条件は $m(T) + y(T) = K + n$ によって定められている。
- $K$ は、区画内に存在する栄養資源の制約下において複製過程で新しく作ることが可能な replicator の鎖の数（環境収容力）を表す。
-
区画の分子をプールに移す（Pool）
- プール内の replicator の割合 $x’(x, \lambda)$ は次の式で表される（詳細は後述する）。
$$ \begin{aligned} x’(x, \lambda) = \frac{\langle m \rangle}{\langle n \rangle} = \frac{\lambda x + K e^{-\lambda} (e^{\lambda x} -1)}{\lambda + K(1 - e^{-\lambda x})} \end{aligned} $$
- ここで $\langle \cdot \rangle$ の期待値操作は、Inoculate ステップのポアソン分布と二項分布のサンプリングによる確率密度関数 $P_\lambda(n, m, x)$ によるものであることを明記しておく（自分が少し混乱したので）。
- また、この $x’(x, \lambda)$ を次の round の $x$ として使うことにも注意する。
-
プールの中の分子を希釈する（Dilution）
- 次の round の Inocurate ステップに使うパラメータ $\lambda$ は、希釈係数 $d > 1$ を用いて次式で更新される。
$$ \lambda = \frac{\langle n \rangle}{d} $$
式の導出
元論文には $P_\lambda(n, m, x)$ がどういう確率密度関数なのか明示的に書かれていないが、これは Selection Dynamics in Transient Compartmentalization に書かれており、
$$ \begin{aligned} P_\lambda(n, m, x) &= \mathrm{Poisson}(n, \lambda) \times \mathrm{Binomial}(n,m,x) \\ &= \frac{\lambda^n}{n!}e^{-\lambda} {}_n C_m x^m (1-x)^{n-m} \end{aligned} $$
である。この確率密度関数を用いて $\langle \bar{n} \rangle, \langle \bar{m} \rangle$ を計算する。
$\langle \bar{n} \rangle$ の計算
元論文の (A4) 式より、 $$ \begin{aligned} \langle \bar{n} \rangle &= \sum_{n > 0} n P_{\lambda}(n, 0, x) + \sum_{n>0} \sum_{m=1}^{n} (K+n) P_{\lambda}(n, m, x). \end{aligned} $$
計算の見通しを良くするために、$y=1-x$ と置いておく。
右辺の第一項は、 $$ \begin{aligned} \sum_{n>0}{n P_{\lambda}(n, 0, x)} &= \sum_{n>0}{ n \frac{\lambda^n}{n!} e^{-\lambda} {}_n C_0 x^{0} y^{n-0}} \\ &= \sum_{n>0} { \frac{(\lambda y)^{n-1}}{(n-1)!} e^{-\lambda y} e^{\lambda y} e^{-\lambda} \lambda y} \\ &= e^{\lambda y} e^{-\lambda }\lambda y. \end{aligned} $$
右辺の第二項は、 $$ \begin{aligned} \sum_{n>0} \sum_{m=1}^{n} (K+n) P_{\lambda}(n, m, x) &= \sum_{n>0} (K+n) \frac{\lambda^n}{n!}e^{-\lambda} \sum_{m=1}^{n} {}_n C_m x^m y^{n-m}\\ &= \sum_{n>0} (K+n)\frac{\lambda^n}{n!}e^{-\lambda} \left( 1- {}_n C_0 x^0 y^{n-0} \right) \\ &= \sum_{n>0} (K-Ky^n+n-ny^n) \frac{\lambda^n}{n!} e^{-\lambda} \\ &= K\left(1 - \frac{\lambda^0}{0!} e^{-\lambda} \right) - K \sum_{n>0}\frac{(\lambda y)^n}{n!}e^{-\lambda y}e^{\lambda y} e^{-\lambda} + \sum_{n>0} \frac{\lambda^{n-1}}{(n-1)!} e^{-\lambda} \lambda - \sum_{n>0} \frac{(\lambda y)^{n-1}}{(n-1)!} e^{-\lambda y} e^{\lambda y} e^{-\lambda }\lambda y \\ &= K(1-e^{-\lambda}) - Ke^{\lambda y} e^{-\lambda }(1-e^{-\lambda y}) + \lambda - e^{\lambda y} e^{-\lambda } \lambda y \\ &= K - Ke^{-\lambda} e^{\lambda y} + \lambda - e^{\lambda y} e^{-\lambda} \lambda y. \end{aligned} $$
先程計算した第一項と第二項を足すと、 $$ \begin{aligned} \langle \bar{n} \rangle &= e^{\lambda y} e^{-\lambda }\lambda y + K - Ke^{-\lambda} e^{\lambda y} + \lambda - e^{\lambda y} e^{-\lambda} \lambda y \\ &= \lambda + K(1 - e^{-\lambda}e^{\lambda (1-x)}) \\ &= \lambda + K(1 - e^{-\lambda x}) \end{aligned} $$ となり、元論文の結果と一致することが確認できた。
$\langle \bar{m} \rangle$ の計算
元論文の (A7) 式より、 $$ \langle \bar{m} \rangle = \sum_{n>0} (K+n) P_{\lambda}(n, n, x) + \sum_{n > 0} \sum_{m=1}^{n-1} m P_{\lambda}(n, m, x). $$
右辺の第一項は、 $$ \begin{aligned} \sum_{n>0} (K+n) P_{\lambda}(n, n, x) &= \sum_{n>0} (K+n) \frac{\lambda^n}{n!} e^{-\lambda} {}_n C_n x^{n} y^{n-n} \\ &= \sum_{n>0} (K+n) \frac{(\lambda x)^n}{n!} e^{-\lambda x} e^{\lambda x} e^{-\lambda}\\ &= Ke^{\lambda x}e^{-\lambda}(1-e^{-\lambda x}) + e^{\lambda x}e^{-\lambda} \lambda x. \end{aligned} $$
右辺の第二項は、 $$ \begin{aligned} \sum_{n > 0} \sum_{m=1}^{n-1} m P_{\lambda}(n, m, x) &= \sum_{n > 0} \frac{\lambda^n}{n!}e^{-\lambda} \sum_{m=1}^{n-1} m \times {}_n C_m x^m y^{n-m} \\ &= \sum_{n > 0} \frac{\lambda^n}{n!}e^{-\lambda} \sum_{m=1}^{n-1} n \times {}_{n-1} C_{m-1} x^{m-1} x y^{n-m}\\ &= \sum_{n > 0} \frac{\lambda^n}{(n-1)!} e^{-\lambda} x \left( \sum_{m=1}^{n-1} [ {}_{n-1} C_{m-1} x^{m-1} y^{n-m} ] + x^{n-1} - x^{n-1}\right)\\ &= \sum_{n > 0} \frac{\lambda^n}{(n-1)!} e^{-\lambda} x \left( \sum_{m=1}^{n} [ {}_{n-1} C_{m-1} x^{m-1} y^{n-m} ] - x^{n-1}\right) \\ &= \sum_{n > 0} \frac{\lambda^n}{(n-1)!} e^{-\lambda} \left\{ x (x+y)^{n-1} - x^{n}\right\}\\ &= \sum_{n>0} \frac{\lambda^{n-1}}{(n-1)!} e^{-\lambda} \lambda x - \sum_{n>0} \frac{(\lambda x)^{n-1}}{(n-1)!} e^{-\lambda x} e^{\lambda x} e^{-\lambda} \lambda x \\ &= \lambda x - e^{\lambda x} e^{-\lambda} \lambda x. \end{aligned} $$ 一行目から二行目は、 $$ \begin{aligned} m \times {}_n C_{m} &= m \frac{n!}{m!(n-m)!} \\ &= n \frac{(n-1)!}{(m-1)!(n-m)!} \\ &= n \times {}_{n-1} C_{m-1} \end{aligned} $$ を使った。 また、四行目から五行目は二項定理を用いた。
第一項と第二項を足すと、 $$ \begin{aligned} \langle \bar{m} \rangle &= Ke^{\lambda x}e^{-\lambda}(1-e^{-\lambda x}) + e^{\lambda x}e^{-\lambda} \lambda x + \lambda x - e^{\lambda x} e^{-\lambda} \lambda x \\ &= Ke^{-\lambda}(e^{\lambda x} - 1) \lambda x + \lambda x \end{aligned} $$ となり、元論文の結果と一致することが確認できた。
シミュレーションの再現実装
基本的には実験方法の通りに実装すれば良い。
replicator が振動しながら生存 ($d = 18$)
実験に用いたパラメータを次の表に示す。
| パラメータ | 値 |
|---|---|
| round 数 $R$ | 200 |
| 区画数 $C$ | 5000 |
| 環境収容力 $K$ | 60 |
| replicator の増殖係数 $\alpha$ | 0.29 |
| parasite の増殖係数 $\gamma$ | 1.5 |
| 希釈係数 $d$ | 18 |
| 区画の総分子数 $n$ の平均を表すパラメータ $\lambda$ の初期値 | 5 |
| 区画の総分子数 $n$ から replicator が取られる確率 $x$ の初期値 | 0.5 |
まず最初に、初期値を変えたときの Maturate ステップの時間発展の挙動を見ておこう。
- replicator が最初から存在しない場合、 replicator が増殖しないのはもちろんのこと、parasite も時間変化量が 0 になるので初期値から変わらない。
- replicator と parasite の数が等しい場合、 $ \gamma > \alpha$ であり、parasite の数 $y$ は $m$ に影響を与えないので、ひたすら $y$ が増殖する。
- replicator の数を parasite よりも多く取った場合 ( $m=30, y=1$ )、replicator と parasite の数は良い感じに増殖するが、環境収容力 $K + n$ ($n = 31$ とした) の制約によって実行は停止される。
- ところでポアソン分布によって $n$ をサンプリングしているため、今回の実験パラメータにおいては $\lambda$ がよほど大きくならない限り $n=31$ は殆ど取られない数ということには注意しておこう。
次に、 200 round 実行して得られた $\langle m \rangle, \langle y \rangle, \lambda, x $ の時系列をプロットしてみよう。
良い感じに振動現象が見られたことが確認できた。 振幅が多少揺らいでいるため、論文のような安定した振動は得られなかった。 これは、区画数が $5000$ と少ないことに起因していると考えられる。
論文には無いが、ついでに相図もプロットした。ロトカ・ボルテラのようなリミットサイクルらしきものが見られて嬉しい。 これを見て、一時的な区画化というのは、区画数 $C \rightarrow \infty$ で実はロトカ・ボルテラのような単純な微分方程式に落とし込むことができるのかもしれないな、と思った。 ただし $C\rightarrow \infty$ というのはあまり現実的では無いので、これを考えることに意味があるかどうかは不明（ただ個人的には純粋に面白いとは思う）。
希釈係数 $d = 19$ の結果
パラメータ設定は $d$ 以外は前と同じ。 元論文では $d=19$ において綺麗な振動現象が見られたので、それを再現するかどうかの確認が目的。
先程と同じように、時系列データと相図を貼っておく。
もう少し round 数を増やせば実はリミットサイクルが小さいながら存在しそうなことが確認できそうな気もするが、なんだか減衰振動っぽいものが見られて、元論文の結果を再現できなかった。 原因としては、やはり区画数が少ない？のかもしれない。
希釈係数 $d = 22$ の結果
パラメータ設定は $d$ 以外は前と同じ。 元論文では $d=22$ において減衰振動が見られたので、それを再現することが目的である。
減衰振動ではあるが、元論文のように綺麗に減衰しなかった。
議論
元論文で課題として挙げられていることを見つけた範囲で書いておく。
- $\lambda \rightarrow \infty$ の場合に相を分離する簡単な式は存在するか？
- より多数の相互作用がある系を考慮するとどうなるか
- 統計物理ではよくある定式化であり、quasi-species や composomes の研究が参考になると書かれている。
- 非平衡統計的な扱い ([Blokhuis et. al., 2018], [Fujii et. al., 2013] など)
- 区画化のダイナミクスを拡張する（[Kamimura et al., 2019], [Yoshiyama et al., 2018] など）