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

日本化学療法学会が公開しているPAT（Practical Antimicrobial TDM）は抗菌薬のTDMを実践するためのアプリケーションで、バンコマイシンのAUCガイドによる投与設計に対応しています。ソフトの詳細についてはリンクページを見ていただくことにして、今回はこのバンコマイシンTDMソフトウェア PATをお手本として、同じようなことができるプログラムをPythonで作ってみたいと思います。
最終的にはベイズ推定まで実装しますが、この記事はベイズ推定記事の前編にあたります。今回は血中濃度シミュレーションの実装を解説し、次回記事でベイズ推定による実装について書きます。

　
　

はじめに

本記事で解説・公開しているプログラムは、個人の学習および技術的な検証を目的としたものです。記事内では参考として日本化学療法学会のWebアプリケーション版PAT（Practical Antimicrobial TDM）の出力結果との比較を行っていますが、本プログラムはPATとは一切関係ありません。
　
　

環境

・Windows11
・Python3.12
・NumPy 2.4 / SciPy 1.17
　
　

プログラムの概要

今回のプログラムではPATで選択可能な母集団薬物動態モデルのうち、vancomycinAdult.PAT2024と同じ母集団パラメータを使用します（2026年6月更新）。
プログラムの機能は以下の通りです。

• 患者情報（性別、年齢、体重、Scr、筋肉量低下の有無等）と投与計画（投与量、点滴時間、投与間隔、繰り返し回数）を入力する
• 2-コンパートメント点滴静注モデルによる血中濃度-時間プロファイルを描画する
• 各投与回のトラフ値、投与初期のAUC（台形則）、定常状態のAUCを算出する
• 算出結果をグラフ上に表示する

これらの機能を持つ2-コンパートメントの点滴静注モデルのクラスを作成します。また、母集団PKパラメータを差し替えれば、他の薬物にも応用できる設計としています。
なお、血中濃度その他パラメータの計算には過去の記事で紹介した2コンパートメントモデルの解析解を使用しています。よろしければご参照ください。
www.yakupro.info

重ね合わせによる反復投与の血中濃度シミュレーション

体内での動態が線形（一次消失）を示すVCMのような薬物の場合、反復投与による任意の時間点における血中濃度は、単回投与時の血中濃度推移を投与間隔毎にずらして足し合わせることで求めることができます。

これは不均等投与を行った場合でも同じで、単回投与時の血中濃度を計算できればそれをベースに同様の手順で全体の推移を求めることができます。
今回のプログラムでは投与条件（投与歴）の設定はコード内でdict（辞書型）を使って以下のように記述することにしました。
下の例の場合、初回に1800mgを点滴時間2hrの12時間間隔で1回投与後、850mgを…（以下略）という具合になります。
keyの名前に制限はなく、重複さえなければ自由に設定可能です。たとえば'd3':[750, 1, 12, 3]のように新たにkeyとvalueを追加することで、いくつでも異なる投与条件を重ね合わせることができます。

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

不均等投与において、投与量や点滴時間が変更されると数式の指数項にかかる係数（A、B）の値も変わるので、keyが変わるタイミングでこれらの係数を再計算する処理を行っています。
なお、単回投与の濃度はどれくらいの経過時間まで計算するのか、については半減期とクラス変数FACTORとの積（self.half_life * self.FACTOR）によって決めるように実装しました。

def calc_single_dose_conc(self):
    # 1回投与量の濃度を計算
    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

半減期の5倍経過後には単回投与の濃度は約3%まで低下（(1/2)^5≒0.03）します。10倍経過後なら約0.1%（(1/2)^10≒0.001）で実質０とみなせるので、以前のプログラム（1998モデル使用）ではFACTOR=10と設定していましたが、今回のプログラムでは以下の理由からFACTOR=30へと変更しています。
　

半減期の算出方法について

あくまで私が確認した範囲においてではありますが、Web版PATの成人モデルで表示される半減期（）には以下のような違いがあるようです。

• vancomycinAdult.1998モデル：（β相（消失相）の半減期）
• vancomycinAdult.PAT2024モデル：（Vc（中央コンパートメント）からの消失速度定数による半減期）

前者の算出方法ではVpからVcへの動き（K21）の影響も含めた終末消失速度定数βを使っているため、Vpの大きなモデルではは間延びした長めの値になりますが、後者はCLとVcだけから決まるK10（=CL / Vc）を使ってが求められているため、単純な中央コンパートメントからの一次消失速度を反映した値（β相半減期よりも短い値）として算出されます。

例えば、患者情報（男性・60歳・60kg・Scr=0.8）のとき、β相の=21.5hrに対して、K10によるは6.44hrとなり、K10を基準としたFACTOR=10の地点（64.4hr）では実際の血中濃度（β相）がまだ約12%残存（(1/2)^3）していることになります。
今回のプログラムではvancomycinAdult.PAT2024モデルを使用するので、半減期もで算出することにしています。
そのため、以前のプログラムで使用していた1998モデルのFACTOR=10からFACTOR=30へと乗数を変更しています。

【追記】・・だったのですが、直近のアップデートでPAT2024モデルでも本来のβ相の半減期が表示されるようになったようです。しかし、算出される値がどういう意味を持つか考察した記録としてこの部分は残すことにしました。また、表示される値だけの違いなのでコードもそのままにしていますが、実際に表示をβ相の半減期に変えたい場合はコンストラクタ（__init__）内のself.half_lifeの式を変更すると確認することができます。*1

プログラムの全体的な動作としては、最初に投与の全区間の濃度を収めることができる長さのnumpy配列を確保しておき、単回投与の濃度を計算した後はスライスで投与間隔ごとに時間ずらしながら、投与回数分だけ足し合わせる処理を繰り返します。
コード全体はやや複雑になるため全容の掲載は割愛しますが、血中濃度シミュレーション部分については以上のようなものになります。

投与各回のトラフ値と定常状態トラフの算出

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

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）という意味になります。
　

定常状態のトラフ値

解析解は以下の式です。
　

　
リンク先の過去記事にある単回投与の点滴終了後の式より、等比級数の公式を使って導出できます。（※ ここでのA・Bは点滴速度k0（=D / T）を含んでいます。A・B・α・βの定義についてはページ上のリンク記事をご参照ください。）
　
最初に自分で実装したメソッドcalc_css_trough()は、最後の投与においてトラフ時刻にあたる濃度を投与間隔ごとに合計する、という方法です。（以前のプログラムではこのメソッドを使っていました）
解析解との比較では、少なくとも小数点以下２桁まで同じ、条件によっては3桁まで同じになることが確認できました（STEP=60、FACTOR=30、t_const=2のとき。重ね合わせ法の精度はこれらの設定値に依存します）。クラスには両メソッドとも残していますが、解析解の方はFACTOR等の設定に依存せず値の代入のみで求められるため、こちらを使うことにしました。

# 定常状態トラフ
解析解 : 14.274635141640047        # calc_css_trough_analytic()
重ね合わせ : 14.274337481807457    # calc_css_trough()

　
患者情報を以下の条件にして、上で例に出した投与条件（投与歴）の'd2'の投与回数を3→5回に変更してプログラムを実行すると、次のようなグラフが描画されます。

PATのトラフ値と比較したところ、条件によっては途中のトラフ値で0.01の差が出る場合がありましたが、概ね一致していることが確認できました。
誤差の発生する原因としては、時間をステップ幅に刻んだことによる誤差や、内部実装の差等が考えられます。
　
患者情報に関しては性別、年齢、体重、Scrのほか、（著明な）筋肉量低下があるか、Scr=0.6とする補正を行うか、Albの入力を行うかを選択できるようにしています。シスタチンCの入力については実装していません。

if __name__ == '__main__':
    """母集団薬物動態モデルはvancomycinAdult.PAT2024を使用"""
    # 患者情報
    sex = 'male'        # 'male' or 'female'
    age = 60
    weight = 60         # kg
    scr = 0.8           # mg/dL
    target_auc = 500    # 目標AUC(μg・h/mL) ベイズ推定（次回記事）で使用
    
    alb = None              # アルブミン[g/dL]（入力なしの場合はNone）
    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 = calc_ccr(sex, age, weight, scr)
    egfr = calc_egfr(sex, age, scr)
    
    CL_pop = calc_cl_pop(ccr, alb, muscle_wasting)
    Vc_pop = 32.1 * (weight / 60) ** 0.532
    Q_fixed = 8.33
    Vp_pop = 58.2 * (weight / 60) ** 0.348

2-コンパートメントの点滴静注モデルの構築に必要なパラメータについては次のモデル式を用いています。

出典：[PAT ver.4.1 操作マニュアル～薬物モデルの解説～]
バンコマイシンTDMソフトウェア PAT｜委員会報告・ガイドライン｜公益社団法人日本化学療法学会

バンコマイシンの成人（モデル名：vancomycinAdult.PAT2024）

PKパラメータ	母集団薬物動態モデル式
CL (L/h)	4.9×(Ccr/120)^0.959（Ccr>167.8の場合は167.8で固定）※
Vc (L)	32.1×(体重/60)^0.532
Q (L/h)	8.33【固定値】
Vp (L)	58.2×(体重/60)^0.348

※ CLはシスタチンC、Alb、筋肉量低下等のいずれの追加情報もない場合
患者集団の詳細についてはPAT ver.4.1操作マニュアルの薬物モデルの解説をご確認ください。個体間変動、残差変動については次回記事で解説しています。
　
　

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

PATでは定常状態（繰り返し投与時）の24時間AUC（とAUC/MIC）のほか、投与初期の0-24hrと24-48hrのAUCも表示されます。
定常状態の24時間AUCは「１日投与量 / CL」で計算できます。しかし、投与初期（0-24hrと24-48hr）の区間AUCは重ね合わせで得た離散的な濃度プロファイルから数値積分で求める必要があります。
離散的な値、というのは大げさにすると下のようなかんじです。

台形公式（台形則）

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

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

def calc_trapezoidal_auc_org(self, begin, end, N=1000):
    """台形公式の自作実装"""
    idx_begin = int(begin * self.STEP)
    idx_end = int(end * self.STEP)
    h = (idx_end - idx_begin) / N
    s = 0
    for k in range(1, N):
        s += self.conc[idx_begin + int(k * h)]
    s *= h
    s += h * (self.conc[idx_begin] + self.conc[idx_end]) * 0.5
    return s / self.STEP

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

数値積分と解析解との比較

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

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

　
SciPyにはquad()という連続関数を受け取り、精度が高くなるよう評価点を自動的に選んで数値積分を行う関数があります。これに単回投与の関数を与えて解析解と比較すると、以下の通り、quad()と解析解との差は約4.6×10⁻⁸で値はほぼ一致する結果となりました。

# AUC（単回投与[0-∞]）: SciPy quad(積分結果, 推定誤差)
(521.329959751323, 7.502977147016151e-06)
# AUC（単回投与[0-∞]）: (A+B) * T
521.3299597971854
# AUC_ss : (A+B) * T * 24 / tau 
1042.6599195943709
# AUC_ss : Dose / CL * 24 / tau
1042.6599195943716

　
次に、投与初期のAUCの数値積分の結果を表にして比べてみました。SciPy quad()は関数から積分を行い、Scipy trapezoid（台形公式ライブラリ）と自作関数はプログラムの重ね合わせで得たnumpy配列の離散値を積分しています。参考までに、PAT（ver. 4.1a）で表示される値も載せています。

患者情報は男性、年齢60、 体重60、 Scr 0.8とし、投与条件は1800mgを2hrで12時間間隔、１回投与としています。
（PAT以外は小数第3位を四捨五入、単位はμg*h/mL）

区間(hr)	SciPy quad	SciPy trapezoid	自作(N=1000)	PAT
0-24	314.77	314.65	314.71	314.43
24-48	111.21	111.16	111.24	111.59

［quad、scipy trapezoid、自作］の中で最も精度が高いSciPy quadが解析解にほぼ一致しているとすると、それぞれ誤差はあるものの、充分に実用範囲内だと言えそうです。自作実装もSciPy quadと遜色ない精度であることが確認できました。予想とは異なり、今回の条件では自作実装がScipy trapezoidより僅かにquadに近い結果となりました。trapezoidは配列の全点（60×24=1440）に対し、自作はN=1000点です。今回はたまたまサンプリング点が良かったためかこのような結果でしたが、条件によって精度は変わる可能性があります。
今回のプログラムでは既に計算済みの離散値に直接適用でき、実装も簡潔なSciPy trapezoid()を使うことにしました。私が試したところでは、trapezoidとPATの値とでは最大で2μg・h/mL程度のずれが生じる場合が確認できました。
　
次回は今回作成したプログラムを使って、ベイズ推定の実装をしていきたいと思います。
www.yakupro.info

参考

クリニカルファーマコメトリクス

今回作成したコードです。

from decimal import Decimal, ROUND_HALF_UP
from scipy.integrate import trapezoid, quad
import japanize_matplotlib
import matplotlib.pyplot as plt
import numpy as np

 
class TwoCompInfModel:
    STEP = 60     # 濃度プロファイルの時間単位を定義　1hr/60STEP=1min/step
    FACTOR = 30   # 半減期の何倍まで濃度計算するか. 10→30に変更する. 詳細は記事(1)を参照
 
    def __init__(self, *params, **kwargs):
        self.CL, self.Vc, self.Vp, Q_FIXED = params
        self.Vss = self.Vc + self.Vp
        self.K10 = self.CL / self.Vc
        self.K12 = Q_FIXED / self.Vc
        self.K21 = Q_FIXED / self.Vp
 
        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    # 以前の1998モデルではβ相の半減期を使用
        # PAT2024モデルではVcからの一次消失速度定数による半減期を使用
        self.half_life = 0.693 / self.K10
 
        tau_times_cycle = sum([x[2] * x[3] for x in kwargs.values()]) # dosage毎のtau*numの合計
        self.end_time = tau_times_cycle + int(self.half_life * 12) # β相の半減期にする場合は12→4に変更. 描画範囲調整のため.
        
        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.Vc * self.alpha * (self.beta - self.alpha))
        self.B = (K0 * (self.K21 - self.beta)) / (self.Vc * 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_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_ss_auc(self):
        # AUCss = 一日投与量 / CL
        return self.dose / self.CL * 24 / self.tau
 
    def calc_auc_mic_ratio(self, mic=1):
        """AUC/MIC比. 今回のプログラムでは未使用"""
        return self.calc_ss_auc() / mic
 
    def calc_trapezoidal_auc(self, begin, end):
        """Scipyライブラリによる実装. グラフ表示にはこちらを使用する"""
        idx_begin = int(begin * self.STEP)
        idx_end = int(end * self.STEP)
        auc = trapezoid(self.conc[idx_begin:idx_end])
        return auc / self.STEP
 
    def calc_trapezoidal_auc_org(self, begin, end, N=1000):
        """台形公式の自作実装"""
        idx_begin = int(begin * self.STEP)
        idx_end = int(end * self.STEP)
        h = (idx_end - idx_begin) / N
        s = 0
        for k in range(1, N):
            s += self.conc[idx_begin + int(k * h)]
        s *= h
        s += h * (self.conc[idx_begin] + self.conc[idx_end]) * 0.5
        return s / self.STEP
 
    def calc_css_trough_analytic(self):
        """定常状態トラフ：解析解. グラフ表示にはこちらを使用する"""
        # A, Bはcalc_exp_coeff()で計算済み（dose / T含む）
        tau = self.tau_list[-1]
        T = self.inf_t_list[-1]
        term_alpha = self.A * (1 - np.exp(-self.alpha * T)) * \
                    np.exp(-self.alpha * (tau - T)) / \
                    (1 - np.exp(-self.alpha * tau))
        term_beta = self.B * (1 - np.exp(-self.beta * T)) * \
                    np.exp(-self.beta * (tau - T)) / \
                    (1 - np.exp(-self.beta * tau))
        return term_alpha + term_beta
 
    def calc_css_trough(self, t_const=2):
        """定常状態トラフ：足し合わせ（トラフ時刻のcをtauごとに合計）"""
        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 f(self, t):
        """単回投与（点滴中・点滴終了後）の関数"""
        if t < self.T:
            return self.A * (1 - np.exp(-self.alpha * t)) + self.B * (1 - np.exp(-self.beta * t))
        else:
            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_quad_auc(self, begin, end):
        return quad(self.f, begin, end)
 
 
def round_v(f, digits='0.01'):
    return str(Decimal(str(f)).quantize(Decimal(digits), ROUND_HALF_UP))
 
 
def disp_trough_values(model, ax, y_offset=0, color='black'):
    MARGIN = 7
    for i, _ in enumerate(model.tau_list[:-1], 1):
        elapsed_time = sum(model.tau_list[:i])
        ax.text(
            elapsed_time - 2, model.get_cmax() + MARGIN + y_offset,
                round_v(model.get_trough(i)), color=color
        )
 
 
def disp_pk_result(model, ax, y_offset, label, color='black'):
    auc_early = (
            f'First   0-24hr AUC: '
            f'{round_v(model.calc_trapezoidal_auc(0, 24))} μg h/mL'
    )
    auc_late = (
            f'First 24-48hr AUC: '
            f'{round_v(model.calc_trapezoidal_auc(24, 48))} μg h/mL'
    )
    pk_line = (
            f'{label} PK値: '
            f'CL={round_v(model.CL)} L/h, '
            f'Vss={round_v(model.Vss)} L, '
            f'半減期={round_v(model.half_life)} hr, '
            f'AUC={round_v(model.calc_ss_auc())} μg h/mL, '
            f'トラフ値={round_v(model.calc_css_trough_analytic())} μg/mL'
    )
    ax.text(0, 1.21 + y_offset, auc_early, transform=ax.transAxes, color=color)
    ax.text(0, 1.16 + y_offset, auc_late, transform=ax.transAxes, color=color)
    ax.text(0.27, 1.21 + y_offset, pk_line, transform=ax.transAxes, color=color)
 
 
def plot_conc(model, ax, color='blue'):
    conc, time = model.get_conc()
    ax.plot(time, conc, alpha=0.8, linewidth=0.8, color=color)
 
 
def set_format(model, ax, ccr, egfr, muscle_wasting, dosage):
    TROUGH_MARGIN = 7
    # タイトル・モデル情報
    ax.text(0, 1.32, 'Result', fontsize=23, transform=ax.transAxes)
    ax.text(0.27, 1.32, 'Model: Vancomycin Adult (PAT2024パラメータ使用)', fontsize=12, transform=ax.transAxes, color='darkblue') # 0.27, 1.27
    ax.text(0.27, 1.26, # 0.66, 1.27
            f'CCr={round_v(ccr, '0.1')} mL/min,  '
            f'eGFR={round_v(egfr, '0.1')} mL/min/1.73m$^{2}$',
            transform=ax.transAxes)
    # 筋肉量低下時の腎機能について警告
    if muscle_wasting:
        ax.text(0.59, 1.26, '腎機能推定誤差が大きい可能性あり. 推定値に注意!', # 0.27, 1.33
                transform=ax.transAxes, fontsize=12, color='red')    # fontsize=11, 
    # トラフ値ラベル
    ax.text(0, model.get_cmax() + TROUGH_MARGIN, 'Trough')
    # 軸・グリッド設定
    ax.set_ylabel('Concentration(μg/mL)')
    ax.set_xlabel('Time after dose(hr)')
    ax.set_xticks(np.arange(0, model.time[-1] + 12, 12))
    ax.hlines([10, 20], -12, model.time[-1], color='red', linewidth=0.8)
    ax.tick_params(width=2, length=5)
    ax.set_xlim(-3, model.end_time)
    ax.set_ylim(-3, model.get_cmax() + 15)
    # 投与歴
    ax.text(0, -0.2, '<投与歴>', transform=ax.transAxes)
    for i, (dose, inf_t, tau, n) in enumerate(dosage.values()):
        line = f'{i+1}. Dose: {dose:4}mg,  Inf: {inf_t:3}h, Tau: {tau:3}h, N: {n:2}回'
        ax.text(0, -0.25 - i * 0.05, line, transform=ax.transAxes)
    
 
def calc_ccr(sex, age, weight, scr):
    """Cockcroft-Gault式によるCcr計算"""
    f = 0.85 if sex == 'female' else 1.0
    ccr = (140 - age) * weight / (72 * scr) * f
    return float(round_v(ccr, '0.1'))
 
 
def calc_egfr(sex, age, scr):
    """JSN eGFRcrによるeGFR計算"""
    f = 0.739 if sex == 'female' else 1.0
    return 194 * scr ** -1.094 * age ** -0.287 * f
 
 
def calc_cl_pop(ccr, alb=None, muscle_wasting=False):
    """
    vancomycinAdult.PAT2024のアルゴリズムに基づくCL_pop計算
    ※ シスタチンCによる計算は未実装. またAlb<2.5でもwarning表示なし
    詳細はPATver.4.1マニュアルを参照
    """
    if alb is not None and muscle_wasting:
        # ケース2-1): Alb入力あり + 筋肉量低下あり
        ccr_capped = min(ccr, 120)
        return 5.15 * (ccr_capped / 120) ** 0.842 * (alb / 4) ** 0.345

    if alb is not None and not muscle_wasting:
        # ケース2-2): Alb入力あり + 筋肉量低下なし
        ccr_capped = min(ccr, 170)
        return 5.54 * (ccr_capped / 120) ** 0.903 * (alb / 4) ** 0.361

    if muscle_wasting:
        # ケース3: Alb/シスタチンCなし + 筋肉量低下あり
        ccr_capped = min(ccr, 145)
        return 4.22 * (ccr_capped / 120) ** 0.859

    # デフォルト: 基本形
    ccr_capped = min(ccr, 167.8)
    return 4.9 * (ccr_capped / 120) ** 0.959
 
 
if __name__ == '__main__':
    """母集団薬物動態モデルはvancomycinAdult.PAT2024を使用"""
    # 患者情報
    sex = 'male'        # 'male' or 'female'
    age = 60
    weight = 60         # kg
    scr = 0.8           # mg/dL
    target_auc = 500    # 目標AUC(μg・h/mL) ベイズ推定（次回記事）で使用
    
    alb = None              # アルブミン[g/dL]（入力なしの場合はNone）
    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 = calc_ccr(sex, age, weight, scr)
    egfr = calc_egfr(sex, age, scr)
    
    CL_pop = calc_cl_pop(ccr, alb, muscle_wasting)
    Vc_pop = 32.1 * (weight / 60) ** 0.532
    Q_fixed = 8.33
    Vp_pop = 58.2 * (weight / 60) ** 0.348
 
    # 投与条件の設定. 不均等投与に応じてdosageを作成
    # [投与量(mg), 点滴時間(hr), 投与間隔(hr), 投与回数]
    dosage = {
        'd1': [1800, 2, 12, 1],
        #'d2': [850, 1, 12, 3]
    }
    
    params = [CL_pop, Vc_pop, Vp_pop, Q_fixed]
    pop_model = TwoCompInfModel(*params, **dosage)
    print('定常状態トラフ')
    print(f'解析解 : {pop_model.calc_css_trough_analytic()}')
    print(f'重ね合わせ : {pop_model.calc_css_trough()}')
    print()
    
    A, B = pop_model.A, pop_model.B
    T = pop_model.T
    tau = pop_model.tau
    print(f'AUC（単回投与[0-∞]）quad : {pop_model.calc_quad_auc(0, float('inf'))}')
    print(f'AUC（単回投与[0-∞]）(A+B)*T: {(A + B) * T}')
    print(f'AUC_ss (A+B)*T*24/tau : {(A + B) * T * 24 / tau}')
    print(f'AUC_ss（Dose/CL*24/tau: {pop_model.calc_ss_auc()}')
    print()
    
    print('*--- 投与初期の区間AUC比較（単回投与時） ---*')
    print(f'quad AUC_0-24  : {pop_model.calc_quad_auc(0, 24)}')
    print(f'quad AUC_24-48 : {pop_model.calc_quad_auc(24, 48)}')
    print(f'trapezoid AUC_0-24  : {pop_model.calc_trapezoidal_auc(0, 24)}')
    print(f'trapezoid AUC_24-48 : {pop_model.calc_trapezoidal_auc(24, 48)}')
    print(f'自作(N=1000) AUC_0-24  : {pop_model.calc_trapezoidal_auc_org(0, 24)}')
    print(f'自作(N=1000) AUC_24-48 : {pop_model.calc_trapezoidal_auc_org(24, 48)}')
    
    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_values(pop_model, ax)
    disp_pk_result(pop_model, ax, 0, '〇母集団平均')

    plt.show()

*1:グラフ描画を見やすくするため、その場合はFACTORとself.end_timeの値も小さくしてください。