血中濃度データに1-コンパートメント急速静注モデルを当てはめる

急速静注した薬物の5点の血中濃度推移のデータが得られたとして*1、これを1-コンパートメント急速静モデルに当てはめ、消失速度定数、半減期、分布容積、クリアランスの薬物動態パラメータを求めてみようと思います。

1-コンパートメント急速静注モデル
回帰直線の傾きと切片を求める
求めた傾きと切片から薬物動態パラメータを求める
参考

f:id:enokisaute:20200307124727p:plain

1-コンパートメント急速静注モデル

急速静注（吸収過程がない）の1-コンパートメントモデルにおける薬物量の変化の微分方程式を組み立て、これを解くと以下の式になります。

$\displaystyle{ log_{10}C = \frac{-k_e}{2.303} \cdot t + log_{10}C_0}$

$C = C_0 \cdot e^{-k_e \cdot t}$

※ $C_0$ ：初期血中濃度、 $t$ ：時間、 $k_e$ ：消失速度定数

これらは前回記事「0次反応と1次反応のグラフをPythonで描く - 薬剤師のプログラミング学習日記」で出てきた1次反応の速度式と同じものです。（上の式は常用対数にしたもの、下の式は指数表示にしたもの）
また、 $C_0(初期血中濃度) = X_0(初期投与量) / Vd(分布容積)$ です。

回帰直線の傾きと切片を求める

血中濃度データを線形回帰することによって、直線y = ax + bのパラメータa（傾き）, b（切片）を求めます。ただし、yは常用対数をとった値です。
上の式はすでにy = ax + bの形になっているので、傾きからはkeが、切片からはVdがわかります。ではやっていきます。

import numpy as np
import matplotlib.pyplot as plt
 
# 5点の血中濃度(μg/mL)データ
xdata = np.array([0.01, 0.5, 1.0, 1.5, 3.5])      # 時間
ydata = np.array([6.1, 4.7, 3.9, 2.8, 1.2])       # 血中濃度
 
# 最小二乗法により回帰直線の傾きと切片を求める
X = np.zeros((len(xdata), 2), float)    # [データ数]行×2列の0行列を作る
X[:, 0] = 1　　　　　　  # 切片項を追加. すべての行の1列目(0から開始)に1を代入
X[:, 1] = xdata         # 2列目にxdataを代入
logy = np.log10(ydata)  # 常用対数をとる
(theta, residuals, rank, s) = np.linalg.lstsq(X, logy, rcond=None)
 
intercept = theta[0, ]    # 切片
slope = theta[1, ]        # 傾き

まず5点のデータについて。xdataとydataのnp.array()内のリスト[]の数字は並び順にそれぞれ対応しています。(x, y) = (0.01, 6.1), (0.5, 4.7), ・・・(3.5, 1.2)です。
次に、回帰直線のパラメータa,bを求める部分はnumpyのlinalg.lstsq関数を使います。関数に渡す行列Xをつくる部分が少しわかりにくいかもしれませんが、コメントを参照してください。Xは1列目はすべて1, 2列目にxdataが入った次のような形をしています。

>>> print(X)
[[1.   0.01]
 [1.   0.5 ]
 [1.   1.  ]
 [1.   1.5 ]
 [1.   3.5 ]]

また、ydataは関数に渡す際に常用対数に変換しています。
np.linalg.lstsq()のところは戻り値が4つありますが、今回必要なのはthetaだけです。私の環境ではthetaは[ 0.07780118 -0.20203139]となりました。
前から順にintercept（切片）、slope（傾き）として値を取り出しています。
線形回帰についてはここでは詳しく触れませんが、こちらの記事で取り上げて書いてみました。
www.yakupro.info

求めた傾きと切片から薬物動態パラメータを求める

続いて、ke, vdを求めます。

x0 = 100  # 初期投与量(mg)
 
k = slope * -2.303                # 傾きからkを計算
vd = x0 / (10 ** intercept)       # 切片からvdを計算
 

print('傾き: ', slope)
print('切片: ', intercept)
print('消失速度定数(hr^-1): ', k)
print('半減期(hr): ', np.log(2) / k)
print('分布容積(L): ', vd)
print('クリアランス(mL/min):', k * vd * 1000 / 60)

ここで得られたinterceptは常用対数値で、 $log_{10}C_0 =$ interceptです。
$C_0$ をx0 / vdに置き換えてvdを求めています。また「**」はべき乗を計算する演算子です。
半減期とクリアランスは公式通りです。クリアランスの *1000 / 60のところはkの単位が $hr^{-1}$ 、vdはLなのでmL/minに変換しているだけです。

傾き:  -0.20203138695229084
切片:  0.7780111812844477
消失速度定数(hr^-1):  0.46527828415112576
半減期(hr):  1.4897475428593312
分布容積(L):  16.67204288392156
クリアランス(mL/min): 129.28565843875018

参考

はじめての薬物速度論-薬物動態の基礎

最後に全コードを載せておきます。

import numpy as np
import matplotlib.pyplot as plt
  
# 5点の血中濃度(μg/mL)データ
xdata = np.array([0.01, 0.5, 1.0, 1.5, 3.5])     # 時間
ydata = np.array([6.1, 4.7, 3.9, 2.8, 1.2])      # 血中濃度

 
# 最小二乗法により回帰直線の傾きと切片を求める
X = np.zeros((len(xdata), 2), float)    # [データ数]行×2列の0行列を作る
X[:, 0] = 1         # 切片項を追加. すべての行の1列目(0から開始)に1を代入
X[:, 1] = xdata     # 2列目にxdataを代入
print(X)
logy = np.log10(ydata)  # 血中濃度の常用対数をとる
(theta, residuals, rank, s) = np.linalg.lstsq(X, logy, rcond=None)
 
x0 = 100                    # 初期投与量(mg)
intercept = theta[0, ]      # 切片
slope = theta[1, ]          # 傾き
 
k = slope * -2.303                  # 傾きからkを計算
vd = x0 / (10 ** intercept)         # 切片からvdを計算
 
print('傾き: ', slope)
print('切片: ', intercept)
print('消失速度定数(hr^-1): ', k)
print('半減期(hr): ', np.log(2) / k)
print('分布容積(L): ', vd)
print('クリアランス(mL/min):', k * vd * 1000 / 60)
 
plt.figure(figsize=(10, 6))
plt.subplots_adjust(wspace=0.35, hspace=0.35)
 
x = np.arange(0, 4.6, 0.1)  # plotのため数値の配列を作成. (開始, 終値, step)
 
# 片対数プロット(左側)
plt.subplot(1, 2, 1)    # (行数, 列数, 何番目のプロットか)
# プロットは生の値をlogスケールで表示する
y1 = np.power(10, theta[0] + theta[1] * x)
 
plt.plot(x, y1, 'r-', label='$log(X0/Vd) - kt / 2.303$')
plt.ylabel('conc[ug/mL] (log scale)')
plt.yscale('log', basey=10)
plt.ylim(1, 8)
 
# ノーマルプロット(右側)
plt.subplot(1, 2, 2)
y2 = x0 / vd * np.exp(-k * x)
plt.plot(x, y2, 'r-', label='$Co * exp(-kt)$')
 
t_half = np.log(2) / k
c2 = x0 / vd * np.exp(-k * t_half)          # 半減期後のC
c3 = x0 / vd * np.exp(-k * t_half * 2)      # (半減期×2)後のC
 
# 半減期毎に補助線を入れる
# 1半減期後まで
plt.vlines(x=t_half, ymin=0, ymax=c2, color='gray', linestyles='dashed')
plt.hlines(y=c2, xmin=0, xmax=t_half, color='gray', linestyles='dashed')
# 2半減期後まで
plt.vlines(x=t_half * 2, ymin=0, ymax=c3, color='gray', linestyles='dashed')
plt.hlines(y=c3, xmin=0, xmax=t_half * 2, color='gray', linestyles='dashed')
plt.ylabel('conc[ug/mL]')
plt.ylim(0, 8)
 
# 座標軸の一覧取得
axs = plt.gcf().get_axes()
for ax in axs:
    # current axsを設定
    plt.sca(ax)
    # 測定データをプロット
    plt.plot(xdata, ydata, 'b.', label='data', alpha=0.6, markersize=8)
    # 横軸ラベルを表示
    plt.xlabel('time[hr]')
    # 横軸の範囲を指定する
    plt.xlim(0, 4.5)
    # グリッドの表示
    plt.grid(which='both', linestyle='dashed')
    # ラベルに応じた凡例を表示
    plt.legend()
 
plt.show()

*1:適当な数値です。