自宅でできる化学②～化学反応解析～

飛べない鳥

１．まえがき

　自宅でできる化学シリーズ２番目の記事になります。2Sセメスターの部内の自主ゼミでした、Pythonを使って化学反応を解析する話を記事化してみました。化学反応でどのように反応物と生成物の濃度が変化するかは、実はパソコン上で簡単にシミュレートできます。Pythonのコードの詳細な説明をする時間も気力もないので、そこは読者のみなさんが補完してくれると幸いです。別に専門的なPythonの知識は必要なく、少しググって勉強すれば分かる程度のものしか使っていないので、問題はないと思います。

２．化学反応論の基礎

　たとえば、
\[ \rm{A + B \rightarrow C} \] という化学反応について考えてみましょう。このときCの生成速度\(v\)は
\[ v=\frac{d\rm C}{dt}=k\rm{[A][B]} \]
と表されます。（\(k\)は反応速度定数、\([\rm X]\)は物質Xの濃度）高校で化学を学んだ方であればご存知でしょう。主な化学反応はこのような１階常微分方程式を解くことで解析が可能です。上の例のように複数の物質が関与する化学反応であれば、連立１階常微分方程式を解くことになります。
　まずは単純な一次反応
\[ \rm{A \rightarrow B} \]
を解析してみましょう。\(k\)を反応速度定数として、
\[ -\frac{d [\rm A]}{dt}=\frac{d [\rm B]}{dt}=k[\rm A] \]
という微分方程式が立ちます。時間が\(0 \rightarrow t\)に変化したとき、Aの濃度が\(\rm{[A]_0 \rightarrow [A]}\)に変化したとします。このとき、
\[ -\frac{1}{[\rm A]}\frac{d [\rm A]}{dt}=k \\ -\int_0^t\frac{1}{[\rm A]}\frac{d [\rm A]}{dt}dt=k\int_0^tdt \\ -\int_{[\rm A]_0}^{[\rm A]}\frac{1}{[\rm A]}d [{\rm A}]=k\int_0^t dt \\ {\rm ln[A]-ln[A]_0}=-kt \\ [{\rm A}]=[{\rm A}]_0{\rm exp}(-kt) \]
となります。このように微分方程式を解析的に解いたときの濃度変化のグラフは以下のようになります。（\([\rm A]_0=1.0\)としています）

　一応このグラフを出力するプログラムを下に載せておきます。

    
      import math
      import copy
      import matplotlib.pyplot as plt
      from tqdm.notebook import tqdm

      t = []
      a = []
      b = []

      t1 = 0 # 時間の初期値
      a1 = 1 # Aの初濃度
      b1 = 0 # Bの初濃度
      dt = 0.1 # 刻み幅
      tm = 100 # 最大時刻

      k = 0.1 # 反応速度定数

      pbar = tqdm(total = tm / dt) # プログラムの進捗を表現
      count = 0

      a0 = copy.copy(a1)

      while t1 <= tm:
          if count % 10 == 0: # プロットする点の数を制限
              a.append(a1)
              b.append(b1)
              t.append(t1)
          t1 += dt
          a1 = a0*math.exp(-k*t1)
          b1 = a0-a1
          count += 1
          pbar.update(1)
          
      pbar.close()
      print("初期条件：[A]0 = " + str(c0))
      plt.plot(t, a, color="r", label="[A]")
      plt.plot(t, b, color="b", label="[B]")
      plt.xlabel("t", fontsize="16")
      plt.ylabel("c", fontsize="16")
      plt.legend()
      plt.show()

　じゃあ、常微分方程式の勉強をすれば化学反応が全部わかるじゃん！　とは残念ながらならないんです。解析的に解ける常微分方程式の形は限られており、解けない常微分方程式はいっぱいあります。世の中の常微分方程式を２種類に分けると、以下のような形をした「線形常微分方程式」
\[ \frac{d^{n}}{d x^{n}} y+p_{1}(x) \frac{d^{n-1}}{d x^{n-1}} y+p_{2}(x) \frac{d^{n-2}}{d x^{n-2}} y+\cdots+p_{n}(x) y=q(x) \]
と、上のように表せない「非線形常微分方程式」に分かれます。線形常微分方程式は解析的に解けることがわかっていますが、非線形常微分方程式は一般に解析的に解くことができません。そこで、多少の誤差を許容しつつもあらゆる常微分方程式に通用する解法が求められるわけです。

３．微分方程式の数値的解法

　ここでは以下の３つの数値的解法を紹介します。
①　Euler法
②　（4次の）Runge-Kutta法
③　Runge-Kutta-Fehlberg法

①　Euler法
　まずは一番簡単なEuler法から始めていきましょう。関数\(f(x)\)の微分の定義は高校で以下のように習うと思います。
\[ \lim_{h\rightarrow 0}\frac{f(t+h)-f(t)}{h}=\frac{df}{dt} \]
ここで、\(h\)が十分小さいならば、 \[ \frac{f(t+h)-f(t)}{h}\approx\frac{df}{dt} \] すなわち、 \[ f(t+h)\approx\frac{df}{dt}h+f(t) \] と表せるはずです。
　初期値\(t_0\)と、刻み幅\(h\)、増分\(df/dt\)を与えたとき、 \[ f(t_0+h)\approx\frac{df}{dt}h+f(t_0) \]

\[ f(t_0+2h)\approx\frac{df}{dt}h+f(t_0+h) \]

\[ \vdots \]
のようになり、漸化式を解くようにして常微分方程式が解けることになります。
　先程の一次反応をEuler法で解くプログラムを作りました。

    
      import math
      import copy
      import matplotlib.pyplot as plt
      from tqdm.notebook import tqdm

      aa = []
      bb = []
      t = []

      t1 = 0 # 時間の初期値
      a2 = 1 # Aの初濃度
      b2 = 0 # Bの初濃度
      dt = 0.1 # 刻み幅
      tm = 100 # 最大時刻

      k = 0.1 # 反応速度定数

      a0 = copy.copy(a2)

      pbar = tqdm(total = tm / dt) # プログラムの進捗を表現
      count = 0

      f = lambda a2: -k*a2 # d[A]/dt

      while t1 <= tm:
          if count % 10 == 0: # プロットする点の数を制限
              aa.append(a2)
              bb.append(b2)
              t.append(t1)
          a2 += dt * f(a2)
          b2 += - dt * f(a2)
          t1 += dt
          count += 1
          pbar.update(1)

      pbar.close()
      print("初期条件：[A]0 = " + str(c0))
      plt.plot(t, aa, color="r", label="[A]")
      plt.plot(t, bb, color="b", label="[B]")
      plt.xlabel("t", fontsize="16")
      plt.ylabel("c", fontsize="16")
      plt.legend()
      plt.show()

　すると以下のようなグラフを得ました。

　解析的に解いたときとグラフの形が変わっていないように見えます。どれほど誤差が少ないかそれぞれのグラフを重ねて確認してみましょう。

    
      plt.scatter(t, a, color="k", label="analytical")
      plt.plot(t, aa, color="r", label="Euler method")
      plt.xlabel("t", fontsize="16")
      plt.ylabel("c", fontsize="16")
      plt.legend()
      plt.show()

　上のプログラムを走らせるとこのようなグラフを得ました。

こうして見ると、このケースにおいてEuler法はかなり正確に結果を再現できていることがわかります。

②　（4次の）Runge-Kutta法
　①のEuler法は精度が十分なように見えます。ここで先程でてきた近似式 \[ f(t+h)\approx\frac{df}{dt}h+f(t) \] について考えてみましょう。\(f(t+h)\)をTaylor展開すると、
\[ f(t+h)=f(t)+f^{\prime}(t) h+\frac{1}{2} f^{\prime \prime}(t) h^{2}+\cdots \] となり、Euler法は第３項以降を無視し、第２項までで打ち切っていることがわかります。これを「\(O(h^2)\)の精度」と表現します。（”\(O\)”はランダウの記号と呼ばれるものです）これだと、時間発展するたびに無視された第３項以降の誤差が蓄積してズレが大きくなってしまいます。先程の例であれば問題ありませんでしたが、より複雑な問題に対しては\(O(h^2)\)の精度では不十分です。そこでRunge-Kutta法が出てきます。このRunge-Kutta法ではなんと\(O(h^5)\)の精度が実現できます。計算方法は以下の通り。

　まず次のような微分方程式を考えます。 \[ \frac{dx}{dt}=f(t, x) \] 初期値として\(t=t_0,\ x=x_0\)を与え、\(t=t_0+h\)のときの\(x=x_1\)を求めたいです。このとき、（4次の）Runge-Kutta法によれば、 \begin{eqnarray*} k_1 &=& h×f\left(t_0, x_0\right) \\ k_2 &=& h×f\left(t_0+\frac{h}{2}, x_0+\frac{k_1}{2}\right) \\ k_3 &=& h×f\left(t_0+\frac{h}{2}, x_0+\frac{k_2}{2}\right) \\ k_4 &=& h×f\left(t_0+h, x_0+k_3\right) \end{eqnarray*} を計算し、これらの加重平均 \[ k=\frac{k_1+2k_2+2k_3+k_4}{6} \] が、時間が\(t\)から\(t+h\)に進んだときの\(x_0\)の変化量を表し、\(x_1=x_0+k\)となります。この計算を繰り返して\(x_2,x_3,x_4,\dots\)を求めていくことで常微分方程式を数値的に解いたことになります。これらの式をうまく変形すればTaylor展開したときに第5項以降を落とした式が出てきます。式変形を記述するのは大変なので、ここでは直感的な説明にとどめます。（ちなみにカッコ書きで（４次の）とつけているのは４次以外のものがあるからです）
直感的には、Runge-Kutta法の計算は以下のように理解できます。

まずEuler法同様に、幅を半分にして傾き\(k_1\)を算出 \[ k_1 = h×f\left(t_0, x_0\right) \]
真ん中の値を\(k_1\)から推定。その後この値を初期値として、\(h/2～h\)までの間の傾き\(k_2\)を計算 \[ k_2 = h×f\left(t_0+\frac{h}{2}, x_0+\frac{k_1}{2}\right) \]
今度は、\(0～h/2\)での傾きが\(k_2\)だったとして、真ん中の値を推定。その上で、この値を初期値にして、\(h/2～h\)までの間の関数の傾き\(k_3\)を計算 \[ k_3 = h×f\left(t_0+\frac{h}{2}, x_0+\frac{k_2}{2}\right) \]
今度は、\(0～h\)での傾きが\(k_3\)だったとして、傾き\(k_4\)を推定 \[ k_4 = h×f\left(t_0+h, x_0+k_3\right) \]
最後に、これらの加重平均 \[ k=\frac{k_1+2k_2+2k_3+k_4}{6} \] を求めて完了！

③　（６段５次の）Runge-Kutta-Fehlberg法
　先述のEuler法・Runge-Kutta法では、刻み幅\(h\)を自分で設定する必要がありました。このとき自分が設定した刻み幅\(h\)で十分な精度が実現できているかどうかは、刻み幅を細かくしたもの（たとえば、\(h/2\)）をあらかじめ計算して、比較する必要があります。しかし、これでは計算量が膨大になってしまうのが良くないところ。そこでRunge-Kutta-Fehlberg法です。Runge-Kutta-Fehlberg法では時間の刻み幅を変更しながら計算するため、Runge-Kutta法より高い\(O(h^6)\)の精度が実現できます。計算方法は以下の通り。

　次のような微分方程式を考えます。 \[ \frac{dx}{dt}=f(t, x) \] 初期値として\(t=t_0,\ x=x_0\)を与え、\(t=t_0+h\)のときの\(x=x_1\)を求めたいです。（６段５次の）Runge-Kutta-Fehlberg法によれば、 \begin{eqnarray*} k_{1} &=& h \times f\left(t_{0}, x_{0}\right) \\ k_{2} &=& h \times f\left(t_{0}+\frac{1}{4} h, x_{0}+\frac{1}{4} k_{1}\right) \\ k_{3} &=& h \times f\left(t_{0}+\frac{3}{8} h, x_{0}+\frac{3}{32} k_{1}+\frac{9}{32} k_{2}\right) \\ k_{4} &=& h \times f\left(t_{0}+\frac{12}{13} h, x_{0}+\frac{1932}{2197} k_{1}-\frac{7200}{2197} k_{2}+\frac{7296}{2197} k_{3}\right) \\ k_{5} &=& h \times f\left(t_{0}+h, x_{0}+\frac{439}{216} k_{1}-8 k_{2}+\frac{3680}{513} k_{3}-\frac{845}{4104} k_{4}\right) \\ k_{6} &=& h \times f\left(t_{0}+\frac{1}{2} h, x_{0}-\frac{8}{27} k_{1}+2 k_{2}-\frac{3544}{2565} k_{3}+\frac{1859}{4104} k_{4}-\frac{11}{40} k_{5}\right) \end{eqnarray*} を計算すると、４次の近似\(x_{1}\) は、 \[ x_{1}=x_{0}+\frac{25}{216} k_{1}+\frac{1408}{2565} k_{3}+\frac{2197}{4104} k_{4}-\frac{1}{5} k_{5} \] ５次の近似 \(x_{1}^{\prime}\) は、 \[ x_{1}^{\prime}=x_{0}+\frac{16}{135} k_{1}+\frac{6656}{12825} k_{3}+\frac{28561}{56430} k_{4}-\frac{9}{50} k_{5}+\frac{2}{55} k_{6} \] で求められます。このとき誤差Rを、 \[ R=\left|x_1^{\prime}-x_1\right| \] とします。
　次に、\(R\)が許容範囲内であるか確認します。許容範囲内であるかどうかは以下の式を計算して求めます。 \[ s=\left(\frac{εh}{2R}\right)^{\frac{1}{4}} \] 　\(\varepsilon\)は許容誤差であり、 \(R>\varepsilon\)ならば、刻み幅を\(h \rightarrow sh\)に修正して再度計算し直します。\(R \leq \varepsilon\)ならば、\(x_1\)が、時間が\(t+h\)のときの\(x\)の値となります。どうやら実際のプログラム上では、刻み幅の修正の仕方をここまで単純にはしないようですので詳しくはこちらをご参照下さい。

４．実際の化学反応の解析

　解法の説明に随分と時間を掛けてしまいました。最後に、Sel’kovモデルという解糖系の反応モデルに対して今回の解法を使ってみることにしましょう。なお当方は生物弱者ですので、解糖系の説明はWikipediaの説明を~~コピペ~~引用することで対応しようと思います。

解糖系（かいとうけい、Glycolysis）とは、生体内に存在する生化学反応経路の名称であり、グルコースをピルビン酸などの有機酸に分解（異化）し、グルコースに含まれる高い結合エネルギーを生物が使いやすい形に変換していくための代謝過程である。ほとんど全ての生物が解糖系を持っており、もっとも原始的な代謝系とされている。（Wikipediaの項目：「解糖系」より引用）

右枠内の反応が、解糖系の律速段階になります。化学反応式で表すと、

F6P + ATP → F1,6BP + ADP

となり、6-ホスホフルクトキナーゼという酵素によって触媒されています。
解糖系での濃度変化がSel’kovモデルで記述される様子をグラフで見てみましょう。
　Sel’kovモデルの中身は以下の通りです。
　ADP、F6Pの濃度をそれぞれ\(x,\ y\)とします。このとき、
①　F6Pは一定の速さbで供給される
②　反応で生成するADPは6-ホスホフルクトキナーゼを活性化する（反応速度：\(ay+x^2 y\)、\(a\)は定数）
③　ADPは一次で分解していく
という条件があります。①、②、③を微分方程式で表現すると、 \[ \frac{dx}{dt}=(ay+x^2y)-x \\ \frac{dy}{dt}=b-(ay+x^2y) \] となります。これをシミュレーションしてみましょう。

      
        import math
        import copy
        import matplotlib.pyplot as plt
        from tqdm.notebook import tqdm

        x = []
        y = []
        t = []

        x1 = 0.5 # ADPの初濃度
        y1 = 0.5 # フルクトース－6－リン酸の初濃度
        t1 = 0 # 時間の初期値
        dt = 0.001 # 刻み幅
        tm = 50 # 最大時刻

        a = 0.06
        b = 0.6

        x_0 = copy.copy(x1)
        y_0 = copy.copy(y1)

        f = lambda x1, y1: a * y1 + x1**2 * y1 - x1 # dx/dt
        g = lambda x1, y1: b - a * y1 - x1**2 * y1 # dy/dt

        pbar = tqdm(total = tm / dt) # プログラムの進捗を表現
        count = 0

        while t1 <= tm:
            if count % 10 == 0: # プロットする点の数を制限
                x.append(x1)
                y.append(y1)
                t.append(t1)
            x1 += dt * f(x1, y1)
            y1 += dt * g(x1, y1)
            t1 += dt
            count += 1
            pbar.update(1)
            
        pbar.close()
        plt.plot(t, x, color="b", label="ADP")
        plt.plot(t, y, color="r", label="F6P")
        plt.legend(loc="best")
        plt.xlabel("t", fontsize="16")
        plt.ylabel("C", fontsize="16")
        plt.show()

　すると、次のグラフを得ました。（初期値を\((x_0,y_0)=(0.5,0.5)\)、\(a=0.06\)、\(b=0.6\)としました）

　ADPとF6Pの濃度がきれいに振動していますね。こうした周期的に濃度が変化する反応を「振動反応」と呼びます。演示実験解説で出てきた「BZ反応」もこの振動反応の１つです。もう１つ、このモデルの面白いところを見せましょう。ADPの濃度\(x\)を縦軸、F6Pの濃度\(y\)を横軸にすると……

      
        plt.plot(x, y, color="b")
        plt.xlabel("x", fontsize="16")
        plt.ylabel("y", fontsize="16")
        plt.show()

　このようなグラフを得ました。初期値\((x_0,y_0)=(0.5,0.5)\)の点から伸びた線が１つの閉曲線上に収束していく様子がわかると思います。このような閉曲線を「リミットサイクル」と呼び、今回のように適切に定数\(a,\ b\)を設定すれば、なんとどのように初期値を選んでも必ずこのリミットサイクル上をぐるぐる回るようになります。試しに初期値を\((x_0,y_0)=(1.0,1.0)\)に変えてみると、

となりました！　面白いですね～～せっかくコードを載せたので、是非読者の皆さんにも試して頂きたいです。

５．あとがき

　高校時代に振動反応の研究をしていたのですが、その中でこういう反応モデルを立ててExcelのマクロでシミュレーションすることをしていました。ここで紹介したRunge-Kutta法は実際の研究でも利用される方法ですので、知っておいて損はないと思います。あ、でもExcelでやるのはおすすめしないです。Pythonの方が圧倒的に早いです。当時はPythonを知らなかったのと、最初にExcelを使用したシミュレーションの本を紹介されたのがあってExcelを使っていました。でも計算時間が長ければPythonよりC言語とかFortranのほうが早くておすすめらしいですね？　数値計算の話は軽く知っている程度なので具体的にどれくらい差がでるか実感したことはないです。化シスの主題科目でこういうことをしたときはPythonを使っていました。時間さえ気にしなければPythonでもいいでしょう。とっつきやすいですし。プログラムは全部Euler法で書きましたが、少し書き換えるだけでRunge-Kutta法にすることができますのでそちらも試してみてはいかがでしょうか。参考文献にRunge-Kutta法のプログラムの例があるので是非。

６．参考文献リスト

（最終閲覧日は全て2020/11/19）
[1] 「オイラー法/ホイン法/ルンゲクッタ法をつかった常微分方程式の数値解析超入門 – Qiita」，
https://qiita.com/Dr_ASA/items/d44560cf1367290aee6d
[2] 「ルンゲ=クッタ法の説明と刻み幅制御 | シキノート」，
https://slpr.sakura.ne.jp/qp/runge-kutta-ex/#rkf45hautop
[3] 「解糖系 – Wikipedia」，
https://ja.wikipedia.org/wiki/解糖系
[4] 藤井雅史（2013）「生化学反応系で見られる振動現象」，
http://kurodalab.bs.s.u-tokyo.ac.jp/class/Summer/2013/Day6/day6.pdf
[5] 神足史人（2018）『Excelで操る!ここまでできる科学技術計算』丸善出版．