【python】抽象化学反応系でヨウ素時計反応をシミュレーションする

ARMSについて
ヨウ素時計反応について
コード
結果
【参考文献】

ARMS（Abstract Rewriting System on Multisets）[1]というものを用いてヨウ素時計反応（説明は後述）をシミュレートしてみます。

本来のARMSは文字列をある規則で書き換えることにより抽象的に化学反応を表すものです。

ですが、今回は他の論文で用いられていた、ベクトルの加法演算で化学反応を表現するという手法で行こうと思います。

ARMSについて

反応系にある各分子の数をベクトルの成分とし、全体の分子構成をベクトルとして記述します。

分子 $\displaystyle\rm{a, b, c}$ とあってそれぞれ1, 2, 3個なら $\displaystyle\rm{（1, 2, 3）}$ というベクトルになります。

次に反応規則についてです。ある反応 $\displaystyle\rm{a+a=b}$ で、分子 $\displaystyle\rm{a}$ が2個消費され分子 $\displaystyle\rm{b}$ が1個生成するとします。

このとき $\displaystyle\rm{a}$ は2個減り $\displaystyle\rm{b}$ は1個増えるので、ベクトルの成分をその分だけ増減させます。そうやって反応規則を表現します。

ヨウ素時計反応について

ヨウ素時計反応とは、2つの水溶液を混ぜ合わせ、最初何も起きない状態から少し経つと突然溶液が青紫色に呈色する反応です。

青紫色の呈色はヨウ素デンプン反応からきています。
ヨウ素時計反応 - Wikipedia

ヨウ素時計反応は以下のような3つの反応で成り立っています。

$\displaystyle\rm{➀IO_3^− + 3HSO_3^− → I^− + 3HSO_4^−}$

$\displaystyle\rm{➁IO_3^− + 5I^− + 6H^+ → 3I_2 + 3H_2O}$

$\displaystyle\rm{➂I_2 + HSO_3^− + H_2O → 2I^− + HSO_4^− + 2H^+ }$

一番上の反応は遅い反応で、2番目は一番遅い反応、3番目は速い反応です。

この反応速度の違いを作るために、リストの要素に定数をかけることにしました。

速い反応ほど一度にたくさんの分子が反応すると解釈したからです。

ちなみに確率を決定してやるほうが簡単ですがわざわざこうした理由はありません。プログラムを作成した後に気が付きました。

コード

そうしてできたのが以下のコードです。

modecule_listが分子数、r_1, r_2, r_3は反応規則を表します。

反応に必要な分子が足りなくなった場合にはその反応をスキップするようにしています。

また、ベクトルでやらなくてもリストでできるのでは？と思ったのでリストでやりました。

import random
import numpy as np
import matplotlib.pyplot as plt
modecule_list = ['a', 'b', 'c', 'd', 'e', 'f', 'g']

X = [1000, 1000, 0, 0, 1000, 0, 10000]
α = 3
β = 1
γ = 6

matrix = [X]

r_1 = [-1, -3, 1, 3, 0, 0, 0]

r_2 = [-1, 0, -5, 0, -6, 3, 3]

r_3 = [0, -1, 2, 1, 2, -1, -1]


for i in range(len(r_1)):
  
    r_1[i] *= α
    r_2[i] *= β
    r_3[i] *= γ
   

for i in range(500):
    j = random.randint(0,3)
    if j <= 1:
        if (X[0] < 1*α) or (X[1] < 3*α):
            continue
        X = [x + y for x, y in zip(X, r_1)]

    if 1 < j <= 2:
        if (X[0] < 1*β) or (X[2] < 5*β) or (X[4] < 6*β):
            continue
        X = [x + y for x, y in zip(X, r_2)]

    if 2 < j <= 3:
        if (X[1] < 1*γ) or (X[5] < 1*γ) or (X[6] < 1*γ):
            continue
        X = [x + y for x, y in zip(X, r_3)]
    
    matrix += [X]
    if X[1] <= 0:
        continue

a_numbers = []
b_numbers = []
c_numbers = []
d_numbers = []
e_numbers = []
f_numbers = []

for i in range(len(matrix)):
    a_numbers += [matrix[i][0]]
    b_numbers += [matrix[i][1]]
    c_numbers += [matrix[i][2]]
    d_numbers += [matrix[i][3]]
    e_numbers += [matrix[i][4]]
    f_numbers += [matrix[i][5]]

plt.plot(a_numbers, label = 'IO3^-')
plt.plot(b_numbers, label = 'HSO3^-')
plt.plot(c_numbers, label = 'I^-')
plt.plot(d_numbers, label = 'HSO4^-')
plt.plot(e_numbers, label = 'H^+')
plt.plot(f_numbers, label = 'I2')
plt.title('iodine clock reaction')
plt.xlabel('steps')
plt.legend()
plt.show()