この記事は現代化学(東京化学同人)の2022年4月号から2024年3月号まで連載された安藤耕司先生の「Pythonによる化学シミュレーション入門」を読んでやってみた記事です。詳しい内容は上記の「Pythonによる化学シミュレーション入門」をご参照ください。今回は、2022年10月号の振動反応の内容です。
目標
振動反応のシミュレーションをします。化学反応には、物質の濃度が時間的もしくは空間的に周期的な変化をするものがあります。これを振動反応と呼びます。振動反応はさまざまな研究がなされ、多くの機構やモデルが考えらていますが今回はそのうちの三つをシミュレーションしてみます。
Lotka-Volterra機構
次のような、Aが中間体X、Yを経てBになる多段階の化学反応を考えます
\displaylines{
A + X \longrightarrow 2X \\
X + Y \longrightarrow 2Y \\
Y \longrightarrow B \\
}
ここで、1段階目と、2段階目の反応は自己触媒反応です。
1段階目、2段階目、3段階目の反応速度定数を$k_1$、$k_2$、$k_3$とするとそれぞれの反応の反応速度は次のようになります。
\displaylines{
v_1 = k_1 A(t) X(t) \\
v_2 = k_2 X(t) Y(t) \\
v_3 = k_3 Y(t)
}
ここで$v_1$、$v_2$、$v_3$は1段階目、2段階目、3段階目の反応速度、$A(t)$、$X(t)$、$Y(t)$は時刻$t$でのAとXとYの濃度。$X(t)$、$Y(t)$は次のようになります。
\displaylines{
\frac{d X(t)}{dt} = v_1 - v_2 = (k_1 A(t) - k_2 Y(t)) X(t)\\
\frac{d Y(t)}{dt} = v_2 - v_3 = (k_2 X(t) - k_3) Y(t)
}
以下、Aを系に一定速度で追加し続ける定常状態を考えます。A(t)は定数となるので次のように文字で置き換えましょう。
\displaylines{
\frac{d X(t)}{dt} = (p_1 - p_2 Y(t)) X(t)\\
\frac{d Y(t)}{dt} = (p_2 X(t) - p_3) Y(t)
}
この微分方程式は残念ながら解けないみたいです。数値解を計算していきます。
・数値解(オイラー法)
オイラー法では$\frac{d y(t)}{dt}=x$を$y(t+h) \approx y(t) + x*h$と近似して計算します。$X(t)$、$Y(t)$は次のようになります。
\displaylines{
X(t+h) = X(t) + (p_1 - p_2 Y(t)) X(t)*h \\
Y(t+h) = Y(t) + (p_2 X(t) - p_3) Y(t)*h
}
import numpy as np
import matplotlib.pyplot as plt
t_max, h = 20.0, 0.01
t, x, y = 0.0, 2.0, 1.0
p1, p2, p3 = 1.0, 1.0, 1.0
tt, xx, yy = [], [], []
while t <= t_max:
tt.append(t)
xx.append(x)
yy.append(y)
t += h
x += (p1 - p2*y)*x*h
y += (p2*x - p3)*y*h
plt.plot(tt, xx, label='[X]')
plt.plot(tt, yy, label='[Y]')
plt.legend(fontsize=14)
plt.ylabel("consentration", fontsize=16)
plt.xlabel("time", fontsize=16)
plt.show()
XとYの濃度の時間変化をプロットしました。周期的に変化していることがわかります。次に、Xの濃度-Yの濃度のグラフをプロットしてみます。ついでにXの初期濃度を変化させています。
import numpy as np
import matplotlib.pyplot as plt
t_max, h = 20.0, 0.01
p1, p2, p3 = 1.0, 1.0, 1.0
x_list = [0.1, 2, 5 ,10]
for i in range(len(x_list)):
t, x, y = 0.0, x_list[i], 1.0
tt, xx, yy = [], [], []
while t <= t_max:
tt.append(t)
xx.append(x)
yy.append(y)
t += h
x += (p1 - p2*y)*x*h
y += (p2*x - p3)*y*h
plt.plot(xx, yy, label=f'X(0)={xx[0]}, Y(0)={yy[0]}')
plt.legend(fontsize=14)
plt.ylabel("[Y]", fontsize=16)
plt.xlabel("[X]", fontsize=16)
plt.show()
XとYの濃度は周期的に変化しているので円を描いて元の値に戻ります。つまり、プロットのような閉曲面となります。Xの初期濃度が変わっても閉曲面を描きますが、円の形や大きさが変化します。
ブリュッセル振動子
\displaylines{
A \longrightarrow X \\
2X + Y \longrightarrow 3X \\
B + X \longrightarrow Y + C \\
X \longrightarrow D \\
}
反応全体としては、AとBからCとDできる反応です。ここで、2段階目の反応は自己触媒反応です。
1、2、3、4段階目の反応速度定数を$k_1$、$k_2$、$k_3$、$k_4$とする次が求まります。
\displaylines{
\frac{d X(t)}{dt} = k_1 A(t) + k_2 X(t)^2 Y(t) - k_3 B(t) X(t) - k_4 X(t)\\
\frac{d Y(t)}{dt} = - k_2 X(t)^2 Y(t) + k_3 B(t) X(t)
}
ここで$A(t)$、$B(t)$、$X(t)$、$Y(t)$は時刻$t$でのAとBとXとYの濃度です。
以下、AとBを系に一定速度で追加し続ける定常状態を考え、A(t)とB(t)が定数となるので、次のように文字で置き換えます。
\displaylines{
\frac{d X(t)}{dt} = p_1 + p_2 X(t)^2 Y(t) - p_3 X(t) - p_4 X(t)\\
\frac{d Y(t)}{dt} = - p_2 X(t)^2 Y(t) + p_3 X(t)
}
・数値解(オイラー法)
\displaylines{
X(t+h) = X(t) + (p_1 + p_2 X(t)^2 Y(t) - p_3 X(t) - p_4 X(t))*h \\
Y(t+h) = Y(t) + (- p_2 X(t)^2 Y(t) + p_3 X(t))*h
}
import numpy as np
import matplotlib.pyplot as plt
t_max, h = 20.0, 0.01
p3_list = [0.5, 2, 3]
for i in range(len(p3_list)):
t, x, y = 0.0, 1.0, 1.0
p1, p2, p3, p4 = 1.0, 1.0, p3_list[i], 1.0
tt, xx, yy = [], [], []
while t <= t_max:
tt.append(t)
xx.append(x)
yy.append(y)
t += h
x += (p1 + p2*(x**2)*y - p3*x - p4*x)*h
y += (-p2*(x**2)*y + p3*x)*h
plt.plot(tt, xx, label=f'[X], p3 = {p3_list[i]}')
plt.plot(tt, yy, label=f'[Y], p3 = {p3_list[i]}')
plt.legend(fontsize=14, loc='upper left', bbox_to_anchor=(1, 1))
plt.ylabel("consentration", fontsize=16)
plt.xlabel("time", fontsize=16)
plt.show()
$p3$の値を色々変化させてプロットしてみます。$p3$が変化したというのは、3段階目の反応の反応速度が変化したか、Bの初期濃度が変化したことに対応します。$p3=0.5$の時は振動してないです。$p3=2$では振動しているように見えますが実際は減衰して行きます。$p3=3$ではしっかり振動してます。
次に、反応速度やAやBの濃度を固定してXの初期濃度を変更してみます。
import numpy as np
import matplotlib.pyplot as plt
t_max, h = 20.0, 0.01
x_list = [0.01, 1, 9]
fig, ax = plt.subplots(1,2,figsize=(15,5))
for i in range(len(x_list)):
t, x, y = 0.0, x_list[i], 1.0
p1, p2, p3, p4 = 1.0, 1, 3, 1.0
tt, xx, yy = [], [], []
while t <= t_max:
tt.append(t)
xx.append(x)
yy.append(y)
t += h
x += (p1 + p2*(x**2)*y - p3*x - p4*x)*h
y += (-p2*(x**2)*y + p3*x)*h
ax[0].plot(tt, xx, label=f'[X], X(0) = {x_list[i]}')
ax[0].plot(tt, yy, label=f'[Y], X(0) = {x_list[i]}')
ax[1].plot(xx, yy, label=f'X(0) = {x_list[i]}')
ax[0].legend(fontsize=14)
ax[0].set_ylabel("consentration", fontsize=16)
ax[0].set_xlabel("time", fontsize=16)
ax[1].legend(fontsize=14)
ax[1].set_ylabel("[Y]", fontsize=16)
ax[1].set_xlabel("[X]", fontsize=16)
plt.show()
Xの濃度-Yの濃度のグラフを見るとわかるやすいのですが、Xの初期濃度が違っていても最終的に同じ閉曲面を描くようになります。この落ち着く先の共通のサイクルを究極サイクルと言います。
オレゴン振動子
最後にオレゴン振動子をシミュレーションします。今まで具体例を出してませんでしたが、オレゴン振動子はBelousov-Zhabotinsky反応をモデル化したものです。Belousov-Zhabotinsky反応は酸性溶液中で臭素酸カリウム、マロン酸、セリウム(IV)塩を混ぜることによって進行する振動反応です。反応式を次に示します。
\displaylines{
\rm{5CH_2(COOH)_2 + 3 BrO_3^- + 3 H^+ \longrightarrow 3CHBr(COOH)_2 + 2HCOOH + 4CO_2 + 5H_2O}
}
セリウム(IV)イオンは触媒として働くので反応式には出てきません。Belousov-Zhabotinsky反応は18個の素課程と21種の化学種が関与する恐ろしく複雑なものですが、その挙動を再現しつつ五つの反応に簡略化したモデルが提案されています1。これが次のオレゴン振動子です。
\displaylines{
A + Y \longrightarrow X \\
X + Y \longrightarrow C \\
A + X \longrightarrow 2X + Z \\
2X \longrightarrow D \\
Z \longrightarrow f Y \\
}
実際のBelousov-Zhabotinsky反応ではAが$\rm{BrO_3^-}$、Cが$\rm{HBr}$、Dが$\rm{HOBr}$や$\rm{BrCH(COOH)_2}$などの生成物、Xが$\rm{HBrO_2}$、Yが$\rm{Br^-}$、Zが$Ce^{4+}$に対応します。モデルなので各文字に化学種を代入しても反応の前後で原子の数等は合いません。反応全体としては、AからCとDできる反応です。ここで、3段階目の反応は自己触媒反応です。
1、2、3、4、5段階目の反応速度定数を$k_1$、$k_2$、$k_3$、$k_4$、$k_5$とする次が求まります。
\displaylines{
\frac{d X(t)}{dt} = k_1 A(t) Y(t) - k_2 X(t) Y(t) + k_3 A(t) X(t) - 2 k_4 X(t)^2\\
\frac{d Y(t)}{dt} = - k_1 A(t) Y(t) - k_2 X(t) Y(t) + k_5 f Z(t) \\
\frac{d Z(t)}{dt} = - k_3 A(t) X(t) - k_5 Z(t)
}
ここで$A(t)$、$X(t)$、$Y(t)$、$Z(t)$は時刻$t$でのAとXとYとZの濃度です。
以下、Aを系に一定速度で追加し続ける定常状態を考え、A(t)が定数とします。ここからかなり複雑ですが、$x(t)=\frac{2 k_4}{k_3 A} X(t)$、$y(t)=\frac{k_2}{k_3 A} Y(t)$、$z(t)=\frac{k_4 k_5}{(k_3 A)^2} Z(t)$、$\epsilon_1 = k_3 A$、$\epsilon_2 = \frac{k_2 k_3 A}{2 k_4}$と置いてやって式を変形します2。こうすることで変数が減っていくばくかましになります。
\displaylines{
\frac{d x(t)}{dt} = \epsilon_1 (q y(t) - x(t) y(t) + x(t)(1 - x(t)))\\
\frac{d y(t)}{dt} = \epsilon_2 (- q y(t) - x(t) y(t) + f z(t)) \\
\frac{d z(t)}{dt} = x(t) - z(t)
}
・数値解(オイラー法)
\displaylines{
X(t+h) = X(t) + \epsilon_1 (q y(t) - x(t) y(t) + x(t)(1 - x(t)))*h \\
Y(t+h) = Y(t) + \epsilon_2 (- q y(t) - x(t) y(t) + f z(t))*h \\
Z(t+h) = Z(t) + (x(t) - z(t))*h \\
}
import numpy as np
import matplotlib.pyplot as plt
t_max, h = 30.0, 1e-5
t, x, y, z = 0.0, 1.0, 1.0, 1.0
e1, e2, q, f = 1/0.1, 1/0.0004, 0.0008, 1.0
tt, xx, yy, zz = [], [], [], []
fig, ax = plt.subplots(3,1,figsize=(8,15))
while t <= t_max:
tt.append(t)
xx.append(x)
yy.append(y)
zz.append(z)
t += h
x += e1*(q*y - x*y + x*(1 - x) )*h
y += e2*(-q*y - x*y + f*z )*h
z += (x - z)*h
ax[0].plot(tt, xx, label='[X]')
ax[1].plot(tt, yy, label='[Y]')
ax[2].plot(tt, zz, label='[Z]')
#plt.plot(xx,yy)
ax[0].legend(fontsize=14)
ax[0].set_ylabel("consentration", fontsize=16)
ax[1].legend(fontsize=14)
ax[1].set_ylabel("consentration", fontsize=16)
ax[2].legend(fontsize=14)
ax[2].set_ylabel("consentration", fontsize=16)
ax[2].set_xlabel("time", fontsize=16)
plt.show()
まとめ
今回は振動反応をシミュレーションしてみました、実際にプロットして見ることで物質濃度が時間に対して周期的に変化していることがわかりやすいですね。振動反応は生体内でも起こっており、グルコースからATPが生じる解糖系がその代表例です。また、心臓の鼓動のリズムを維持しているもの振動反応です。ただ面白いだけでなく生物にとって、とっても重要な反応です。
連載記事目次
・「Pythonによる化学シミュレーション入門」をやってみたよ~単分子一次反応の反応速度論~
・「Pythonによる化学シミュレーション入門」をやってみたよ2~色々な反応の反応速度~
-
『アトキンス 物理化学』P.W.Atkins著、千原秀昭・中村亘男訳、東京化学同人 ↩