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

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

1-コンパートメント点滴静注繰り返し投与ー最高・最低血中濃度とAUCー

1-コンパートメント点滴静注繰り返し投与の定常状態における最高・最低血中濃度とAUCを求めてみます。

f:id:enokisaute:20200307125912p:plain


定常状態における平均血中濃度と最高・最低血中濃度

定常状態における最高・最低血中濃度(それぞれCss_max, Css_minとします)は次の式で表せます。

\displaystyle{ C_{\rm ss\_\max} = \frac{k_0}{Vd \cdot k_e} \cdot \frac{1-e^{-k_eT}}{1-e^{-k_eτ}}}

 \displaystyle{ \begin{eqnarray*} C_{\rm ss\_\min}  & = & \frac{k_0}{Vd \cdot k_e} \cdot \frac{1-e^{-k_eT}}{1-e^{-k_eτ}} \cdot e^{-k_e(τ-T)} \\ \\ & = & C_{\rm ss\_\max} \cdot e^{-k_e(τ-T)} \end{eqnarray*}}

1-コンパートメント点滴静注繰り返し投与のグラフ描画」のt≥Tの式において、定常状態下(n=∞)と考えるとe^{-nkτ}≒0となります。Css_maxは点滴終了時(t=T)の血中濃度なので後ろの項のe^{-k_e(t-T)}は1となり、上のようになります。
また、Css_minはt=τ~~のときです。
次に、定常状態における平均血中濃度(以後Css_aveとします)ですが、これは定常状態において1投与間隔のAUC(血中濃度-時間曲線下面積)を投与間隔で割ったものとなります。

\displaystyle{ C_{\rm ss\_ave} =\frac{1}{τ} \cdot \int_0^τ C_{\rm ss}dt}

今回は前回記事
www.yakupro.info

で作成したクラスを継承して、これら3つの値を計算する機能を加えたいと思います。
Css_max, Css_minについては数式をそのままコーディングしますが、Css_aveの\int_0^τ C_{ss}dtを求める部分はscipyの数値積分によって求めます。

AUC(血中濃度-時間曲線下面積)とは

ここでAUCについて少し触れておきます。AUCは文字通り血中濃度曲線と時間軸で囲まれた部分の面積(単位は濃度×時間)のことで、その意味としては、体内に入った総薬物量と捉えられます。
定常状態に達した後の1投与間隔のAUCを長方形と考えると、投与間隔(横)×平均血中濃度(高さ)であるので、AUCがわかれば平均血中濃度が求まることになります。

ちなみに、薬物動態が線形*1である場合は1-コンパートメントモデルでも2-コンパートメントモデルでも単回投与時の区間0~∞のAUCと、繰り返し投与時の1投与間隔のAUCは同じになります。
f:id:enokisaute:20200523102439p:plain
この1-コンパートメント点滴静注モデルの単回投与時の\rm{AUC_{0-∞}}は0~点滴終了時のAUC+点滴終了時~∞のAUCで求められ、\frac{\rm{X_0}}{Vd \cdot k_e}=\frac{C_0}{k_e}となります(これは急速静注モデルでのAUCともまったく同じ)。
なので、上で書いた\int_0^τ C_{ss}dt*2は簡単に\frac{C_0}{k_e}で求めることができますが、ここはプログラミング学習ということで数値積分のコードを書いてみようと思います。

クラスを継承する

1-コンパートメント点滴静注繰り返し投与のグラフ描画」で実装したBloodConcクラスの機能を受け継いだまま、新たな機能をもったBloodConc2クラスを作成します。

クラス名横の()に親クラスの名前を書くことで継承できますが、今回私は別ファイル(repeat_div.py)として定義している親クラスを継承しているため、「モジュール名.クラス名」となっています(モジュールのインポート(下のコード1行目)が必要です)。
同じようにする人は前回書いたBloodConcクラスのコードを同一ディレクトリ内にrepeat_div.pyとして保存しましょう。なお、このファイル名は変更可能ですが、それに合わせてインポートとモジュール.クラス名の部分も適宜変更してください。

import repeat_div
 
class BloodConc2(repeat_div.BloodConc):
    def __init__(self, x0, vd, ke, T, interval):
        super().__init__(x0, vd, ke, T, interval)
        t_half = np.log(2) / ke
        # 投与間隔の倍数のうち、半減期の7倍以上の最小値を定常状態到達時間とする
        self.css_t = (t_half * 7 // interval + 1) * interval

コンストラクタ(def __init__())をオーバーライドして新たに定常状態到達時間(css_t)という変数を持たせてみることにしました。
super().__init__()を書くことで親クラスのコンストラクタをそのまま引き継ぎ、下に新たなインスタンス変数の初期化を追加しています。

数値積分でAUCを求める

scipyは高度な科学技術計算を行うためのライブラリです。Anacondaをインストールした人は別途インストールする必要なく、そのまま下のように書いて使うことができます。

簡単な例として、y=x^2を0~3の範囲で積分してみます。
f:id:enokisaute:20200307130111p:plain

なお、数式から公式を使って解いた値は9となります。

from scipy import integrate
 

s = integrate.quad(lambda x: x ** 2, 0, 3)  # (被積分関数, 積分の下限値, 上限値)

print(s[0], s[1])  # (9.000000000000002, 9.992007221626411e-14)

integrate.quad()のところで積分計算をしています。
戻り値は二つあり、積分値と推定誤差がタプルで返されます。また、quad()の引数の被積分関数のところはlamda式を使って関数を定義していますが、これは

def function(x):
    return x ** 2

という意味です。lamdaを使えばちょっとした関数を作るときに、わざわざ名前をつけることなく簡潔に書くことができます。

ではこのintegrate.quadを使って定常状態下の投与間隔区間のAUCを求め、平均血中濃度を返すメソッドを作成します。

    # 平均血中濃度
    def compute_cssave(self):
        # 点滴開始~点滴終了までの区間で定積分
        auc1 = integrate.quad(
            lambda t: (self.k0 / self.vd / self.k) * ((1 - np.exp(-self.k * self.T))
            / (1 - np.exp(-self.k * self.interval))
            * np.exp(-self.k * (t + self.interval - self.T)) + (1 - np.exp(-self.k * t))),
            0, self.T)

        # 点滴終了~次回時までの区間で定積分
        auc2 = integrate.quad(
            lambda t: (self.k0 / self.vd / self.k) * (1 - np.exp(-self.k * self.T))
            / (1 - np.exp(-self.k * self.interval)) * np.exp(-self.k * (t - self.T)),
            self.T, self.interval)

        return (auc1[0] + auc2[0]) / self.interval

コメントにあるとおり、点滴開始~点滴終了時(0~T)までと、点滴終了時~次回点滴時(T~τ)までを別々に積分し、それを足すことによって投与間隔区間のAUCを求めています。なお、積分しているのは以下の式です。

<点滴開始~点滴終了>
\small \displaystyle{ C_{ss}=\frac{k_0}{Vd \cdot k_e}\{ \frac{1-e^{-k_e T}}{1-e^{-k_eτ}} e^{-k_e(t+τ-T)}+(1-e^{-k_e t})\}}

<点滴終了~次回投与時まで>
\small \displaystyle{ C_{ss}=\frac{k_0}{Vd \cdot k_e} \cdot \frac{1-e^{-k_e T}}{1-e^{-k_eτ}} e^{-k_e(t-T)}}

また、上でも説明しましたが、AUCは\frac{C_0}{k_e}なので、ここはself.k0 * self.T / self.vd / self.kとしても同じように求めることができます。
returnのところで、求めたAUCの合計を投与間隔で割りCss_aveとして返しています。

Css_maxとCss_minについては、一番上で示した数式そのままです。

    # 最高血中濃度
    def calc_cssmax(self):
        return (self.k0 / self.vd / self.k) * (1 - np.exp(-self.k * self.T)) / \
               (1 - np.exp(-self.k * self.interval))
 
    # 最低血中濃度 
    def calc_cssmin(self):
        return self.calc_cssmax() * np.exp(-self.k * (self.interval - self.T))


今回は1日の投与量を固定(1200mg)してこれを1~3回に分割投与した際の3つのグラフを描画してみました(ページ上の図)。
平均血中濃度はどれも同じですが、投与回数が多いと最高ー最低血中濃度の差が小さくなっていることがわかります。

color:  blue
投与量: 1200, 投与回数: 1
定常状態到達時間:  72.0
平均血中濃度:  18.033688011112044
最高血中濃度:  33.417729436898156
最低血中濃度:  7.79496101671918
------------------------------
color:  red
投与量: 600, 投与回数: 2
定常状態到達時間:  72.0
平均血中濃度:  18.03368801111204
最高血中濃度:  23.981820514791366
最低血中濃度:  12.851539405530925
------------------------------
color:  green
投与量: 300, 投与回数: 4
定常状態到達時間:  72.0
平均血中濃度:  18.033688011112044
最高血中濃度:  19.90196072839641
最低血中濃度:  16.165415293827678
------------------------------


最後に、今回のコードを載せておきます。

import numpy as np
import matplotlib.pyplot as plt
import repeat_div
from scipy import integrate
 
class BloodConc2(repeat_div.BloodConc):
    def __init__(self, x0, vd, ke, T, interval):
        super().__init__(x0, vd, ke, T, interval)
        t_half = np.log(2) / ke
        # 投与間隔の倍数のうち、半減期の7倍以上の最小値を定常状態到達時間とする
        self.css_t = (t_half * 7 // interval + 1) * interval
 
    # 定常状態到達時間
    def get_csst(self):
        return self.css_t
 
    # 平均血中濃度 
    def compute_cssave(self):
        # 点滴開始~点滴終了までの区間で定積分
        auc1 = integrate.quad(
            lambda t: (self.k0 / self.vd / self.k) * ((1 - np.exp(-self.k * self.T))
            / (1 - np.exp(-self.k * self.interval))
            * np.exp(-self.k * (t + self.interval - self.T)) + (1 - np.exp(-self.k * t))),
            0, self.T)
 
        # 点滴終了~次回時までの区間で定積分
        auc2 = integrate.quad(
            lambda t: (self.k0 / self.vd / self.k) * (1 - np.exp(-self.k * self.T))
            / (1 - np.exp(-self.k * self.interval)) * np.exp(-self.k * (t - self.T)),
            self.T, self.interval)
 
        return (auc1[0] + auc2[0]) / self.interval
 
    # 最高血中濃度
    def calc_cssmax(self):
        return (self.k0 / self.vd / self.k) * (1 - np.exp(-self.k * self.T)) / \
               (1 - np.exp(-self.k * self.interval))
 
    # 最低血中濃度
    def calc_cssmin(self):
        return self.calc_cssmax() * np.exp(-self.k * (self.interval - self.T))
 
 
if __name__ == '__main__':
    vd = 40                    # 分布容積(L)
    t_half = 10                # 半減期(hr)
    k = np.log(2) / t_half     # 消失速度定数(hr^-1)
    T = 3                      # 点滴時間(hr)
    days = 4                   # 何日投与するか
 
    # 1回の投与量(mg), 投与間隔(hr)を変えて投与する. 1日の投与量は同じにする
    x0_1, interval_1 = 1200, 24      # 1200
    x0_2, interval_2 = 600, 12       # 600
    x0_3, interval_3 = 300, 6        # 300
 
    # BloodConc2(投与量, 分布容積, 消失速度定数, 点滴時間, 投与間隔)
    bc1 = BloodConc2(x0_1, vd, k, T, interval_1)    # blue
    bc2 = BloodConc2(x0_2, vd, k, T, interval_2)    # red
    bc3 = BloodConc2(x0_3, vd, k, T, interval_3)    # green
 
    # 描画期間の投与回数
    n1 = int(days * 24 / interval_1)
    n2 = int(days * 24 / interval_2)
    n3 = int(days * 24 / interval_3)
 
    conc1, conc2, conc3 = [], [], []
 
    # key: グラフの色, value: (インスタンス, 投与回数, 血中濃度)
    bc = {'blue': (bc1, n1, conc1), 'red': (bc2, n2, conc2), 'green': (bc3, n3, conc3)}
 
    period = 24 * days
    time = np.linspace(0, period, period / 0.1, endpoint=False)
 
    for color, value in bc.items():
        _bc, n, conc = value
        # 各クラス毎に1~n回目の投与までの血中濃度を加えていく
        for i in range(1, n + 1):
            conc.extend(_bc.get_section_conc(i))
        # 1回投与量, 1日投与回数
        dose = int(_bc.k0 * _bc.T)
        num = int(24 / _bc.interval)
        plt.plot(time, conc, c=color, alpha=0.8, label=str(dose) + 'mg×' + str(num) + '回/日')
        # 定常状態到達時間, 平均血中濃度, 最高血中濃度, 最低血中濃度の表示
        print('color: ', color)
        print('投与量: {}, 投与回数: {}'.format(dose, num))
        print('定常状態到達時間: ', _bc.get_csst())
        print('平均血中濃度: ', _bc.compute_cssave())
        print('最高血中濃度: ', _bc.calc_cssmax())
        print('最低血中濃度: ', _bc.calc_cssmin())
        print('-' * 30)
 

    # 軸ラベルを表示
    plt.ylabel('血中濃度[μg/mL]')
    plt.xlabel('時間[hr]')
    # 目盛りを24時間毎にする
    plt.xticks(np.arange(0, time[-1] + 24, 24))
    plt.xlim(0,)
    plt.ylim(0,)
    plt.grid(linestyle='dashed')
    plt.legend()
    plt.show()

参考文献

1. 徹底解説 薬物動態の数学ー微積分と対数、非線形ー
2. 4ステップ 臨床力UPエクササイズ4 TDM領域
3. https://www.pmda.go.jp/files/000206738.pdf

*1:薬物動態に関する速度 (吸収速度や代謝速度など )が投与量に比例すること(参考文献3 10.用語一覧)

*2:この積分区間の下限の0はn回目の投与開始時を表します