Pythonでバンコマイシンの血中濃度プロファイルをシミュレーションする

日本化学療法学会がホームページで公開しているバンコマイシンTDMソフトウェアで
PAT（Practical AUC-guided TDM for vancomycin）というWebアプリケーションがあります。名前が示すように、AUCガイドによる投与設計を実践するためのTDMソフトです。ソフトの詳細についてはリンクページを見ていただくことにして、今回はこのPATをお手本に、同じようなことができるプログラムをPythonで作ってみたいと思います。最終的にはベイズ推定まで実装しますが、最初に言っておくと、今回はベイズの話は出てきません。そこまでに必要なプログラム作成のお話となります。

　
　

環境

・Windows10
・Python3.12
・NumPy 1.26
・SciPy 1.13
Pythonはdictが順序を保証するようになった3.7以降であれば特に問題ないと思います。他にも必要なライブラリはありますが、そんなに珍しいものはないと思います。
　
　

プログラムの概要

• 年齢、体重、Scr、性別、著明な筋肉量低下があるか、Scr=0.6とする補正を行うか等の患者情報・設定を入力し、それに基づいてバンコマイシン（以下VCM）のクリアランス（CL）を算出する。
• 成人患者のVCM母集団PKパラメータを用いて、血中濃度-時間プロファイルを描画する。ここで使用する母集団PKパラメータはPATと同じものです。
• トラフ値や投与初期におけるAUCを算出する。
• グラフ描画の他、Figureには算出した結果を表示できるようにする。

これらの機能を持つ2-コンパートメントの点滴静注モデルのクラスを作成します*1。別の薬物の母集団PKパラメータがあれば、同様のことができるものを作りたいと思います。
なお、今回プログラム中で使用する計算式はこちらで書いた記事のものを使っています。よろしければ参考にしてください。
www.yakupro.info

重ね合わせにより反復投与の血中濃度を求める

体内での分布が線形であり、一次消失を示すVCMのような薬物の反復投与の経過時間における血中濃度は、単回投与時の血中濃度を投与間隔毎にずらして重ね合わせた合計値として求めることができます。

これは不均等投与を行った場合でも同じで、単回投与時の血中濃度を計算できれば同様の方法で求めることができます。
投与条件（投与歴）の設定はPATと同じような形式で、コード内にdictを使って以下のように記述することにしました。下の例だと、初回に1500mgを点滴時間1.5の8時間間隔で1回投与後、1000mgを…（以下略）という具合です。
keyの名前に制限はありません、重複さえなければOKです。また、...,'d3':[750, 1, 12, 3]のように新たにkeyとvalueを追加することで、いくつでも異なる投与条件を重ねることができます。

# 投与条件の設定. 不均等投与にも対応
# [投与量(mg), 点滴時間(hr), 投与間隔(hr), 投与回数]
dosage = {
        'd1': [1500, 1.5, 8, 1],
        'd2': [1000, 1, 12, 4]
}

不均等投与の場合、投与量や点滴時間が変わると指数項の係数A、Bの値も変わるので、keyが変わるタイミングで再計算しています。
なお、単回投与の濃度の経過時間はどれくらいまで計算するのか、については（消失相の）半減期とクラス変数FACTORとの積によって決めるようにしています。通常、半減期の4-5倍の時間も経てば定常状態に達すると言われますが（濃度は(1/2)^5≒0.03となる）、正確さを考慮すると5倍では全く短かかったので、10倍とすることにしました。
最初に投与の全区間の濃度を収めることができる長さのnumpy配列を確保しておき、投与条件の1回投与の濃度を計算した後は、スライスで投与間隔時間ずらしながら、投与回数分合計を繰り返します。
コードそのものはやや複雑なので、ここでは一つ一つ取り上げませんが、計算方法の概要は以上のようなものになります。

投与回数毎のトラフ値と定常状態のトラフ濃度を求める

各投与回数ごとのトラフ値を求めるメソッドは次のように作成しました。

def get_trough(self, num_doses):
    # 指定された投与回数時のトラフ値を返す. 範囲を超える場合は最後のトラフを返す
    elapsed_time = np.sum(self.tau_list[:num_doses])
    return self.conc[int(elapsed_time * self.STEP)]

全区間の血中濃度はもう算出されているので、投与回数num_dosesを指定したときに、トラフ値となる経過時間がわかればOKです。このために各投与間隔時間のリストをインスタンス生成時に作ることにしました。例えば、self.tau_listが[8, 12, 12, 12, ...]だったとすると、投与回数が4回目のトラフ時間は投与開始44時間後です。
クラス変数のSTEPは、コード中のコメントにもあるように濃度プロファイルの時間単位（刻み幅）を定義します。今回設定した60なら、1stepで1分（1hr / 60step）という意味になります。
Figureへの表示のコードについては、枝葉の部分になるので省略します。

定常状態のトラフ値については、次のようにして求めることにしました。
考え方は重ね合わせと同じで、最後の投与のトラフ値の時間となる部分の濃度を合計する、というものです。メソッドはcalc_css_trough()です。
　
患者情報を以下の条件にして、上で例に出した投与条件（投与歴）の'd2'の投与回数を4→6にした場合、次のようなグラフが描画されました。

5回目投与のトラフ値がPATで出される値と0.01だけ違いますが、プログラムは正しく計算できているようです。6回目投与の結果は一致しているので、原因としては、どこかの段階で丸め誤差が発生しているのか、時間ステップ（刻み幅）の影響等が考えられます。重ね合わせ期間をFACTOR=13にして延ばしたり、クラス変数のSTEPを80以上に増やすとこの値も一致するようになりますが、ここはプログラムが合っているという確認ができたので、特に変更は加えずにこのままいくことにしました。

 
def calc_vancomycin_cl(sex, age, weight, scr, muscle_wasting):
    f = 0.85 if sex == 'female' else 1.0
    ccr = ((140 - age) * weight) / (72 * scr) * f   # Cockcroft-Gault
    ccr = float(round_v(ccr, '0.1'))  # PATはここで丸めてる
    if muscle_wasting and ccr > 85:
        print("Renal function estimation may be inaccurate. ")
        cl = 3.51
    else:
        cl = 0.0478 * ccr
    return ccr, cl
 
if __name__ == '__main__':
    # 患者情報の入力
    age = 60
    weight = 60
    scr = 0.8
    sex = 'male'      # 'male' or 'female'
    muscle_wasting = False   # 著明な筋肉量低下があるか
    scr_round_up = False     # Scr<0.6の場合にScr=0.6とする補正を行うか
    scr = 0.6 if scr_round_up and scr < 0.6 else scr

    ccr, CL_pop = calc_vancomycin_cl(sex, age, weight, scr, muscle_wasting)
    egfr = calc_egfr(sex, age, scr)

患者の腎機能がCCr≤85mL/minの場合はVCMのCLは次の式で求めます。ただし、CCr＞85で、かつ、著明な筋肉量低下がある場合はCL=3.51とします。
　
これが母集団平均のCLとなります。ついでに、他のパラメータについては以下のとおりです。患者情報に関しては年齢、体重、Scr、性別のほか、著明な筋肉量低下があるか、Scr=0.6とする補正を行うかを入力できるようにしています。

	母集団平均（θ）
CL(L/h)	0.0478 × CCr (CCr ≤ 85mL/min) 3.51 ( CCr > 85mL/min)
K12(/hr)	0.525
K21(/hr)	0.213
Vss(L)	60.7

ちなみに、もしWeb版PATをお使いになる場合は、リンクページの「注意事項」や「利用規約」、「マニュアル」等を熟読の上、サーバーに過剰な負荷をかけないようご注意ください。
　
　

台形公式による数値積分で投与初期のAUCを求める

PATでは定常状態（繰り返し投与時）の24時間AUC（とAUC/MIC）のほか、0-24hrと24-48hrのAUCも表示されます。定常状態のAUCは積分すべき関数から解析的に（数式を使って*2）求めることができますが、重ね合わさった濃度プロファイルは直接関数から求めたのではなく、時間ごとの離散的な値を算出しただけなので、区間[0-24]や[24-48]の定積分は数値計算によって求めることになります。
離散的な値、というのは大げさにすると下のようなかんじです。

台形公式（台形則）

一言でいうと、定積分

の値を小さく分割した台形の面積の和で近似的に求める方法です。
積分区間を分割して、分割された小区分の和を求める方法といえば、シンプルなものでは区分求積法（長方形近似）がありますが、こちらは被積分関数を定数（１点）で近似しているのに対して、台形公式では一次関数（直線）で近似しています。ちなみに、二次関数で近似する方法にはシンプソンの公式というものがあります。
概念やわかりやすい図については
台形公式 - Wikipedia
にありました。
ここでは天下りではありますが、区間幅が等間隔の場合の公式を引用します。

これをそのまま実装したものがこちら。

def calc_trapezoidal_auc(self, begin, end, N=1000):
    begin *= self.STEP
    end *= self.STEP
    h = (end - begin) / N
    s = 0
    for k in range(1, N):
        s += self.conc[begin + int(k * h)]
    s *= h
    s += h * (self.conc[begin] + self.conc[end]) * 0.5    # Nが充分大きければ, 消しても同じ
    return s / self.STEP

参考ページにもあるように、分割数Nが充分大きい場合は、数式の第1項は全体の和の中での寄与が相対的に小さくなり、計算結果にはほとんど影響しなくなります。コード中でいうと、return文の一つ上のステートメントです。
ライブラリでは、scipy.integrateモジュールにtrapezoid()関数があります。

解析解との比較

AUCの値はTDMソフトの目的から考えて、高い精度が要求されるものではありませんが、どの程度であるのか少々気にはなるところです。
単回投与時は被積分関数がはっきり与えられており、この関数の区間[0, ∞]におけるAUCは解析的に求めることができます。
*3
または

これらはもちろん、すべて同じ値になります。
またSciPyにはquad()という数値積分を計算する関数があり、単回投与の関数を与えて解析的な値と比較すると、以下のとおり、quadはかなり精度が高いことがわかります。

# quad 0-∞ (積分結果, 推定誤差)
(376.71972554717047, 3.6446181184146553e-07) 
# (A+B)*T
376.71972554712266 
# (A+B)*T に１日投与回数を乗じたもの. これがPATでグラフ上部に表示される値
1130.159176641368

他の関数の結果とも比べて表にしてみました。SciPyのquad()のみ関数から積分を行い、trapezoidと自作関数は今回のプログラムで得たnumpy配列の離散データを積分しています。参考までに、PATで表示される値も載せています。
投与条件は年齢=60、 体重=60、 Scr=0.8、男性で1500mgを1.5hrで8時間間隔、１回投与です。（PAT以外は小数第3位を四捨五入、単位はμg*h/mL）

区間(hr)	SciPy quad	SciPy trapezoid	自作(N=1000)	PAT
0-24	289.25	289.16	289.19	289.02
24-48	62.99	62.96	63.01	63.36
0-∞	376.72	-	-	-

SciPy quadが（ほぼ）真の値だとすると、他3つはそれなりに誤差はあるものの、充分実用範囲内だと言えると思います。自作も全然悪くないですね。と確認できたところで、今回のプログラムでは大人しくSciPy trapezoid()を実装に使うことにしました。私が試したところでは、trapezoidとPATの値とでは最大で2程度のずれが生じる場合もありました。
　
　
今回はこれで以上です。下に今回作成したコードを載せておきます。
次回は今回作成したプログラムを使って、ベイズ推定の実装をしていきたいと思います。

import matplotlib.pyplot as plt
import numpy as np
from decimal import Decimal, ROUND_HALF_UP
from scipy.integrate import trapezoid
 
 
class TwoCompInfModel:
    STEP = 60    # 濃度プロファイルの時間単位を定義　1hr/60STEP=1min/step
    FACTOR = 10  # 半減期の何倍まで計算するか
 
    def __init__(self, *params, **kwargs):
        self.CL, self.Vss, self.K21, self.K12 = params
        self.Vd1 = self.Vss / (1 + self.K12 / self.K21)
        self.K10 = self.CL / self.Vd1
 
        x = self.K12 + self.K21 + self.K10
        y = self.K21 * self.K10
        self.alpha = 0.5 * (x + np.sqrt(x ** 2 - 4 * y))
        self.beta = x - self.alpha
        self.half_life = 0.693 / self.beta  # np.log(2) / self.beta
 
        tau_times_cycle = sum([x[2] * x[3] for x in kwargs.values()])  # dosage毎のtau*numの合計
        self.period = int(tau_times_cycle + self.half_life * self.FACTOR)
        self.conc = np.zeros(self.period * self.STEP)
        self.time = np.linspace(0, self.period, self.period * self.STEP, endpoint=False)
        self.total_n = 0  # 投与回数を保持
        self.tau_list = []
        self.inf_t_list = []
        self.create_dosing_conditions_list(kwargs)
        self.calc_conc(kwargs)
 
    def create_dosing_conditions_list(self, dosage):
        for d in dosage.values():
            _, inf_t, tau, cycle = d
            self.tau_list += [tau for _ in range(cycle)]
            self.inf_t_list += [inf_t for _ in range(cycle)]
        self.tau_list = np.array(self.tau_list)
        self.inf_t_list = np.array(self.inf_t_list)
 
    def calc_conc(self, dosage):
        for d in dosage.values():
            self.dose, self.T, self.tau, self.cycle = d
            self.calc_exp_coeff()
            self.calc_overlapping_conc()
 
    def calc_exp_coeff(self):
        K0 = self.dose / self.T
        self.A = (K0 * (self.K21 - self.alpha)) / (self.Vd1 * self.alpha * (self.beta - self.alpha))
        self.B = (K0 * (self.K21 - self.beta)) / (self.Vd1 * self.beta * (self.alpha - self.beta))
 
    def calc_during_inf(self):
        t = np.linspace(0, self.T, num=int(self.T * self.STEP), endpoint=False)
        return self.A * (1 - np.exp(-self.alpha * t)) + self.B * (1 - np.exp(-self.beta * t))
 
    def calc_post_inf(self, begin, end):
        t = np.linspace(begin, end, num=int((end - begin) * self.STEP), endpoint=False)
        return self.A * (np.exp(-self.alpha * (t - self.T)) - np.exp(-self.alpha * t)) \
               + self.B * (np.exp(-self.beta * (t - self.T)) - np.exp(-self.beta * t))
 
    def calc_single_dose_conc(self):
        c1 = self.calc_during_inf()
        c2 = self.calc_post_inf(self.T, self.half_life * self.FACTOR)
        c = np.concatenate([c1, c2], axis=0)
        return c
 
    def calc_overlapping_conc(self):
        c = self.calc_single_dose_conc()
        length = len(c)
        loop_cnt = 0
        for i, _ in enumerate(self.tau_list, self.total_n):
            elapsed_time = np.sum(self.tau_list[:i])
            start_idx = int(elapsed_time * self.STEP)
            end_idx = start_idx + length
            self.conc[start_idx:end_idx] += c
            loop_cnt += 1
            if loop_cnt == self.cycle:
                self.total_n += self.cycle
                break
  
    def get_conc(self):
        return self.conc, self.time
 
    def get_cl(self):
        return self.CL
        # CL = Dose / AUC
        # return self.dose / ((self.A + self.B) * self.T)
 
    def get_cmax(self):
        # トラフ値textの座標取得のため
        return np.max(self.conc)
 
    def get_trough(self, num_doses):
        # 指定された投与回数のトラフ値を返す. 範囲を超える場合は最後のトラフを返す
        elapsed_time = np.sum(self.tau_list[:num_doses])
        return self.conc[int(elapsed_time * self.STEP)]
 
    def calc_auc(self):
        # AUC = Dose / CL
        return self.dose / self.CL * 24 / self.tau
 
    def calc_auc_mic_ratio(self, mic=1):
        return self.calc_auc() / mic
 
    def calc_trapezoidal_auc(self, begin, end):
        begin *= self.STEP
        end *= self.STEP
        auc = trapezoid(self.conc[begin:end])
        return auc / self.STEP
 
    def calc_trapezoidal_auc_org(self, begin, end, N=1000):
        begin *= self.STEP
        end *= self.STEP
        h = (end - begin) / N
        s = 0
        for k in range(1, N):
            s += self.conc[begin + int(k * h)]
        s *= h
        s += h * (self.conc[begin] + self.conc[end]) * 0.5
        return s / self.STEP
 
    def calc_css_trough(self, t_const=1.7):
        css_range = self.half_life * self.FACTOR * t_const
        c = self.calc_post_inf(self.T, css_range)
        n = css_range // self.tau
        tr_idx = (np.arange(1, n) * self.tau - self.T) * self.STEP
        return np.sum(c[tr_idx.astype(np.int32)])
 
 
def round_v(f, digits='0.01'):
    return str(Decimal(str(f)).quantize(Decimal(digits), ROUND_HALF_UP))
 
 
def disp_trough(TwoCompInfModel, ax, add_pos=0, color='black'):
    margin = 7
    for i, t in enumerate(TwoCompInfModel.tau_list[:-1], 1):
        etime = sum(TwoCompInfModel.tau_list[:i])
        ax.text(etime - 2, TwoCompInfModel.get_cmax() + margin + add_pos,
                round_v(TwoCompInfModel.get_trough(i)), color=color)
 
 
def disp_text(TwoCompInfModel, ax, add_pos, situation, color='black'):
    ax.text(0, 1.21 + add_pos, 'First   0-24hr AUC: ' +
            round_v(TwoCompInfModel.calc_trapezoidal_auc(0, 24)) + ' μg h/mL',
            transform=ax.transAxes, color=color)
    ax.text(0, 1.16 + add_pos, 'First 24-48hr AUC: ' +
            round_v(TwoCompInfModel.calc_trapezoidal_auc(24, 48)) + ' μg h/mL',
            transform=ax.transAxes, color=color)
    ax.text(0.27, 1.21 + add_pos, situation + ' PK値: ' +
            'CL=' + round_v(TwoCompInfModel.get_cl()) + ' L/h, ' +
            'Vd=' + round_v(TwoCompInfModel.Vss) + ' L, ' +
            '半減期=' + round_v(TwoCompInfModel.half_life) + ' hr, ' +
            'AUC=' + round_v(TwoCompInfModel.calc_auc()) + ' μg h/mL, ' +
            'トラフ値=' + round_v(TwoCompInfModel.calc_css_trough()) + ' μg/mL',
            transform=ax.transAxes, color=color)
 
 
def plot_conc(TwoCompInfModel, ax, color='blue'):
    conc, time = TwoCompInfModel.get_conc()
    ax.plot(time, conc, alpha=0.8, linewidth=0.8, color=color)
 
 
def set_format(TwoCompInfModel, ax, ccr, egfr, warning, dosage):
    margin = 7
    ax.text(0, 1.29, 'Result', fontsize=21, transform=ax.transAxes)
    ax.text(0.27, 1.27, 'Model: Vancomycin (Adult) ', fontsize=12, transform=ax.transAxes)
    ax.text(0.50, 1.27,
            'CCr=' + round_v(ccr, '0.1') + ' mL/min,  ' +
            'eGFR=' + round_v(egfr, '0.1') + ' mL/min/1.73m$^{2}$',
            transform=ax.transAxes)
    ax.text(0, TwoCompInfModel.get_cmax() + margin, 'Trough')
    ax.set_ylabel('Concentration(μg/mL)')
    ax.set_xlabel('Time after dose(hr)')
    ax.set_xticks(np.arange(0, TwoCompInfModel.time[-1] + 12, 12))
    ax.hlines(10, -12, TwoCompInfModel.time[-1], color='red', linewidth=0.8)
    ax.hlines(20, -12, TwoCompInfModel.time[-1], color='red', linewidth=0.8)
    ax.tick_params(width=2, length=5)
    ax.set_xlim(-3, TwoCompInfModel.period - TwoCompInfModel.half_life * 7)
    ax.set_ylim(-3, TwoCompInfModel.get_cmax() + 15)
    ax.text(0, -0.2, '<投与歴>', transform=ax.transAxes)
    for i, d in enumerate(dosage.values()):
        dose, inf_t, tau, n = d
        dosage = f'{i+1}. Dose: {dose:4},  Inf_time: {inf_t:3}, Tau: {tau:3}, N: {n:2}'
        ax.text(0, -0.25 - i * 0.05, dosage, transform=ax.transAxes)
    if warning:
        ax.text(0.50, 1.33, 'Warning: Renal function estimation may be inaccurate.',
                transform=ax.transAxes, fontsize=11, color='red')
 
 
def calc_vancomycin_cl(sex, age, weight, scr, muscle_wasting):
    f = 0.85 if sex == 'female' else 1.0
    ccr = ((140 - age) * weight) / (72 * scr) * f   # Cockcroft-Gault
    ccr = float(round_v(ccr, '0.1'))  # PATはここで丸めてる
    if muscle_wasting and ccr > 85:
        print("Renal function estimation may be inaccurate. ")
        cl = 3.51
    else:
        cl = 0.0478 * ccr
    return ccr, cl
 
 
def calc_egfr(sex, age, scr):
    f = 0.739 if sex == 'female' else 1.0
    return 194 * scr ** -1.094 * age ** -0.287 * f
 
 
if __name__ == '__main__':
    # 患者情報の入力
    age = 60
    weight = 60
    scr = 0.8
    sex = 'male'    # 'male' or 'female'
    muscle_wasting = False # 著明な筋肉量低下があるか
    scr_round_up = False   # Scr<0.6の場合にScr=0.6とする補正を行うか
    scr = 0.6 if scr_round_up and scr < 0.6 else scr
 
    ccr, CL_pop = calc_vancomycin_cl(sex, age, weight, scr, muscle_wasting)
    egfr = calc_egfr(sex, age, scr)
 
    # PATで使用されているものと同じ母集団薬物動態パラメータを使用
    Vss_pop = 60.7       # Vss = Vd1 + Vd2 = (1 + k12 / k21) * Vd1
    K12_fixed = 0.525    # 中央→抹消コンパートメントへの消失速度定数(/hr)
    K21_pop = 0.213      # 抹消→中央コンパートメントへの消失速度定数(/hr)
 
    params = [CL_pop, Vss_pop, K21_pop, K12_fixed]
 
    # 投与条件の設定. 不均等投与にも対応
    # [投与量(mg), 点滴時間(hr), 投与間隔(hr), 投与回数]
    # Set dosing conditions. 'Key' names cannot be the same.
    # Format: [dose(mg), infusion time(hr), dosing interval(hr), num of doses]
    dosage = {
        'd1': [1500, 1.5, 8, 1],
        'd2': [1000, 1, 12, 4]
    }
 
    pop_model = TwoCompInfModel(*params, **dosage)
 
    fig = plt.figure(figsize=(12.5, 7))
    fig.subplots_adjust(right=0.95, left=0.08, top=0.75, bottom=0.22)
    ax = fig.add_subplot(111)
 
    set_format(pop_model, ax, ccr, egfr, muscle_wasting, dosage)
    plot_conc(pop_model, ax)
    disp_trough(pop_model, ax)
    disp_text(pop_model, ax, 0, '〇母集団平均')
 
    plt.show()
 

*1:表示部分のコードはクラスとは別です

*2:後述するDose/CLならそもそも与えられた値のみで計算できます

*3:A/α＋B/βとなっているものもありますが、その場合はA, Bの定義が異なります。A, Bの定義についてはページ上部のリンク記事をご参照ください