第3回 2019/6/28

• 講師：黒田 真也
• 教科書
  • システムバイオロジー 岩波書店（現代生物化学入門８）
  • システム生物学入門 -生物回路の設計原理-
• 参考書
  • ファイマン物理学<1>力学
  • やさしい機械制御
  • 力学系入門
  • 非線形制御
  • 物理法則はいかにして発見されたか

簡単な生化学反応の特性

※前回同様、なぜそうなるのか？の部分を、『岩波教科書 現代生物科学入門第８巻「システムバイオロジー」第４章 生命現象の動的特性 付録A～E pdfファイル』を参考に数式を通して学んだ。

なお、以下を意識すべきである。

C-7)単純な生化学反応における周波数応答の解析

単純な生化学反応を考える。 1. \(A\overset{k_f}{\underset{k_b}{\rightleftarrows}}B\)

このとき、微分方程式は、以下で表される。

したがって、一般解を求めると

となるので、初期条件として \(t=0\) のとき \(A=1\) とすると、

これは、時定数が \(\tau = \frac{1}{k_{f}+k_{b}}\) となっただけで、前回暑かった1次反応のシステムと同じ挙動を示すことがわかる。

1. \(F(t) \longrightarrow A\overset{k_f}{\underset{k_b}{\rightleftarrows}}B\)

続いて、周波数応答解析を行う。以下のように、システムに正弦波刺激 \(F(t)\) を加える。

すると、微分方程式は、以下のようにかける。

ここで、\({(\mathrm{C.} 17)} + {(\mathrm{C.} 18)}\) より、\(\dot{A}=F-\dot{B}\) であり、\({(\mathrm{C.} 18)}\) の両辺を \(t\) で微分すると、\(\ddot{B}=k_{f} \dot{A}-k_{b} \dot{B}\) となるので、これらを用いて、

ここで、\(F=\hat{F} e^{i \omega t}, B=\hat{B} e^{i \omega t}\) と複素数表示して式 \((\mathrm{C.}21)\) に代入すると、

よって、

なので、gainを求めると、

これから、

• \(\omega^2\gg\left(k_f + k_b\right)^2\) で \(\mathrm{gain} = \frac{k_f}{\omega^2}\)
• \(\omega^2\ll\left(k_f + k_b\right)^2\) で \(\mathrm{gain} = \frac{k_f}{\omega\left(k_f + k_b\right)}\)

となるので、これらの好転を求めて、

以上で、カットオフ周波数が求まった。

D 安定性と振動

D-1)バネの運動方程式とのアナロジー

以下の図のように、質量 \(m\) の物体が、ばね定数 \(k\) のばね、および粘性摩擦係数 \(c\) を持つダッシュポットで支えられている系を考える。（ダッシュポットは物体の速度 \(\dot{x}\) に比例して \(c\dot{x}\) という抵抗力（摩擦力）を生じる。）

今、物体へ外力 \(F(t)\) が作用するとき、この系の運動方程式は、以下で与えられる。

なお、それぞれの力は以下の通りである。

ここで、式 \((\mathrm{D.}1)\) の両辺を \(m\) で割って、\(\gamma = c/m, \omega_0^2 = k/m\) とおくと、（\(\omega_0\) を固有角振動数という。また、\(\gamma\) は見かけの摩擦係数とみなせる。）

このとき、外力 \(F\) が \(F(t) = \dot{F}\sin\omega t\) で与えられるとき（強制振動）を考える。

（入出力関係が線形のシステムでは）正弦波の入力を与えた場合、出力も入力と同じ周波数を持った正弦波となるはずなので、ここでも同様にして \(F = \hat{F}e^{i\omega t},x = \hat{x}e^{i\omega t}\) と複素数で解を仮定する。

簡単のため、\(m=1\) とすると、

すると、gainは、

より、

これを横軸に \(\omega\)、縦軸にgainをとってグラフを書くと、ある角周波数 \(\omega_r\) において一つの極値を取るグラフが描ける。

# 以下のプログラムで描ける。
import numpy as np
import matplotlib.pyplot as plt
X = np.linspace(0,1,1000)
w0 = 1; gamma = 1
gain = lambda x:1/np.sqrt((w0**2-x**2)**2 + gamma**2*x**2)
Y = np.vectorize(gain)(X)
plt.plot(X,Y)
plt.xlabel("$\omega$")
plt.ylabel("gain")
plt.show()

そこで、\(\omega_r\) を求めるため、gainを \(\omega\) の関数とみて、\(\omega\) について微分すると、

よって、\(\omega\neq0\) のとき、

となるので、gainは \(\omega_r=\sqrt{\omega_0^2-\frac{\gamma^2}{2}}\) のときに極大値

をとる。

特に、摩擦係数 \(\gamma\) が非常に小さい場合（\(\gamma\ll\omega_0\)）には、gainは \(\omega=\sqrt{\omega_{0}^{2}-\frac{\gamma^{2}}{2}} \cong \omega_{0}\) で極大値

を取ることがわかる。このように、外部から振動が与えられる場合に、「その外力の角振動数が系の固有角振動数に近づくとgainが急激に大きくなること」を共振現象といい、gainが最大となるときの角振動数 \(\omega_r\) を共振角振動数という。

D-2)概念的な生化学反応と力学系

\(x,y\) を分子、\(a,b,c,d\) を速度定数とし、以下の連立微分方程式によって表される反応を考える。

この式は、以下のように行列を用いて表せる。

ここで、

とおいて、行列 \(A\) の固有値を \(\lambda_1,\lambda_2\) とすると、\(\lambda_1,\lambda_2\) は固有方程式の解である。

ここで、\(\lambda_1\neq\lambda_2\) のとき、それぞれの固有値に対する固有ベクトルを

とすると、\(\mathbf{v}_{1},\mathbf{v}_{2}\) は一次独立だから、\( (\mathrm{D.}22)\) の一般解は

と表すことができる。このとき、固有値 \(\lambda_1,\lambda_2\) の取り得る値（正負、実虚）によって式 \((\mathrm{D.}26)\) の挙動が大きく変わる。それをまとめると、以下の表になる。

なお、以下のプログラムによって実際にプロットすることも可能である。

import numpy as np
import numpy.linalg as LA
import matplotlib.pyplot as plt

def plot_concentration(a,b,c,d, xmin=0,xmax=10,c1=1,c2=2):
    t = np.linspace(xmin,xmax,1000)
    A = np.array([
        [a,b],
        [c,d]
    ])

    eig_value, eig_vector = LA.eig(A)
    lambda1, lambda2 = eig_value
    v1,v2 = eig_vector
    alpha1,_ = v1; alpha2,_ = v2

    func = lambda t: c1 * np.e ** (lambda1*t) * alpha1 + c2 * np.e ** (lambda2*t) * alpha2
    x = np.vectorize(func)(t)

    plt.plot(t,x, color="#6da3fc")
    plt.title("$\lambda_1={:.2f},\ \lambda_2={:.2f},\ a={:.2f}$".format(lambda1,lambda2,a))
    plt.xlabel("t")
    plt.ylabel("oncentration of a reactant x")
    plt.hlines(0,xmin,xmax, color="#80160e")
    plt.show()

# eigen values are Real numbers.
plot_concentration( 1, 0, 1, 2, xmin=0,  xmax=3)  # ① lambda_1 > 0, lambda_2 > 0
plot_concentration(-1, 0,-1,-2, xmin=0,  xmax=3)  # ② lambda_1 < 0, lambda_2 < 0
# eigen values are Complex numbers.
plot_concentration(0, -1, 1, 0, xmin=0,  xmax=20) # ③ a = 0, Re[lambda] = 0
plot_concentration(2,  1,-5,-1, xmin=10, xmax=20) # ④ a > 0, Re[lambda] > 0
plot_concentration(-2,-1, 5, 1, xmin=10, xmax=20) # ④ a < 0, Re[lambda] < 0

D-4)ヌルクライン

ヌルクラインは、微分方程式系を解析するのに有効な方法の一つである。以下の微分方程式

の場合、ヌルクラインとは、\(\dot{x}_{j}=0\) つまり、\(f_{j}\left(x_{1}, \cdots, x_{n}\right) = 0\) によって定まる点の集合（つまり、平衡曲線）のことである。

まずは、簡単のため、以下のような生化学反応モデル（総和保存あり）のヌルクラインを考える。

\(a>0,b>0,x+y=C(\mathrm{const.})\) とすると、\(x>0,y>0,C>0\)

の直線がヌルクラインである。

D-5)ヌルクラインと生化学反応

続いて、一般的な経路を想定した（総和保存がない）生化学反応を考える。

なお、このときヌルクラインの平衡解が安定化不安定かどうか、また発散振動、減衰振動、収束するのか等に関してはヌルクラインのみでは判定できず、

の \(\det J,\mathrm{tr}J\) の値、及び \(J\) の固有値 \(\lambda_1,\lambda_2\) の関係によって決まる。

D-6)平衡点の分類

これらも、以下のプログラムによってプロットができる。

import numpy as np
import numpy.linalg as LA
import matplotlib.pyplot as plt

def plot_nurcline(a,b,c,d, xsize=10, ysize=10):
    J = np.array([
        [a,b],
        [c,d]
    ])

    det = np.linalg.det(J)
    tr = a + d
    eig_value, eig_vector = LA.eig(J)
    lambda1, lambda2 = eig_value

    xmin = -xsize; xmax = xsize
    ymin = -ysize; ymax = ysize
    X,Y = np.mgrid[xmin:xmax+1, ymin:ymax+1]
    dX = a*X + b*Y
    dY = c*X + d*Y

    plt.figure(figsize=(6,6))

    plt.quiver(X, Y, dX, dY)

    x = X[:,0]
    plt.plot(x,x*lambda1, label="$y={:.2f}x$".format(lambda1), color="#80160e")
    plt.plot(x,x*lambda2, label="$y={:.2f}x$".format(lambda2), color="#6da3fc")

    plt.title("$\det J = {:.2f},\ ".format(det) +\
              "\mathrm{tr}" + "J = {:.2f},\ ".format(tr) +\
              "1/4(\mathrm{tr}J)^2 = " + "{:.2f}$".format(1/4*(tr)**2))
    plt.xlim(xmin,xmax)
    plt.ylim(ymin,ymax)
    plt.grid()
    plt.legend()
    plt.show()

plot_nurcline( 1, 0, 1, 2) # ① detJ>0, trJ<0, detJ<1/4(trJ)^2
plot_nurcline(-1, 0,-1,-2) # ② detJ>0, trJ>0, detJ<1/4(trJ)^2
plot_nurcline( 1, 2, 3, 4) # ③ detJ<0
plot_nurcline( 1,-2, 3,-2) # ④ detJ>0, trJ<0, detJ>1/4(trJ)^2
plot_nurcline(-1, 2,-3, 2) # ⑤ detJ>0, trJ>0, detJ>1/4(trJ)^2
plot_nurcline( 1,-2, 3,-1) # ⑥ detJ>0, trJ=0, detJ>1/4(trJ)^2