薬剤師のプログラミング学習日記

プログラミングやコンピュータに関する記事を書いていきます

0次反応と1次反応のグラフをPythonで描く

薬物動態でよく出てくる0次、1次の速度過程のグラフを描いてみます。

0次反応と1次反応の速度式

0次反応と1次反応における薬物量変化の収支式は以下の微分方程式で表すことができます。

0次反応: \displaystyle{ \rm \frac{dX}{dt} = -k \cdot X ^0 = -k}

1次反応: \displaystyle{ \rm \frac{dX}{dt} = -k \cdot X ^1 = -k \cdot X}

※ X:薬物量、t:時間、k:消失速度定数

それぞれの消失速度について、0次では残存薬物の量に依存せず(Xの0乗)常に一定、1次ではその時の薬物量(Xの1乗)に比例します。

次に、これらの微分方程式を解くと

0次反応: \rm X = -kt + X_0

1次反応: \rm lnX = -kt + lnX_0

X_0はt=0のときの薬物量

となります。

各反応の半減期

薬物の量や濃度が半分に減少するのに要する時間を半減期といい、上の式のXにX_0 /2を代入してtについて解くと得られます。

0次反応: \displaystyle{ \rm t_{1/2} = \frac{X_0}{2k}}

1次反応: \displaystyle{ \rm t_{1/2} = \frac{ln2}{k} = \frac{0.693}{k}}


速度式を関数として定義する

# 0次反応
def zero_order_reaction(x_0, k, t):
    return -k * t + x_0
 
# 1次反応
def first_order_reaction(x_0, k, t):
    return x_0 * np.exp(-k * t)

Pythonで関数を作るときには、defの後に関数名を書き、それに続くカッコ内に引数を書きます。また、returnを使うと関数内で処理した結果を戻り値として関数の外に返すことができます。
この場合、引数は初期薬物量(X_0)、消失速度定数(k)、時間(t)、戻り値として残存薬物量Xを返すようにしています。
Xを返すので、1次反応の速度式は指数表示(X = X_0 \cdot e^{-k \cdot t})にしています。

この関数の引数に次の値を渡してグラフを描画しました。

time = np.linspace(0, 10, 100)  # (始点, 終点, 分割数) 0~10を100分割
x0 = 100        # t=0のときの薬物量
t_half = 2.5    # 半減期(0次反応においては100→50までの時間)
 
# 各反応における消失速度定数kを計算する
k0 = x0 / (2 * t_half)
k1 = np.log(2) / t_half

消失速度定数kは各反応の半減期の式から求めています。

f:id:enokisaute:20200510192219p:plain

各反応における消失のパターン

0次反応では時間に対して直線的に減少します(=速度が一定)が、1次反応では薬物量が減ってくるにしたがって変化の割合が緩やかになっていきます(=速度は薬物量に依存)。
また、半減期については0次反応では濃度に比例しています(100→50のときは2.5、50→25のときは1.25、・・・)が、1次反応では濃度に依存せず一定(100→50は2.5、50→25は2.5、・・・)となっています。

以下が全体のコードです。プログラムの実行方法がわからないという方はこちらの記事を参考にしてみてください。
Pythonプログラミングの始め方 - 薬剤師のプログラミング学習日記


import matplotlib.pyplot as plt
import numpy as np
 
# 0次反応と1次反応のグラフを描く
 
# 0次反応
def zero_order_reaction(x_0, k, t):
    return -k * t + x_0
 
# 1次反応
def first_order_reaction(x_0, k, t):
    return x_0 * np.exp(-k * t)
 
  
time = np.linspace(0, 10, 100)    # (始点, 終点, 分割数) 0~10を100分割
x0 = 100        # t=0のときの薬物量
t_half = 2.5    # 半減期(0次反応においては100→50までの時間)
 
# 各反応における消失速度定数kを計算する
k0 = x0 / (2 * t_half)
k1 = np.log(2) / t_half
 
# x0(初期薬物量), k(消失速度定数), time(時間)を渡してyを求める
y0 = zero_order_reaction(x0, k0, time)
y1 = first_order_reaction(x0, k1, time)
 
plt.plot(time, y0, color='blue', label='0次反応')
plt.plot(time, y1, color='red', label='1次反応')
# t=2.5における濃度の点に補助線を入れる
plt.hlines(y=50, xmin=0, xmax=2.5, color='gray', linestyles='dashed')
plt.vlines(x=2.5, ymin=0, ymax=50, color='gray', linestyles='dashed')
# x,y軸の範囲を指定する
plt.xlim(0, 10)
plt.ylim(0, 100)
# 軸ラベルを表示
plt.xlabel('時間')
plt.ylabel('残存薬物量')
# 軸の目盛り間隔を指定する. x軸は2.5, y軸は10
plt.xticks(np.arange(0, 10.5, 2.5))
plt.yticks(np.arange(0, 110, 10))
# グリッドの表示
plt.grid(linestyle='dashed')
# ラベルに応じた凡例を表示
plt.legend()
# グラフを表示
plt.show()


参考文献

・徹底解説 薬物動態の数学ー微積分と対数、非線形ー
・薬物動態解析入門 はじめての薬物速度論
・みんなのPython 第3版