Pythonで数理モデルシミュレーション

2021年2月6日 15:58

Noteなるものを初投稿してみます．

今日はpythonを使って数理モデルシミュレーションをしてみました．

数理モデルってのは，生物の挙動を数式によって表現して，将来的にどのように変化していくのかとか，その変化を制御したらどうしたらいいのかといったことを解析するためのツールです（説明下手ですね．．）．

今回は細胞内プロセスであるセントラルドグマ（DNA > mRNA > protein）を例にモデルを書いてシミュレーションしてみました．
MatlabというMathwaorks社が提供してるソフトが数理モデル解析によく使われているのですが，今回はpythonで実施してみました．
使用例はMatlabの例にもなっているモデルで，以下HP参照．

このモデルでは以下のプロセスが含まれています．
1. DNAからmRNAへの転写
2. mRNAからproteinへの翻訳
3. DNAとproteinとの複合体の形成
4. mRNAの分解
5. proteinの分解
最終産物のproteinがDNAと複合体を形成することで転写が抑制される，いわゆる”フィードバック制御”が作用するわけですね．

それでは早速pythonで書いていきましょう．
環境はこちらの通り．
・Windows 10
・Python 3.7.4
・PyCharm（無料版）

できあがったコードはこんな感じ．
微分方程式（ODE）の積分にはscipyのintegrate内のLSODAを使いました．
（いろいろ試してみたら，LSODAはstiffなODEでもけっこう頑張って積分してくれそう）

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt


def ode_f(t, y, k):
   """
   y[0]: DNA
   y[1]: mRNA
   y[2]: protein
   y[3]: DNA-protein complex
   """
   dy = [0]*4
   dy[0] = -k[2] * y[0] * y[2] + k[3] * y[3]
   dy[1] = k[0] * y[0] - k[4] * y[1]
   dy[2] = k[1] * y[1] -k[2] * y[0] * y[2] + k[3] * y[3] - k[5] * y[2]
   dy[3] = k[2] * y[0] * y[2] - k[3] * y[3]

   return dy

def integration(k,t_span, y0, method="LSODA", atol=1e-10, t_eval=None):
   sol = integrate.solve_ivp(lambda t, y: ode_f(t, y, k), t_span, y0,
                             method=method,
                             atol=atol,
                             t_eval=t_eval
                             )
   return sol

if __name__ == "__main__":
   t_span = (0, 10)
   
   y0 = [50.,#DNA
         0.,#mRNA
         0.,#protein
         0.#DNA-protein complex
         ]  # Initial value vector

   k = [0.2,#k1
        20.,#k2
        0.2,#k3
        1.,#k3r
        1.5,#k4
        1.#k5
        ]  # parameter vector

   sol = integration(k, t_span, y0)
   result = np.vstack([sol.t,sol.y])

   plt.plot(result[0], result[1], label = "DNA")
   plt.plot(result[0], result[2], label = "mRNA")
   plt.plot(result[0], result[3], label = "protein")
   plt.plot(result[0], result[4], label = "DNA-protein complex")
   plt.title("Gene regulation")
   plt.xlabel("Time")
   plt.ylabel("Species Amount")
   plt.legend(fontsize=7)

   plt.show()

結果は以下の通り．

MatlabのHPに載っているFigと同様の結果となりました，めでたしめでたし👏

この記事が気に入ったらサポートをしてみませんか？