108

拡散異方性によるチューリングパターン

Turing Patterns by Anisotropic Diffusions
岩本凌, 昌子

浩登

関西学院大学理工学部

Ryo Iwamoto, Hiroto Shoji
School of Science and Technology, Kwansei Gakuin Univ.

1

はじめに
自然界では、様々なスケールで一様な状態から自発的に形成される現象が見られる。これまで、化学反応で形成されるパターン[1], 生物の形態形成 [2] などでは様々な研究が行われ、その数理的メカニズムの解析研究が行われてきた。

その中でも、化学反応では Chlorite‐Iodide‐Malonic Acid (CIMA) 反応 [3] 、生物では魚の体表模様 [2, 4] や指の形態形成 [5] は、1952年アラン・チューリング [6] が提唱した2種以上の構成要素の相互作用を表す微分方程式系 (反応項) とそれぞれの構成要素の空間拡散の項がカップルした
いわゆる反応拡散方程式において自発的に形成されるパターンとの対応が議論されてきた。
これまで、このチューリングの提唱した機構については数理解析的においても、数値計算をメインにした研究においても、たくさんの研究が行われてきた [1, 2] 。2次元システムにおいて、この機構の数値計算を行うと、式中のパラメータによって縞模様や斑点模様が生成されることが知られている [2] 。また、縞模様が生成されるようにパラメータを設定し、初期分布を変えて数値計算を行うと、毎回異なる方向の縞が形成されることが知られている。しかし、この機構のコンセプトを導入してモノづくりの仕組みを考えるためには、形成される縞の方向は初期分布によって変わらない仕組みを考えたい。

著者らはこれまで、方向によって流れの速度が加速もしくは減速するという異方性を拡散項に導入したモデルについて報告を行ってきた [7] 。本研究では、よりシンプルな拡散異方性である方向によって拡散レートが異なる異方的な拡散を導入したモデルを作成し、モデル解析を行った。形成される縞の方向性を定量化する指標を設定し、形成されるパターンの数理的解析を行った。

2チューリングパターン

1952年、イギリスの数学者アランチューリングは、卵のような一様な状態から生体のような複雑な構造が自発的に形成するのはなぜかを考えた。そして、周期構造形成の数学的メカニズムの1つとして、次のような数学的メカニズムを提案した。2変数の連立偏微分方程式\frac{}=D_{u}(\frac{\partil^{2}u\partilx^{2}\partil^{2}v\partilx^{2}+ \frac{\partil^{2}u\partilv,\partily^{2}\partily^{2})+f(u,v) \frac{\partilv\partilu\partil}{\partil}=D_{v}(\frac{}+\frac{})+ g(u,v)(1)

109を考える。それぞれ変数 u,についての変化式を表していて、右辺第1項は拡散項を表し、vD_{v} は拡散係数である。そして右辺第2項は u とvの多項式で構成されていて、

D_{u},f(u, v) と g(u, v)

は反応項とよばれる。このような拡散と反応をカップルした形の方程式を反応拡散方程式という。
チューリングは、これらの方程式において、拡散係数の比が十分大きく、反応項がある条件を満たせば、空間的に非一様な周期解が安定に存在することを示した [6] 。
具体的には、 \overline{u}_{\backslash }\overline{v} を、反応項 f(\overline{u}, \overline{v})=0 、 g(\overline{u}, \overline{v})=0 を満たす時間空間変化しない解がただ一つであるとする。 u=\overline{u}+ae^{\lambda t+iqx_{\ovalbox{\tt\small REJECT}}}v=\overline{v}+ be^{\lambda t+iqx} を方程式 (1) に代入し、平衡点 (\overline{u}, \overline{v}) において線形化すると、\lambda(\begin{ar ay}{l a b \end{ar ay})=(\begin{ar ay}{l} -q^{2}D_{u}+f_{u} f_{v} 9u -q^{2}D_{v}+g_{v} \end{ar ay})(\begin{ar ay}{l a b \end{ar ay})(2)が得られる。ここで、 f_{u}, f_{v}, g_{u}, g_{v} は f, g の u, v に関する (\overline{u}, \overline{v}) での偏微分係数である。波数 q は空間の変調を表すパラメータである。(2) が a , わが自明でない解を持つためには、\lambda^{2}+P(q)\lambda+Q(q)=0 ,(3)でなければならない。ここで、P(q)=q^{2}(D_{u}+D_{v})-f_{u}-g_{v},Q(q)=D_{u}D_{v}q^{4}-q^{2}(D_{u}g_{v}+D_{v}f_{u})+f_{u9v}-f_{v9u},を表す。方程式 (1) において、拡散項がない場合の時間変化しない解の線形安定性は、式(3) でq=0 を考え、\lambda^{2}-(f_{u}+9_{v})\lambda+f_{u}g_{v}-f_{v9u}=0 .

(4)

ここでは2次方程式 (3) の解が共に負、つまり f_{u}+9v<0 とん9。一ん仇 >0 を満たすとする。
(3) の固有値 \lambda の1つが正になる条件、すなわち不安定性の十分条件は Q(q) である。
Q(q) の最小値を求めるため、 Q(q) を q で微分すると最小値を得る q は、q_{c}^{2}= \frac{D_{u}g_{v}+D_{v}f_{u} {2D_{u}D_{v} ,であり、その時の Q の値は、Q(q_{c})=- \frac{(D_{u}g_{v}+D_{v}f_{u})^{2} {4D_{u}D_{v} +f_{u}g_{v}-f_{v} g_{u},である。よって、ある有限の q で Q(q) が負になるようにパラメータを設定すれば、一様解は有限の波数を持った変調に対して不安定になる。これを拡散誘導不安定性あるいはチューリング不安定性とよび、これらのメカニズムから形成されるパターンをチューリングパターンと呼んでいる[2] 。
反応項に、例えばFitzHugh‐ 南雲方程式 [2]\{ begin{ar ay}{l f(u,v)=u- ^{3}-v, g(u,v)=\gam a(u-\alphav-\beta), \end{ar ay}を用いると、 D_{u}=1.0\cross 10^{-4}, D_{v}=1.0\cross 10^{-2}, \alpha=0.50 , î
安定性が起こる \beta の範囲は、(5)=-0.099<\beta<0.099 と算出できる。

26.0とすると、チューリング不110
(a)

(b)

図1: 方程式 (1) と (5) の2次元数値計算で得られたチューリングパターン(濃度プロット) : (a) 縞模様 (\beta=0.00), (b) 斑点模様 (\beta=0.08)

この範囲のパラメータ \beta を与え、平衡点 (\overline{u}, \overline{v}) に少しのノイズを加えた初期分布を使い数値計算を行うと、1次元空間では u, v の同位相の1次元周期的なパターンが形成される。
同様に2次元空間において数値計算を行うと、パラメータ \beta によって図1のように縞模様や斑点模様が生成されることが知られている。

また、偏微分方程式の初期値問題として知られるように、初期分布に与えるノイズの仕様を変えて分布をそれぞれ数値計算すると、その都度異なる方向の縞が形成されることが知られている[2] 。実際に縞模様が得るようにパラメータ \beta=0.00 に固定して、初期分布を変えたものを与える。
すると、得られる縞パターンの方向性はその都度異なった。

3異方性を導入したチューリングモデル
(a)

(b)

\frac{\partil^{2}u\partilx^{2}+\frac{\partil^{2_L}{\partily^{A} \prime}
図2: 異方的拡散の概念図 (a) 通常の拡散 (b) 方向性のある拡散本研究では、図2(a) のように濃度の高いところが等方的に広がっていく通常の拡散ではなく、図2(b) のような拡散しやすい方向性があるという拡散異方性を考えた。
またこの異方性は u 、vそれぞれ違う方向性をとれるものとして、それぞれ u の異方性の強さ a 、v の異方性の強さ
a>1 , or b>1 だと、そそして の異方性の強さ
れぞれ x 軸方向により広がりやすくなる。それに対してそれぞれの異方性の強さ a<1 , or b<1vb と定義した。それぞれ u 、だと、それぞれ y 軸方向により広がるという異方的な拡散を考え、下のように導入したモデルを

111
111
解析した。

\{ frac{} =D_{v}(b+2-b)\frac{)\partil}{\partil}+9(u,v)\frac{\partilu} {\partilv,\partil \partil }=D_{u}(a\frc{\partil^{2}u{\frac{} ,\partil \partil^{2}vx^{\Gam a}\partilx^{2}+(2-a\frc{\partil^{2}u{\partily^{2}, y^{2} v) +f(u,v)(6)(0<a<2,0<b<2)

数値計算では、式(5) で用いた FitzHugh‐ 南雲方程式 [2] の反応項を主に用いて以後の解析結果を示す。

縞模様が形成されるように、 D_{u}=1.00\cross 10^{-4}, D_{v}=1.00\cross 10^{-2}, \alpha=0.50, \gamma=26.0 、そして\beta=0.00 とパラメータの値は固定して計算し、 u の異方性の強さ a と v の異方性の強さ b をそれぞれコントロールパラメータとして様々に変化させて得られるパターンを比較した。
数値計算スキームについては、オイラー法を用い、格子空間のサイズは 100\cross 100 の正方領域において数値計算を行った。

初期分布は、平衡点 (\overline{u}, \overline{v}) の値に各点に空間相関のない微小摂動を加えたものを与えた。
領域の境界条件は、広い領域での様子を観察するために周期境界条件を課した。空間分割
1.0\cross 1.00^{-2}, dt=1.00\cross 10^{-4} として、値の発散が起きないように設定した。

dx=

このチューリング機構から形成されるパターンについて観察すると、初期に線形作用によって特定のモードが成長して初期パターンが形成される。初期パターンには、ローカルに成長したモードのズレであるディフェクトなどが存在し、これらを解消すべく、引き続き、非線形項の作用によって長時間かけて特定の縞の方向性が揃ってくるのがわかる。しかし、その非線形性の作用については形成されるディフェクトなどがパターン毎に違ってきたりして縞の方向性が揃う時間は変わってくる。
本研究では、このような非線形項による作用を考えるのではなく、初期に形成されやすいモードにおけるパターン形成を考えた。そのため、数値計算についてはその初期モード形成時問を考慮して、基準の時間を設定して、その時に形成されているパターンを比較した。本研究では、式(6) と (5) モデルでは t=100 の時点で得られたパターンを比較した。2次元プロットについては、
uの濃度プロットを示す。

vの分布は u の分布と同位相、同周期であるので省略した。
(c)

図3:

uuの分布の濃度プロット、(a)a=b=0.8, (b)a=b=1.0, (c)a=b=1.2.の異方性の強さ a と v の異方性の強さ b ともに1.0にした場合、ともに等方的な拡散である。

この時得られたパターンは、これまでの通常の拡散と同様なので、図3(b) のように縞の方向は定まらパターンが得られた。

112
まず、の異方性の強さ a 、 v の異方性の強さ b を同じにしながら、強くしたり (a=b>1) 、弱(a=b<1) してみた。図3はその結果である。図 3(a) は a=b=0.8 の時、また、図3(b) は a=b=1.2 の時の結果をそれぞれ示す。これらだけでなく、異方性の強さ a と b を同じに保ちながら、定義域の範囲 (0.0<a(=b)<2.0) 様々な値をとって数値計算を行って見ても、図のuくしたりように縞の方向は定まらなかった。

図4:

a=1.0 で固定し b を変化させた

次に、uuの分布の濃度プロット、(a)b=0.8, (b)b=1.0, (c)b=1.2の異方性の強さ a を等方的に、つまり a=1.0 に固定したまま、vの異方性の強さ b を変化させていった。その結果が、図4である。もちろん、図 4(b) のように v の異方性の強さが b=1.0の場合には方向性の定まらない縞が形成された。しかし、は等方的拡散にしたまま (a=1.0) 、 vが垂直方向により広がりやすくなるという異方的拡
uの異方性の強さを (b<1) にした時、つまり散を導入した場合、図 4(a) のように、水平方向の縞が得られた。図4(a) は、vを示しているが、今度は、vb=0.8 の時の結果b<1 の範囲で値を変化させてもやはり水平方向の縞が得られた。の異方性の強さを (b>1) にした時、つまりvが水平方向により広がりやすくなるという異方的拡散を導入した場合、図 4(c) のように垂直方向の縞が得られた。図 4(c) は、時の結果を示しているが、この場合もb=1.2 のb>1 の範囲で値を変化させてもやはり垂直方向の縞が得られた。

これらを合わせて考えると、uの拡散を等方的、つまり異方性の強さ a=1.0 にして、vに異方性を導入した場合、拡散しやすい方向と垂直な方向に縞が形成されることが分かった。

(a)

図5:

-b=1.0 で固定し方、v(b)a(c)を変化させた u の分布の濃度プロット、の異方性の強さ b を等方的に、つまり ó=(a)a=0.8, (b)a=1.0, (c)a=1.21.0に固定したまま、uの異方性の強さ a を変113化させていった。その結果が、図5である。もちろん、図5(b) のように u の異方性の強さが a=1.0の場合には方向性の定まらない縞が形成された。しかし、の異方性の強さを (a<1) にした時、つまりuvは等方的拡散にしたまま (b=1.0) 、拡散を導入した場合、図5(a) のように、垂直方向の縞が得られた。図5(a) は、果を示しているが、uが垂直方向により広がりやすくなるという異方的
a=0.8 の時の結a<1 の範囲で値を変化させてもやはり垂直方向の縞が得られた。

今度は、 u の異方性の強さを (a>1)\ovalbox{\t \small REJECT} こした時、つまり u が水平方向により広がりやすくなるという異方的拡散を導入した場合、図 5(c) のように水平方向の縞が得られた。図5(c) は、 a=1.2の時の結果を示しているが、この場合もa>1 の範囲で値を変化させてもやはり水平方向の縞が得られた。

これらを合わせて考えると、vの拡散を等方的、つまり異方性の強さb=1.0 にして、uに異方性を導入した場合、拡散しやすい方向に縞が形成されることが分かった。

4

縞模様の方向解析
縞の方向を定量的に議論するために、どちらの方向に明確に縞模様が形成されたか、もしくは方向性の定まらない縞が得られたかを判定するために指標を作成した。以下のように離散フーリエ変換(DFT) を用いた指標を作成した。まず、 u(x , 勿の2次元 DFT の q_{x}, q_{y} に関する分布、 F(q_{x}, q_{y})を次のように定義する。

|F(q_{x}, q_{y})|=| \sum_{x=0}^{N-1}\sum_{y=0}^{N-1}u(x, y)\exp(\frac{-i2\pi(q_ {x}x+q_{y}y)}{N})|.
i は虚数単位、N は格子数とする。

数値計算によって得られた u (or v) の分布を2次元 DFT を行い、次のように規格化したものを|F'(q_{x}, q_{y})| とする。

|F'(q_{x}, q_{y})|= \frac{|F(q_{x},q_{y})|}{\sum_{q_{x}=0}^{N-1}\sum_{q_{y}=0}^ {N-1}|F(q_{x},q_{y})|}.
この |F'(q_{x}, q_{y})| の濃度プロットを次のように用いて、どちらの方向に縞が形成されたか、またその縞の方向は定まっているのかどうかを次の観点から判定した。

(i): ピーク点が明確に2点に分かれているかどうか。
(ii): (i) を満たす時、 |F'(q_{x}, q_{y})| の2つのピークが明確かどうか。
それぞれ、(i) は |F'(q_{x}, q_{y})| の分布を目で見て判断し、(ii) は(i) を満たす時 |F'(q_{x}, q_{y})| の分散の値を比較した。つまり、2点ピークを持つ分布の |F'(q_{x}, q_{y})| の分散の値の閾値を用意しておき、それより閾値より分散の値が大きければ、縞の方向が定まっていると判定し、それ以下であると方向は定まっていないと判定した。

この判定基準を使って、異方性の強さ a 、b を0.05ずつ変化させたとき、どの方向の縞が得られたかを示したのが図6である。それぞれ、図 4(a), (b), (c) で得られたパターンの DFT 分析結果を示している。図 4(a), (c) はそれぞれ2点に値のピークが見られるが、(b) は特定の2点が定まらず、中心から一定の距離の点に分布が集まっている。また、詳細は省くが、図6(a), (c) の分布の分散を計算すると、両方とも設定した閾値を超えていた。これらのことから、図 6(a), (c) はそれぞれのピークの法線ベクトル方向の縞が形成されていると判断し、図4(a) は水平方向の縞が形成

114され、図4 (c) では垂直方向の縞が形成されたが、図4 (b) では、縞は特定の方向には形成されていないと判断した。

0\emptyset
0\infty0\alpha00”n\infty00\prime.
00\prime.

00\prime f0 \ovalbox{\t \smal REJECT}
00t001om

om
0\propto

0\infty ’om

図6: それぞれ図4(a), (b), (c) で得られたパターンの DFT 分析結果。
この判断指標を用いて、異方性の強さ a, b を様々に変えて見た時の結果が、図7である。

\blacksquare :垂直方向\bullet :水平方向▲:定まらない

図7: 縞方向の相図

注意として、異方性の強さはの強さを0.3<a<1.7 、0.0<a<2.0 、0.0<b<2.0 で理論的には定義できるのだが、こ0.3<b<1.7 以外の範囲に設定してしまうと、縞の間隔が非常に狭いパターンが得られ、 u の濃度プロットが判別できないことや、上で設定した方向判定指標の効能が効かない。そのため、本研究では 0.3<a<1.7 、 0.3<b<1.7 の範囲で縞の方向性解析を行い、判定を行った。

115
図7を見ると、

a=b またはaと b にほとんど差が無い場合は、いくら異方性の強さを強くしても縞の方向が定まらないという結果が得られた。しかし、異方性の強さについて物質問で差があれば縞の方向性が定まることがわかった。つまり、a<b のときは垂直方向の縞に、a>b のときは水平方向の縞となった。

5

数理解析
ここまでは、数値解析によって得られる結果であった。ここからその結果を、どちらの方向に縞が形成されるかを次のように仮説をもとに数理的に考えていく。
法線ベクトル (k_{x}, k_{y}) 方向の縞模様の分布 u と v を

u=\overline{u}+Acos(k_{x}x+k_{y}y)e^{\lambda t},
v=\overline{v}+Bcos(k_{x}x+k_{y}y)e^{\lambda t},とする。ここで、得られる縞模様は最も不安定な方向の縞が得られると仮定して考える。そのために、これらを (6) に代入し、最も不安定な縞の方向分布を探索する。
平衡点 (\overline{u}, \overline{v}) の周りで線形化し、式(3) のように整理すると、\lambda^{2}+P(k_{x}, k_{y})\lambda+Q(k_{x}, k_{y})=0

P(k_{x}, k_{y})=D_{u}ak_{x}^{2}+D_{u}(2-a)k_{y}^{2}+D_{v}bk_{x}^{2}+D_{v}(2-b) k_{y}^{2}-f_{u}-g_{v}
Q(k_{x}, k_{y})=(D_{u}ak_{x}^{2}+D_{u}(2-a)k_{y}^{2})(D_{v}bk_{x}^{2}+D_{v}(2- b)k_{y}^{2})(D_{u}ak_{x}^{2}+D_{u}(2-a)k_{y}^{2})g_{v}-(D_{v}bk_{x}^{2}+D_{v}(2-b)k_{y}^{2} )f_{u}+f_{u}g_{v}-f_{v}g_{u}が得られる。さP(k_{x}, k_{y})はチューリング不安定性がおこる条件より必ず正なので、定数項 Q(k_{x}, k_{y})ar ow>の値が最も小さい値をとるモード、つまり最も不安定なモードがパターンを形成すると考え、のモードがどちらの方向の縞模様になるかを Q(k_{x}, k_{y}) の 3D プロットから推定した。

図8は、図4(それぞれ DFT 解析結果が図6) で得られたパターンに対して、縦軸 Q(k_{x}, k_{y}) の3D プロットである。

図8 (b) (図 4(b) ならびに図 6(b) ) のように異方性の強さ a 、b 共に1.0としたとき、(k_{x}, k_{y}) の値に関わらず、極小点が原点から一定に距離になっていて、最も不安定な方向というのが定まっていないことがわかる。

これに対して、図8 (a), (c)(それぞれ図4(a),(c) ならびに図6(a), (c) ) のように、a、b に差がある場合はピークとなる (k_{x}, k_{y}) の点が存在し、最も不安定な縞の方向というのが定まることがわかる。さらに、その方向が図 4(a), (c) ならびに図6(a), (c) で示されている方向と一致していることがわかる。
これらの一致から、拡散異方性の影響により不安定な縞の方向制御が行われ、縞の方向が定まると考えられる。

116
(a)

(b)

(c)

図8: 定数項 Q の (k_{x}, k_{y}) の依存性、(a) b=0.8, (b)b=1.0, (c)b=1.26

モデルの一般性
ここまで、FitzHugh‐ 南雲方程式を用いて議論してきた。このモデルだけでなく、我々は他の反応項モデルで、チューリング機構を示すものについても、これらの事実が成り立つことを検証してきている。

例えば、下の Schnakenberg モデル [8]\{ begin{ar ay}{l f(u,v)=\gam a(\alpha-u+u^{2}v), g(u,v)=\gam a(\beta-u^{2}v), \end{ar ay}(7)を反応項に用いて、チューリング不安定性が起こり、また縞模様が形成されるようなパラメータを
D_{u}=1.00, D 。 =20.0, \alpha=5.00\cross 10^{-2}, \beta=1.55, \gamma=1.00x10^{4} を固定し、異方性の強さ a, b を変化させると数値計算を行った。すると、 u の分布と v の分布が反位相になることは、FitzHugh‐
南雲方程式を用いたモデルと異なるが、得られる縞の方向性については、異方性の強さが u と vで同程度なら、縞の方向は揃わず、また u と v との間で異方性の強さに差があれば縞の方向は特定の方向に形成されるという同じような縞の方向性の結果が得られた。また、数理解析についても同様に得られる縞の方向性と一致することが確認できている。

7

まとめ拡散項にシンプルな形で異方性を導入するというモデルを作成し、形成される縞の方向について解析を行なってきた。このモデルでは、u、vの両方に同程度の異方性の強さを導入しても縞の方向は偏微分方程式本来の性質である初期値依存性があり、縞の方向は定まらない。しかし、u,v間で異方性の強さにある程度以上差があれば、形成される縞の方向は決まってくる。
この数理的メカニズムについて、最も不安定なモードがパターンを形成するという仮定のもと検証すると、形成される縞のパターンの方向性とその不安定モードの方向性が一致することから、異方性を導入することでその方向制御が行われるモデルであると考えられる。

117
参考文献
[1] ‐.Wave and Patterns in Biological and Chemical Excitable Media, edited by B. Krisky andH. Swinney, (1991) North‐Holland, Amsterdam.

[2] Murray, J. D. (2003) Mathematical Biology Springer.
[3] V. Castets, et al (1990) Phys. Rev. Lett, 64,2953[4] S. Kondo and R. Asai. (1995) , Nature 376,765.
[5] Rushikesh Sheth, et al Hox Genes Regulate Digit Petterning by Controlling the Wavelengthof a Turing‐Type Mechanism Science 2012: 338 no. 6113 pp. 1476‐1480[6] Turing, A. M. (1952) Phil. Trans. Roy. Soc. 237: 37‐72.
[7] Shoji, H. et al (2003a) Dev. Dyn. 226: 627‐633.
[8] J. Schnackenberg. Simple chemical reaction systems with limit cycle behaviour. J. Theor.
Biol., 81:389‐400, 1979.