目的のない勉強会

しがないニートです

【細胞の理論生物学】振動現象の数理モデル_その1_【ブラッセレータ】

www.youtube.com

しょうもない文章と導入

 昔、振動という現象に取り憑かれたことがあった。僕は、生き物で、生き物は止まっているように見えるけど、実は流れの中で生きている。そういうことを真摯に考え始めたときだったと思う。僕はまだ二十歳にもなっていなくて、なんとなく人に会うのが嫌で、大学の図書館で本を読んで、でも夕方になると唯一の友人(教室では僕より寡黙だったけど、僕の前では口数が多くなるような人間)と大学の奥まった校舎の誰も来ないソファで、意味のあるようなないような話をして、古い地方病院にあるような茶色の革張りの長椅子に腰掛けて、オレンジ色から徐々に青みがかりやがて光が失せていく、その日の終わりを告げる空を窓から時折眺めながら、自分もどこかに沈み込んでいくような、そんな日々だった。暗くなって文字が見えなくなるまでプリゴジンの本なんか読んでいた。

導入

   ベロウソフ・ジャボチンスキー反応は、旧ソ連の科学者であるベロウソフが発見し、ジャボチンスキーが確立した、ある化学反応である。細胞のミトコンドリア内で起こるクエン酸回路における酸化還元反応の実験をおこなっていたとき、容器の中の溶媒の色が振動して変化するのを彼らは発見した。クエン酸回路は数十の多種多様な触媒を含むやや複雑な反応であるので、ブリュッセルの研究者たちが重要なところだけを抽出しモデル化した*1。そうして世に現れたのがラッセレータ(Brusselator)という数理モデルだ。ブリュッセルとオシレータ(振動)を合わせてブラッセレータということらしい。

ラッセレータとその定式化

 ブラッセレータは二種類の成分 X, Y が、外部から流入してくる A, B という成分と反応する変哲もないモデルである。A,B の濃度は系の中で一定であるという仮定をおく。

定式化

 定式化(Formulation)とは、絵(現実やモデル)を数式として定める作業である。成分 X, Y の時間における変化を定式化するとこうなる。


\frac{dX}{dt} = A-(B+1)X+X^2Y \tag{1.13} \\
\frac{dY}{dt} = BX - X^2Y

 そして、A = 1, B = b として、反応速度定数 a, b を導入するとこうなる。


\frac{dX}{dt} = 1-(b+1)X+aX^2Y \tag{3.2} \\
\frac{dY}{dt} = bX - aX^2Y

ヌルクラインと固定点

 成分 X, Y は反応によって時間変化する。その変化はある時にはめちゃくちゃ大きいかもしれないし、小さいかもしれない、あるいは変化はないかもしれない。時間変化がない $\frac{dX}{dt}=0$の瞬間はいろんなバランスで訪れるような気がする。実際、数式上でもそうらしく、時間変化がない点をつなげて現れる線をヌルクライン(Nullcrine)と呼ぶ。

 $\frac{dY}{dt}=0$ のぶんのヌルクラインも描くことができるから、2本の線が現れるはずである。そして、その線が交われば、その点は X, Y どちらの時間変化もゼロである。そうした特殊な点を固定点(Fixed point)*2と呼ぶ。

ヌルクライン

 Xから始めよう。$\frac{dX}{dt}=0$とすれば式(3.2)はこうなる。Yについて整理する。


1-(b+1)X+aX^2Y = 0  \\
Y = \frac{(b+1)X-1}{aX^2} \tag{3.4}

 次にYについて、$\frac{dY}{dt}=0$とすれば式(3.3)はこうなる。Yについて整理する。


bX - aX^2Y = 0 \tag{3.5}\\
Y = \frac{b}{aX}

 テキストでは丁寧に$Y_f$とか$Y_g$と丁寧にそれぞれ置いている。また、$\frac{dY}{dt}=0$は$X=0$のときも成り立つので


X = 0 \tag{3.6}

もヌルクラインの一つである。数式を見てイメージしてみればXがなければ、Yもその材料がないのだから、作られるはずがない。ともあれ、合計3本のヌルクラインが得られた。

ヌルクラインの可視化

 パラメータ a=b=1 とした。式(3.3)を赤, (3,4)を青で可視化するとこうなる。

Rコード

library(tidyverse)

# Nullcrine
Yf = function(X,a,b){((b+1)*X - 1) / (a*X^2)} # (3.4)
Yg = function(X,a,b){b / (a*X)} # (3.5)

# Parameter
a = 1
b = 1

# Data
X_vec = seq(0,5, by=0.01)

df = tibble(
  X = X_vec,
  Yf = Yf(X_vec,a,b),
  Yg = Yg(X_vec,a,b)
)


library(ggdark)
# Visualization
ggplot(data=df)+
  geom_point(aes(x=X,y=Yf), size=0.5, col="red")+
  geom_point(aes(x=X,y=Yg), size=0.5, col="blue")+
  coord_cartesian(ylim = c(0, 5))+
  labs(title = "Nullcrines of the Brusselator (a=1, b=1)",
       x = "X",
       y = "Y")+
  dark_theme_gray()

固定点を求める

 固定点は2本の曲線の交点である。式(3.3)と (3,4)よりこうなる。


\begin{aligned}
\frac{(b+1)X-1}{aX^2} &= \frac{b}{aX} \\
(b+1)X-1 &= bX \\
bX + X -1 &= bX \\
X &= 1
\end{aligned}

 X=1 のとき、$Y = \frac{b}{aX}$より、$Y = \frac{b}{a}$となる。

 よって、固定点(X,Y)= (1, b/a) である。

   a=b=1 の時、(X,Y)= (1, 1)であるのでこうなる。

続く

 固定点を求めた後は、その周りの挙動を調べることになる。線形安定性解析のことである。

*1:どこをどう抽象化したのだろう?

*2:平衡点、不動点とも呼ぶと思う。なんでこんな複数の呼び方があるの?概念が導入された背景によると思うけど、一体。