Skip to content

Latest commit

 

History

History
367 lines (237 loc) · 21.1 KB

bio_turing_pattern.md

File metadata and controls

367 lines (237 loc) · 21.1 KB

Title: チューリングパターンとは Date: 2019-04-04 Category: 生物学 Tags: 生物学, Bioinfo Slug: turing_pattern Authors: Naoto Yamaguchi Summary: チューリングパターンの解説

動機

システム生物学の授業が書くことになったきっかけです。 hatenablog に書いたところ数式を上手く表示できず書き途中であったため、こちらで書きました。hatenablog の方に、意外とアクセスがあったので、折角訪問してくださった方にちゃんと見ていただきたいと思いました。hatenablog の方は放置していてひどいです。こちらのサイトも見た目等変なところが多いですがご容赦下さい。そのうち直します。

参考

実装は、以下のサイトを参考にしました https://ipython-books.github.io/124-simulating-a-partial-differential-equation-reaction-diffusion-systems-and-turing-patterns/

チューリングパターンとは？

誤解を恐れずに言えば、こういう模様(=波)のことです。

こういう模様、パターンのことをチューリングパターンと言います。チューリング波、反応拡散波という言い方もあります。

そして重要なのは、以下の 2 点

  1. このパターンはある条件を満たす化学反応システムが自発的に生み出す周期的なものであるということ
  2. 模様(=波)の仕組みは、反応拡散方程式で数学的に表せる！！ ということ
  3. この模様は、生物によく見られるものであり、一部の生物では、確かに反応拡散系によってこの模様が作られている！！ ということ

遺伝子だけで決まる訳ではない

point①&③

これは化学反応によって自発的に起こるパターンで、生物にも見られる！！ ※ 画像はイメージです。
こちらからお借りしました。チューリング波（反応拡散波）を理解したい

実際の生き物に見られる模様のパターン

ことの始まりは、イギリスの数学者アランチューリングが 1952 年に発表した論文で、「生物の模様は波によって作られる」という仮説。です

生物の体の表面では化学反応が起きており、反応を活性化する因子と抑制する因子の広がる速さが異なることで「波」が生じ、模様が作られる

https://www.jst.go.jp/pr/jst-news/pdf/2015/2015_04_p08.pdf という考え方。

1952 年に、イギリスの代表的な数学者でコンピュータ科学の生みの親でもあるアラン・チューリングが、「2 つの仮想的な化学物質が、ある条件を満たして互いの合成をコントロールしあうとき、その物質の濃度分布は均一にならず、濃い部分と薄い部分が、空間に繰り返しパターン（反応拡散波）を作って安定する」ことを、数学的に証明した。1970 年代に数人の数学者がチューリングの方程式を 2 次元でシミュレーションしたところ（チューリングの時代はコンピュータがなかった）、方程式の定数（仮想的な化学物質の性質）を少し変えるだけで、シマウマのストライプ模様もキリンの網目模様も、豹の斑点模様も作り出せることを発見した。

https://www.brh.co.jp/seimeishi/journal/011/to_1.html 天才か...

簡単に例えてみると... 2 つの赤い物質(A)と青い物質(B)があって、この 2 つの関係が、

  1. 赤い物質が増えると、自分(赤い物質)をより増やすようになる（=赤が自分を「促進」）

  2. 赤い物質が増えると、青い物質を増やすようになる(=赤が青を「促進」)

  3. 逆に青い物質が増えると、赤い物質を減らすようになる(=青が赤を「抑制」)

    ような関係で、

  4. 赤い物質より青い物質の方が周りに広がりやすい

    とき、ランダムに赤青がある状態からスタートして、ある程度時間が経つと、安定した赤青模様のパターンができる！ということである。

    上記のように初期状態やパラメーターによって色々な模様ができます。

白黒でも良いですが、画像に合わせて赤青にしました。

この関係は、食物連鎖の一部に例えてみるとわかりやすいでしょう全然現実世界に沿っていないとは思いますが...

  1. ネズミが増えるとネズミたちは、協力して生活するようになり、どんどんその数を増やしていきます(ネズミが自分たちの増加を「促進」)
  2. ネズミが増えると、ネズミを食べる猫たちは沢山食料を得られるので、猫たちの数も増えていきます(ネズミが猫の増加を「促進」) 3.　猫が増えると、ネズミは今まで以上に食べられてしまうので、数を減らしてしまいます。(猫の増加がネズミの増加を「抑制」)

このような関係のイメージです。

ポイント 2 模様の仕組みは数学的に表せる！

この 1~4 の関係を数学的に表すと、一例として、

$$ \frac{\partial u}{\partial t} = 0.6p - q - p^3 + 0.0002\Delta u$$

$$ \frac{\partial v}{\partial t} = 1.5p - 2q + 0.01\Delta u$$

と表すことができます。これは反応拡散方程式と言われるものです。式の意味は後ほど説明します。

反応拡散方程式とは？？

位置付けとしては、今回の式は、2 変数の反応拡散方程式

反応拡散方程式 : $$ ∂u/∂t = F(u) + D∂^2u/∂x^2 $$

反応方程式 : $$ ∂u/∂t = F(u)$$

拡散方程式 : $$ ∂u/∂t = D∂^2u/∂x^2 (D は拡散係数)$$

1 つめは反応項、2 つめを拡散項と呼ぶことにします。 F(u)は 2 物質の化学反応による物質の濃度の変化を表す式が入ります。例えば、酵素の反応速度を表すミカエリス・メンテンの式などが入ることもあるでしょう。

拡散項は、2 物質が空間を広がっていくことを考えるときに使います。

2 つの化学物質 u と v があって、その二つの濃度が以下のような数式（反応拡散方程式）で表わされるとき、u と v の濃い部分と薄い部分がいい感じに模様のパターンを作って安定化するということである。

キリンやシマウマは違うっぽい

バーゼル市の動物園へ行って、シマウマやキリンの模様を詳しく観察してみた。すると、それらは反応拡散波そのものではないことに気がついた。反応拡散波は、化学反応が作り出す定常波であるから、その波長（模様の間隔）は、一定に保たれなければならない。しかし、シマウマやキリンの模様は皮膚に固定されているらしく、体が大きくなるにつれて模様の間隔が広がっていってしまう。

近藤先生が、実験で証明したのは、タテジマについてです。シマウマやキリンについてはわかりません。近藤先生の物語はこちら https://www.terumozaidan.or.jp/labo/technology/15/index.html 生物が模様を作る仕組みを 25 年かけて明らかにした。それを我々は学ぶことができている。感謝

反応拡散系(=反応系 + 拡散系) 化学反応と分子の拡散を組み合わせた反応システムのこと。

http://www.fbs.osaka-u.ac.jp/labs/skondo/ozaki/what%20is%20RD%202(outline).htm

チューリングマシンとは？（ついでに）

todo 情報科学の勉強しましょう。形式言語理論を学ぶとわかります。

というかアラン・チューリングって誰だよ？

めっちゃすごい人。計算機科学の父。映画、「イミテーションゲーム」はおすすめです。

セルオートマトン、ライフゲームって何？

チューリングは原論文（The Chemical Basis of Morphogenesis）の中で、反応拡散波が皮膚の模様形式だけでなく、「形態形成全般にはたらく基本的なメカニズム」であろうと述べている。

つまり骨格なども形成しているのではという考え方もあるようです。

反応拡散方程式の導出(2 変数、1 次元)

三浦岳、発生の数理を参考

考える系について

まずは一次元の棒状の組織を考えてから二次元に拡張する。横の長さ 1, 縦の長さ$$ dy $$として、$$ dy << 1 $$で縦方向の分布は無視する。この棒状の組織の中の細胞が 2 種類の分子、活性因子(activator)と抑制因子(inhibitor)を産生し、これらの分子が細胞の分子の産生、分解をコントロールしつつ、近傍の細胞に拡散している、という状態を考える。

離散化

このままでは空間的な分布の状態を考えるのは難しいので、この細長い棒を横方向に$$ dx $$ずつ区切って、$$ 1/ dx $$個の小さな組織片に分けて考える。($$ dx << 1 $$) 各組織片内での活性因子、抑制因子の分布はほぼ一様であるとして、

  • 組織片の中での活性因子、抑制因子の相互作用
  • ある組織片と、両隣りの組織片との間での活性因子、抑制因子のやりとり

の 2 つの要素を考え、時間変化をみる。時間は連続的なものだが、便宜上、微小区間$$ dt $$に区切って考える。時刻$$ m * dt $$($$m$$は整数)に置ける、それぞれの組織片の中の活性因子、抑制因子の濃度について、左から数えて n 個目の組織片の中の活性化因子の濃度を$$ p(n,m) $$、抑制因子の濃度を$$ q(n, m) $$とする。

反応項

反応項は 2 種類の分子の相互作用によって決まり、様々なパターンがある。今回は一例として、チューリングパターン形成に置ける反応項を考える。今回考えている、チューリングパターンを生成する反応拡散系における反応系では、上で述べたように、

  • 活性化因子は、自身の産生を促進
  • 活性化因子は、抑制因子の産生を促進
  • 抑制因子は、活性因子の産生を抑制

という相互作用をしていると考える。

これを数式で表すと、 $$ f(p, q) $$, $$ g(p, q) $$ を活性因子、抑制因子の濃度の変化率として、

$$ f(p, q) = ap - bq (a, b は正の定数)$$

$$　 g(p, q) = cp -dq (c, d は正の定数)$$

変化量を考えると $$ f(p(n, m), q(n, m)) \times dt $$

$$ g(p(n, m), q(n, m)) \times dt $$

となる。

今回の例では、2 種類の物質で考えたが、他にも 1 種類、3 種類での反応系を考えることができて、それらの相互作用(促進か、抑制か)に応じて、様々なダイナミクスが生じる。具体例としては、

  • 活性因子は、ある程度以上の濃度になると産生が抑えられる(-p3) ∂u/∂t の反応項 = ∂v/∂t の反応項

拡散項

次に、隣の組織片との間の相互作用を考え、拡散における濃度変化を考える。 $$n$$番目の組織片の活性化因子の濃度を$$p(n, m)$$とすると、右隣の組織片の活性因子の濃度は、$$p(n+1, m)$$となる。単純な拡散の場合、物質の移動量は濃度勾配と境界の長さに比例する(Fick の法則)ので、微小時間 dt の間に右隣の組織へと移動していく活性化因子の量は、

$$ d_p(p(n+1, m) - p(n, m)) / dx \times dy \times dt $$

で表される。$$d_p$$は活性化因子の拡散係数である。これによる$$n$$番目の組織片の濃度変化は、この量を組織片の面積で割って、

$$ d_p(p(n+1, m) - p(n, m)) / dx^2 \times dt $$

となる。同様に左隣の物質とのやりとりを考えると、左右合わせて、$$n$$番目の組織片の濃度変化は、

$$ d_p(p(n+1, m) + p(n, m) - 2 \times p(n, m)) / dx^2 \times dt $$

となる。

同様に抑制因子の濃度変化は、

$$ d_q(q(n+1, m) + q(n, m) - 2 \times q(n, m)) / dx^2 \times dt (d_q は抑制因子の拡散係数)$$

これが拡散項となる

反応拡散方程式

領域が十分広いという境界条件のもとで、$$dt, dx$$を無限に 0 に近づけると、離散ではない連続の方程式を得ることができる。連続関数として、時刻$$t$$における位置$$x$$の活性因子、抑制因子の濃度をそれぞれ$$u(x, t), v(x, t)$$とし、$$dt, dx$$を無限に 0 に近づけると、

$$ ∂u(x, t) / ∂t = f(u(x, t), v(x, t)) + d_p∂^2u(x, t) / ∂x^2$$

$$ ∂v(x, t) / ∂t = g(u(x, t), v(x, t)) + d_q∂^2u(x, t) / ∂x^2$$

とそれぞれの濃度変化を表すことができる。それぞれ第一項が反応項であり、第二項が拡散項である。これは一般の反応拡散方程式の表記

2 次元チューリングパターンの実装

ここでは、反応項として、

$$f(p, q) = 0.6p - q - p^3 $$

$$ g(p, q) = 1.5p - 2q $$

を考える。(ちなみにこれは神経細胞の発火モデルである、FitzHugh-南雲モデルと類似している)
これは、

  1. 活性因子が、自らの産生を促進する$$(0.6p)$$
  2. 活性因子は、抑制因子の産生を促進する$$(1.5p)$$
  3. 活性因子は、ある程度以上の濃度になると産生が抑えられる$$(-p^3)$$
  4. 抑制因子は、活性因子の産生を抑制する$$(-q)$$
  5. 抑制因子は、何もないと減衰していく$$(-2q)$$

という 5 つの条件を備えている。

拡散項については、 $$ d_p = 0.0002 , d_q = 0.01$$ とおく。
$$ d_q $$の方が$$ d_p $$よりも大きい値でないと、安定なパターンが形成されないことがわかっている。意味合いとしては、物質$$v$$の方がより拡散しやすいということである。

濃度変化は以下のように表せる。

$$ \frac{\partial u}{\partial t} = 0.6p - q - p^3 + 0.0002\Delta u$$

$$ \frac{\partial v}{\partial t} = 1.5p - 2q + 0.01\Delta u$$

これを元にして、2 次元でのパターン形成の様子を描画してみる。

$$ \frac{\partial u}{\partial t} = f(u, v) + d_p\frac{\partial^2 u(x, t)}{\partial^2 x}$$

$$ \frac{\partial u}{\partial t} = f(u, v) + d_p\Delta u$$

$$ \frac{\partial v}{\partial t} = g(u, v) + d_q\frac{\partial^2 v(x, t)}{\partial^2 x}$$

$$ ∂v(x, t) / ∂t = g(u(x, t), v(x, t)) + d_q∂^2u(x, t) / ∂x^2$$

実装

1. パッケージのインストール

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

2. チューリングパターンのシミュレーションに、以下の偏微分方程式を用いる。範囲は-1~1

$$ ∂u / ∂t = aΔu + u - u^3 - v + k $$ $$ τ ∂v/∂t = bΔv + u -v $$

u は色素？濃度を表す。v もその濃度。v は u は抑制する。u は v を促進する。初期条件では、u,v は互いに独立したランダムな値とする。 Neumann 境界条件

3. 以下のパラメータを定める

a = 2.0e-4 b = 1e-2 tau = 1 k = 0

4. 時間と空間を離散化して範囲に分ける。

size = 100 dx = 2. / size T = 54.0 dt = .001 n = int(T / dt

5. u と v を初期化。

U = np.random.rand(size, size)
V = np.random.rand(size, size)

def laplacian(Z):
    Ztop = Z[0:-2, 1:-1]
    Zleft = Z[1:-1, 0:-2]
    Zbottom = Z[2:, 1:-1]
    Zright = Z[1:-1, 2:]
    Zcenter = Z[1:-1, 1:-1]
    return (Ztop + Zleft + Zbottom + Zright - 4 * Zcenter) / dx ** 2

7.行列を表示する関数を定義

def show_patterns(U, ax=None):
    ax.imshow(U, cmap=plt.cm.copper, interpolation='bilinear', extent=[-1, 1, -1, 1])
    ax.set_axis_off()

    U = np.random.rand(size, size)
    V = np.random.rand(size, size)

    def laplacian(Z):
        Ztop = Z[0:-2, 1:-1]
        Zleft = Z[1:-1, 0:-2]
        Zbottom = Z[2:, 1:-1]
        Zright = Z[1:-1, 2:]
        Zcenter = Z[1:-1, 1:-1]
        return (Ztop + Zleft + Zbottom + Zright - 4 * Zcenter) / dx ** 2

    fig, axes = plt.subplots(3, 3, figsize=(8, 8))
    step_plot = n // 9
    # 偏微分方程式を有限微小変化methodでシミュレーション

    for i in range(n):
        # uとvのラプラシアンを計算
        deltaU = laplacian(U)
        deltaV = laplacian(V)
        # uとvの値をグリッドの中でとる
        Uc = U[1:-1, 1:-1]
        Vc = V[1:-1, 1:-1]
        # 値を更新する
        U[1:-1, 1:-1], V[1:-1, 1:-1] = \
            Uc + dt * (a * deltaU + 0.6 * Uc - Uc**3 - Vc + k), \
            Vc + dt * (b * deltaV + 1.5 * Uc - 2.0 * Vc) / tau
        # Neumann条件: 端の微分は0?
        for Z in (U, V):
            Z[0, :] = Z[1, :]
            Z[-1, :] = Z[-2, :]
            Z[:, 0] = Z[:, 1]
            Z[:, -1] = Z[:, -2]

        # 9つの異なる時間に状況を分けてプロット
        if i % step_plot == 0 and i < 9 * step_plot:
            ax = axes.flat[i // step_plot]
            show_patterns(U, ax=ax)
            ax.set_title(f'$t={i*dt:.2f}$')
fig , ax = plt.subplots(1, 1, figsize=(8, 8))
show_patterns(U, ax=ax)

感想

今回のお話は、システム生物学という分野で扱われがちなお話です。私が興味を持っているバイオインフォマティクスにシステム生物学が含まれるかどうかはわかりません。僕が所属する学科は、バイオインフォマティクスとシステム生物学という 2 つの学問分野を柱にしています。システム生物学もなかなか面白いなと思いました。今後こちらのブログも充実させていくので、是非チェックして見て下さい！ todo あとで追記