バンコマイシンTDM（２）－ベイズ推定プログラムを実装するー

前回記事では患者情報と母集団PKパラメータを用いてバンコマイシン（VCM）の血中濃度推移をシミュレーションするプログラムを作りました。今回はそのプログラムを拡張し、ベイズ推定により対象患者個人のPKパラメータを推定して推奨投与量を算出するプログラムを作っていきたいと思います。

2026年6月　修正・追記

この記事の旧バージョンには、以下の2点の誤りがありました。

１．残差変動の誤差モデルの誤り
血中濃度の残差変動（個体内変動）に指数誤差モデルを使用していましたが、通常この部分には指数誤差モデルは用いられません。PKパラメータの個体間変動には指数誤差モデル、残差変動には対象となる患者やデータに応じて付加誤差モデル、比例誤差モデル、または両者を組み合わせた混合誤差モデルが通常用いられます。（詳細は本文中で改めて説明します）。

２．理論的根拠のない正則化パラメータλの導入
旧版では目的関数の事前分布項にハイパーパラメータλ=1.5を乗じていましたが、本来のMAP推定の定式化にはこのような任意の係数は登場しません。[１]の誤りのままお手本としたPATの結果に近づけようとした結果、λによる後付け調整を行っていました。本来は理論式の通りのもの（係数なし）を最小化します。
当時は自分のこのあたりの理解が及んでおらず、このような間違いをしていました。修正後もやや曖昧な部分も残りますが、今回は以下を修正・更新して再度試みました。

• 残差変動部分の誤差モデルと目的関数の書き直し
• 理論的根拠のない正則化パラメータλの削除
• コード修正および図（グラフ）の差し替え

理論的に正しいと判断した部分はそのまま残しています。
　
　

前置き（と、おことわり）

本記事で解説・公開しているプログラムは、個人の学習および技術的な検証を目的としたものです。記事内では参考として日本化学療法学会のWebアプリケーション版PAT（Practical Antimicrobial TDM）の出力結果との比較を行っていますが、本プログラムはPATとは一切関係ありません。
　
記事の趣旨としては、Pythonプログラミングによるベイズ推定の実装に主眼を置いています（理論的な部分はやや粗があるかもしれません）。臨床の現場ではTDMソフトを使っているけれど、ベイズ推定についてはよくわかっていない、とかこれから勉強してみたいと思っている薬剤師や学生の方にとって、少しでも理解の助けになるところがあれば幸いです。

ちなみに、Pythonのバージョン等の環境、VCMの濃度シミュレーションのプログラム部分については前回と同様です。
www.yakupro.info

TDMにおけるベイズ推定

では本題に入っていきます。

目的

ある患者の特定の投与条件下において、実際に得られたデータ（実測血中濃度）に基づいて薬物動態パラメータ（CLやVd）を個別に推定し、その推定値を用いて次回以降の投与量・投与間隔を最適化すること
　

誤差モデルの設定

いくつか前提があるのでそれらについて確認していきます。

【１】 推定したい薬物動態パラメータはある特定の確率分布を持つと仮定します。このPKパラメータが持つ分布やその代表値（平均や分散）は、対象の患者が属する（と想定された）母集団から解析されたものを使用します。

今回のプログラムではPATで選択可能な母集団薬物動態モデルのうち、
vancomycinAdult.PAT2024と同じ母集団パラメータを使用します。

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

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

母集団薬物動態モデル式 （【】内は個体間変動 ）
CL (L/h)=4.9×(Ccr/120)^0.959 【36.9%】
（Ccr>167.8の場合は167.8で固定）
Vc (L)=32.1×(体重/60)^0.532 【33.3%】
Q (L/h)=8.33 【固定値】
Vp (L)=58.2×(体重/60)^0.348 【72.8%】
残差変動 18.1%

患者集団の詳細についてはPAT ver.4.1操作マニュアルの薬物モデルの解説をご確認ください。

例として、クリアランス（CL）はモデル式の基本形が次のようになっています。
（※ シスタチンC、Alb、筋肉量低下等のいずれの追加情報もない場合）

（Ccr>167.8の場合は167.8で固定）【36.9%】
Ccr：mL/min

ここで、式の中に組み込まれている Ccr（クレアチニンクリアランス）は、個体差を説明する因子で共変量と呼ばれます（Vc、Vpでは体重が共変量として使用されています）。
あらかじめ患者個々の腎機能（Ccr）に応じてCLの母集団平均​（その集団で最も代表的なCL値。Typical Value）を調整しておくことで、説明しきれないばらつきを小さくし、予測の精度を高めています。
しかし、同じCcrの患者でもCLは人によって異なっており、未知のばらつきが残ります。このばらつきのことを個体間変動 と呼びます。
　
個体間変動は患者ごとの母集団平均からのずれを表し、平均=0、分散=の正規分布に従うと仮定します。このとき正規分布の性質から、の約68.3%が1標準偏差内（）に入り、約95.4%がの範囲に、約99.7%がの範囲に入ります。

また、CLやVd（分布容積）といった各個体のPKパラメータは正の値しか取らず、右に裾を引く形で幅広く分布（対数正規分布）していると考えられることから、指数誤差モデルが適用されることが多くなります。
指数誤差モデルでの個体のは、を母集団平均とすると以下のように表されます。

　
下にCcr=83.3mL/minのときのCLの対数正規分布のプロットと、指数誤差モデル（青線）と比例誤差モデル（緑線）におけるCL値の比較を示します。

この場合の±2ωのCLを計算してみると
（後述）より、
・-2：3.45×exp(-2 × 0.357) ≈ 1.69
・+2：3.45×exp(+2 × 0.357) ≈ 7.05
となるため、Ccr=83.3の患者を100人集めたとき、約95人のCLは1.7～7.0L/hの範囲に収まると解釈できます。
　
ベイズ推定の計算を行うには分布の分散であるが必要になりますが、マニュアルでは個体間変動の大きさがCV％（変動係数：標準偏差 / 平均 × 100）で書かれているため、ここで指数誤差モデルにおけるωとCVの関係についても触れておきます。
導出は割愛しますが、対数正規分布の平均と分散よりCVは以下のように表されます。
　
　
⇔　
　
が小さい場合はと近似できるため ですが、今回はVpのCVが70%を超えており、その他のCL、Vcについては30%を超える値となっています。CV=0.3程度であれば近似とのずれは2%くらいですが、一貫性の観点から今回のプログラムではCL、Vc、VpのCVに対して上記の変換を行います。
　
　
【２】 PKパラメータの変動が確率分布を持つと仮定したのと同様に、個体内での変動や測定誤差等にも確率分布を仮定します。薬物血中濃度の測定値は、コンパートメントモデルにより推定される薬物血中濃度にこれらの誤差が加わって、実際の測定値となっていると仮定します。このようにすることで、も確率変数として扱うことができます。
このような誤差による変動をまとめて残差変動（個体内変動）とし、この誤差を表す確率変数をとすると、は平均0、分散の正規分布に従うと仮定します。

以前の記事内容ではここでも指数誤差モデルを採用してしまいましたが、この残差変動についてはPKパラメータのように観測値が指数関数的に幅広い範囲で分布するわけではないため、通常は指数誤差モデルは用いられません。
この残差変動のσもPATのパラメータを使用しますが、今回は以下の比例誤差モデルを採用することとします。
　

　
この場合のCVとσの関係についてもみてみます。
なので、比例誤差モデルにおけるは平均、分散 の正規分布に従うことになります。
（確率変数の線形変換の公式 より）

よって、このモデルにおける変動係数CVは
　

　
であり、比例誤差モデルではCVをσとしてそのまま使えます。
　
　

ベイズの定理

では次に実際にどのような計算を行えば良いのかを見ていきます。以前ナイーブベイズというアルゴリズムで文書分類を行う記事を書いたことがありましたが、考え方としては基本的には同じです。
www.yakupro.info

ベイズの定理を再掲します。
D：データ、θ：パラメータとすると

• P(θ)：事前分布(Prior Distribution)データを見る前のθが元々持っている分布
• P(D|θ)：尤度(Likelihood)あるパラメータθを仮定したときに、データDが得られる確率です。そのパラメータの尤もらしさを表します。
• P(D)：周辺尤度(marginal likelihood)どのようなθであってもデータの得られる確率を表す定数（正規化定数）です。
• P(θ|D)：事後分布(Posterior Distribution)と呼ばれる、データを見た後のθの分布。つまり、実際の観測データによって事前分布が修正されて得られる分布です。

今回のケースでは、データ（実測濃度）を観測した後に個人のCLやVdなどのパラメータを推定したいのでした。つまり、P(θ|D)を最大にするようなθが得られればよいのであって、具体的なP(θ|D)のそのものの値は今回は必要ではありません。従って、P(D)は無視することができ（θに依存しないため）、次の式のように分子の部分だけを評価することになります。*1

この最適化は最大事後確率（maximum a posteriori：MAP）推定といい、本記事でのTDMにおける「ベイズ推定」や「ベイジアン法」とは、このMAP推定のことを指します。また、この事後確率を最大化させるパラメータをMAP推定量といいます。
次は各項について一つずつみていきます。

・尤度
尤度（Likelihood）は、ある特定のパラメータを仮定したときの観測データの発生確率（もっともらしさ）を表します。
前述の比例誤差モデルでは、実測値は以下の正規分布に従います。

また、ここではコンパートメントモデルによる血中濃度の予測値であり、このはの関数です。（CL, Vdなどはηから算出される値です。）

これらを踏まえた上で、互いに独立な実測血中濃度データがN個与えられたとき、これらN個のデータ点の同時確率（密度）は次のようになります。
　

　
これは正規分布の式と同じですが、観測データではなくパラメータの関数であるため尤度関数といい、以下のように明示的にで表します。
　

　
このように全ての観測時点の確率（尤度）をかけ合わせることで、観測されたデータを最も矛盾なく説明できるを探し出すための尤度関数が定義されます。
　
・事前分布
今回の場合、事前分布は観測データが与えられる前のの確率分布であり、, , の3つあります。
例として、CLではなので
　

他のηVc, ηVpについても同様で、これらは独立して発生すると仮定するため各、、の積となります。
　
・事後分布
事後分布は次のように表されます。

事後分布を最大化するパラメータを探すことが目標でしたが次に示す理由から、この式の対数を取り、さらに-2をかけたものを考えます。
（※-1、-2のどちらをかけても目的関数の値が変わるだけで、推定されるηの値は同じになります。）

• 対数を取ると小さな数の積を和とすることができ、扱いやすくなる
• 対数関数は単調増加でありマイナスをかけることは大小を逆転させるため、最大化問題を最小化問題へ変換することができる

では尤度関数の自然対数をとって-2をかけます。（指数部分は整理せず）
　

　
最小化に関係のない定数部分を無視して整理すると、第１項部分はが残ります。しかし、この項があると目的関数が複雑になり収束が不安定になりやすい、この項を含まない近似としてMAP推定を行うことが標準として定着している等の理由*2から、本記事でもこの項を除いて目的関数を定義します。
　
また、各それぞれの事前分布にも同様にを適用すると
　

　
この第１項はηの値に関係のない定数項なので無視します。今回パラメータは３つ（CL、Vc、Vp）あるので、最終的に以下のようになります。
　

　
これらをまとめて、最小化する関数をOBJ（Objective Function）として表すと以下のようになります。
　

　
ある患者の実際に得られたデータ（実測血中濃度）から、薬物動態パラメータ（CLやVd）を推定するという問題は、ベイジアン法によって、この目的関数を最小化するパラメータηを探すという問題に帰着しました。
　
Pythonで実装した関数名はobjective()としました。この関数では、与えられたηからPKパラメータを算出し、それを元に生成したインスタンスでC_pred（予測血中濃度）を計算し、目的関数の値を返すという処理を行います。

def objective(eta, pop_params, variability, dosage, obs_times, c_obs):
    eta_cl, eta_vc, eta_vp = eta
    CL_pop, Vc_pop, Vp_pop, Q_fixed = pop_params
    omega_CL, omega_Vc, omega_Vp, sigma = variability
 
    # 個体パラメータの算出
    CL_i = CL_pop * np.exp(eta_cl)
    Vc_i = Vc_pop * np.exp(eta_vc)
    Vp_i = Vp_pop * np.exp(eta_vp)
    
    # 事前分布項
    log_prior = eta_cl ** 2 / omega_CL ** 2 \
              + eta_vc ** 2 / omega_Vc ** 2 \
              + eta_vp ** 2 / omega_Vp ** 2
    
    param_i = [CL_i, Vc_i, Vp_i, Q_fixed]
    model = TwoCompInfModel(*param_i, **dosage)
    
    # 尤度項（比例誤差モデル）
    log_likelihood = 0
    for c, t in zip(c_obs, obs_times):
        c_pred = model.get_conc_at_time(t)
        log_likelihood += (c - c_pred) ** 2 / (c_pred *sigma) ** 2 
        
    return log_likelihood + log_prior

実測濃度情報の入力

お手本としたPATではベイズ推定を行う際には「実測濃度(μg/mL)、採血タイミング、n回投与後」の3項目を入力します。作成するプログラムでも同様の方式にしたいと思います。
以下のようにこれら3項目をタプルとしたリストで与えると、投与開始後の実測血中濃度を得た時間とその濃度に変換し、目的関数と共に最適化ソルバーに渡します。
例えば、3回投与後（4回目直前）のトラフの実測濃度が10μg/mL、4回投与後の投与終了1時間後のピーク濃度が25μg/mLの場合は以下のように与えます。

# 実測血中濃度情報の入力
# (実測濃度[μg/mL]、採血タイミング、n回投与後)
cobs_info = [
    (10, 'trough', 3),
    (25, 'c1', 4)
]
# pop_modelは母集団平均のインスタンス
obs_times, c_obs = create_cobs_info_list(pop_model, cobs_info)    
print(f'採血時間:{obs_times}, 実測血中濃度:{c_obs}') # [36, 38.0] [10, 25]

「c1」の意味は「投与終了1時間後」という意味です。もちろん「c1.5」や「c2」などの入力でもOKです。このへんの仕様もPATのルールに倣いました。
関数create_cobs_info_list()は与えられたn回投与後の'trough'や'c1'等の文字列を投与開始後の経過時間に変換する関数です。例えば、投与間隔が[12, 12, 12, 12, 12]で2回目以降の点滴時間が1時間のとき、3回目投与後の'trough'なら36(hr)、4回目投与後の'c1'なら38(hr)、といった具合です。

Scipyによる目的関数の最小化

最小化すべき目的関数は決まりましたが、これを最小にするようなηを数式変形だけで直接（解析的に）求めることはできません。今回はこの非線形最適化問題を数値的に解くためにScipyのoptimize.minimize()を使いました。

def calc_bounds(omegas, multiplier=3):
    return [(-omega * multiplier, omega * multiplier) for omega in omegas]
 
# 個体間変動（CV→ωに変換）
omega_CL = np.sqrt(np.log(1 + 0.369 ** 2))
omega_Vc = np.sqrt(np.log(1 + 0.333 ** 2))
omega_Vp = np.sqrt(np.log(1 + 0.728 ** 2))
 
# 残差変動（CV=σ）
sigma = 0.181
 
variability = [omega_CL, omega_Vc, omega_Vp, sigma]
initial_eta = [0, 0, 0]
multiplier = 3.0  # ±3ωに設定
bounds = calc_bounds([omega_CL, omega_Vc, omega_Vp], multiplier)
 
result = minimize(objective,
              initial_eta,
              method='Powell',
              bounds=bounds,
              args=(params, variability, dosage, obs_times, c_obs))

最適化アルゴリズムはいろいろ試しましたが、最終的にはPowellを選びました。数値勾配を用いるL-BFGS-BやSLSQPでは今回の実装では数値微分の精度が安定しないケースがありましたが勾配不要のPowellでは結果が安定しており、またPATの結果に最も近い推定値が得られたためこれを採用することにしました。なお、ηの初期値についてはすべて0（母集団平均から個体差なしの状態）とし、変数の境界制約（bounds）は±3ω以内とすることにしました。

AUCガイドによる推奨投与量の計算

推定されたηからPKパラメータのCL・Vc・Vpを計算し、その患者の個別化モデルを構築します。この推定モデルを用いて、目標とする定常状態のAUCを達成するための１日投与量を算出します。
定常状態において、AUCは投与量に対して線形（比例）関係なので以下の計算で求めることができます。
　

　
（※ ）
　
コードでは目標AUCはデフォルトで500μg・h/mLに設定していますが、自由に変更できるようにしています。
　
算出された1日投与量をどのような投与回数で振り分けるかは実際のTDMソフトでは複雑なアルゴリズムが使われていると推定されますが、今回はプログラムの簡便さを優先して以下のルールベースで行うことにしました。

• Ccr ≦30のときは24時間ごと、それ以外は12時間ごとの投与とする
• １回投与量は100mg単位となるように四捨五入で丸める
• 点滴時間は投与量によらず一律１時間とする
• トラフ値の上限確認は行わない

あくまでAUCガイドを目安とした推奨投与量の試算として行います。

テスト症例による実行結果

では実際にプログラムを動かしてみます。以下の例はPAT ver. 4.1の操作マニュアル記載のモデルケースを使用しています。

• 症例：60歳 男性、体重 60 kg、血清クレアチニン値 0.8 mg/dL
• 投与計画：初回負荷投与量として 1800mg（点滴時間2hr）、その12時間後の２回目の投与から1回850mg（点滴時間1hr）を12時間ごとに3回投与
• 血中濃度測定：3回目の投与後（4回目投与直前）のトラフ値に 10 μg/mL、4回目の投与終了後1時間後に25 μg/mL

実測濃度情報の入力方法については上の項目で示した通りです。患者背景と投与条件の入力方法につきましては、前回記事に詳しく書きましたので、よろしければご参照ください。
www.yakupro.info

上記の患者背景、投与条件、実測濃度2点をコード中に書いてからプログラムを実行すると、まず最初は母集団平均による時間-濃度プロファイルが描画されます。
グラフ右下に"Bayesian estimation"のボタンを作成しているので、これをクリックすると目的関数の最小化が実行され、ベイズ推定PK値によるグラフ描画や推奨投与量(AUC for 500)の計算が行われます。

以下は私の実行環境（Python3.12, scipy1.17.1, Windows11）での出力結果です。*3（画像はページトップと同じ）

Estimated eta_cl: 0.134264989552591
Estimated eta_vc: -0.0009923348823843337
Estimated eta_vp: 0.12538711569031916
Estimated CL: 3.949 L/h
Estimated Vc: 32.068 L
Estimated Vp: 65.975 L
 message: Optimization terminated successfully.
 success: True
  status: 0
     fun: 0.20849831945277508
       x: [ 1.343e-01 -9.923e-04  1.254e-01]
     nit: 4
   direc: [[ 0.000e+00  0.000e+00  1.000e+00]
           [ 0.000e+00  1.000e+00  0.000e+00]
           [-9.761e-03 -2.723e-03  2.217e-02]]
    nfev: 223

funは最小化された目的関数の値（Objective Function Value）で、xはそのときのパラメータの値です。PATでベイズ推定の際にグラフ右上に表示される”OFV”はこのfunと同じものを指しています（実装の違いにより数値は一致するとは限りません）。今回のプログラムでも同様に表示されるようにしました。
他の出力項目については参考情報ですが、以下の通りです。

• nit（Number of iterations）：オプティマイザが実行した反復回数
• nfev（Number of function evaluations）：目的関数が評価された回数
• direc（direction matrix）：Powellが使用した探索方向（通常は参照不要）

　
比較したのはWeb版PAT ver. 4.1a（2026年6月現在）です。
ベイズ推定によるPKパラメータの推定値はPATとほぼ一致する結果となりました。

項目	PAT ver. 4.1a	本プログラム	差
CL（L/h）	3.93	3.95	+0.02
Vd/Vss（L）	99.29	98.04	-1.25
半減期（hr）	5.66	5.63	-0.03
AUC（μg・h/mL）	432.13	430.51	-1.62
トラフ値（μg/mL）	12.02	11.94	-0.08

推奨投与量（for AUC 500）も同じに計算され1000 mg１日２回となり、この投与量によるAUCssはPATでは508.39μg・h/mL（トラフ14.14μg/mL）に対して、本記事のプログラムでは506.48μg・h/mL（トラフ14.05μg/mL）でした。
推定CLが得られれば推奨1日投与量 = 目標AUC × CLから直接計算できるため、症例に応じてAUCへの適合を優先するか調剤しやすい量を優先するかを考慮して最終的な投与量を決定することができます。

事前分布を使う意味

今回導出された目的関数の尤度の項は実測とモデルによる予測の差（誤差）の大きさを返します。ということは、予測と実測が近ければ近いほど、目的関数の値は小さくなります。また、事前分布の項はCL等の個人のPKパラメータが母集団平均に近ければ近いほど、目的関数の値は小さくなり、反対に離れれば離れるほど、目的関数の値は大きくなります。
　
目的関数を最小化するパラメータηとは、尤度の項についてみれば実測値に沿ったシミュレーションカーブを描けるような、平均から離れたパラメータが良いわけですが、そうすると今度は事前分布の項が大きくなってしまい最小化からは遠ざかってしまいます。
これは見方を変えると、事前分布の項があることでパラメータが極端な値を取ることにブレーキをかけていることになります。

もし事前分布（=正則化項）がなかったら

下に事前分布の項がない場合に、実測値に対するフィッティングがどのように変わるかを試したグラフを載せます。
先ほどのテスト症例において、入力データはまったく変えずに事前分布項のあり・なしで比較してみました。

上図が事前分布ありの場合、下図が事前分布なしの場合の最適化の結果です。
下図では実測データ点に対してぴったりと追従*4することによってOFV=0となり、過剰適合（オーバーフィッティング）していることが見て取れます。Vssは「あり」の場合の97.71から「なし」では47.09と半分近くの不自然な値になり、CLも3.96→4.25となってこれにより推奨投与量にも影響を与えています。

MAP推定は最尤推定（尤度関数だけを使った推定）に過学習を防ぐための工夫である正則化（ペナルティ）を加えたものとみることができます。そしてこの正則化の強さは、ωとσの相対的な大きさによって決定されます。具体的には、ωが大きいほど正則化が弱くなり、σが大きいほど相対的に正則化が強くなります。

正則化については、以前記事で書いたことがありますのでよろしければご参照ください。
www.yakupro.info

2-コンパートメントモデルのような複雑なモデルを、実測値1点でフィッティングするにはモデルの表現力に対して、データの数が少なすぎるということが言えます。この場合は事前分布の項の影響が相対的に大きくなり、パラメータは母集団平均に引き寄せられることになります。
反対に実測データの数が多い場合を考えてみます。
実測データの数が多ければ、目的関数の最小化には尤度の寄与が大きくなり、平均から離れた値でも取りやすくなります。したがって、実測データが新しく得られるたびにその患者固有の推定値に近づいていくことになります。しかし、単に実測データが多ければそれで良いのかというとそうではありません。

採血タイミングが推定に与える影響

採血タイミングがトラフなのかピークなのかによって*5血中濃度が持つ情報には違いがあり、推定PK値にもそれが反映されます。

ピーク値は投与後の早い時間の血中濃度です。つまり薬物の初期の体内分布を反映するため、分布容積（2-コンパートメントモデルでは主にVc、1-コンパートメントモデルではVd）の情報が得られます。一方トラフ値は投与直前、多くの薬物が体内から排出された後の血中濃度であるためCLの情報が得られます。反対に、分布容積に関する情報は少なくなります。（少ないと言っても、含まれていないわけではありません）

どちらか一方だけの採血タイミングでは、分布容積、CL両方のパラメータを正確に推定することは無理があるということです。
今回のモデルには推定すべき個体差（）が3つ（・・）あります。前述したように採血タイミングによって得られる情報の種類が異なるため、異なる時間点の実測データが複数あることでより推定PK値の精度が高まります。
推定するのに情報が不足しているPK値は、結局事前分布に依存した推定値（＝母集団平均に近い値）となり、患者個別の情報を反映した結果とは言えなくなります。

TDMソフトのベイズ推定を使う上で薬剤師が知っておきたいこと

まとめになります。

• 事前分布は結果に大きな影響を与えています。対象患者が事前分布（母集団）から大きく離れていると思われる場合には推定値の信頼性が低下する可能性があるため、そのことを踏まえておく必要があります。

• トラフ値だけでなくピーク値も用いた2点で推定を行う方がAUC算出の正確性は高まります。トラフ値からだけの推定を行う場合は、算出されるAUCには2点推定よりも大きな誤差が含まれる場合があることを考慮しておかなければなりません。

参考

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

月刊薬事2022年05月号 [雑誌] (特集:抗菌薬TDM臨床実践ガイドライン2022はここに注目! 新旧抗菌薬の使い方)

母集団薬物デ-タの解析 (統計科学選書 6)

最後に今回作成したプログラムです。

from decimal import Decimal, ROUND_HALF_UP
from scipy.integrate import trapezoid
from scipy.optimize import minimize
import japanize_matplotlib
import matplotlib.pyplot as plt
import matplotlib.widgets as wg
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)       
        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):
        """台形公式の自作実装. 記事(1)で使用"""
        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ごとに合計）. 記事(1)で使用"""
        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 get_conc_at_time(self, time_hr):
        return self.conc[int(time_hr * self.STEP)]
 
    def convert_timing_to_hours(self, timing, num_doses=0):
        # num_dosesは全投与回数を超えないように
        if type(timing) is str:
            if timing.casefold() == 'trough':
                return np.sum(self.tau_list[:num_doses])
            else:
                t = float(timing[1:])  # delete 'c'
                time = np.sum(self.tau_list[:num_doses - 1]) + \
                       self.inf_t_list[num_doses - 1] + t
                return time
        return float(timing)
 
 
def create_cobs_info_list(model, arg):
    obs_times, c_obs = [], []
    for info in arg:
        c, timing, n_doses = info
        t = model.convert_timing_to_hours(timing, n_doses)
        obs_times.append(t)
        c_obs.append(c)
    return obs_times, c_obs
 
  
def calc_bounds(omegas, multiplier=3):
    return [(-omega * multiplier, omega * multiplier) for omega in omegas]
  
 
def objective(eta, pop_params, variability, dosage, obs_times, c_obs):
    eta_cl, eta_vc, eta_vp = eta
    CL_pop, Vc_pop, Vp_pop, Q_fixed = pop_params
    omega_CL, omega_Vc, omega_Vp, sigma = variability
 
    # 個体パラメータの算出
    CL_i = CL_pop * np.exp(eta_cl)
    Vc_i = Vc_pop * np.exp(eta_vc)
    Vp_i = Vp_pop * np.exp(eta_vp)
    
    # 事前分布項
    log_prior = eta_cl ** 2 / omega_CL ** 2 \
              + eta_vc ** 2 / omega_Vc ** 2 \
              + eta_vp ** 2 / omega_Vp ** 2
    
    param_i = [CL_i, Vc_i, Vp_i, Q_fixed]
    model = TwoCompInfModel(*param_i, **dosage)
    
    # 尤度項（比例誤差モデル）
    log_likelihood = 0
    for c, t in zip(c_obs, obs_times):
        c_pred = model.get_conc_at_time(t)
        log_likelihood += (c - c_pred) ** 2 / (c_pred *sigma) ** 2
        
    return log_likelihood + log_prior
 
 
def calc_recommended_dose(tdd_mg, ccr):
    # 投与間隔決定
    tau = 24 if ccr <= 30 else 12
    n_per_day = 24 // tau  # 1日投与回数
    
    # 1回投与量を100mg刻みで丸める
    single_dose = int(tdd_mg / n_per_day / 100 + 0.5) * 100
 
    return single_dose, tau, n_per_day
 
 
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 disp_recommended_dose(model, ax, new_dose, tau, 
                       target_auc=500, color='purple'):
    text = (
        f'●ベイズ推定 投与量 (for AUC {target_auc}): '
        f'{round_v(new_dose, "1")} mg, '
        f'{round_v(tau, "1")} h 毎に, '
        f'AUC={round_v(model.calc_ss_auc())} μg h/mL, '
        f'トラフ値={round_v(model.calc_css_trough_analytic())} μg/mL'
    )
    ax.text(0.27, 1.06, text, 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, '腎機能推定誤差が大きい可能性あり. 推定値に注意!',
                transform=ax.transAxes, fontsize=12, color='red')
    # トラフ値ラベル
    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
 
    # 投与条件の設定. 不均等投与に応じて作成
    # [投与量(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)
    
    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, '〇母集団平均')
    ax2 = fig.add_axes([0.8, 0.05, 0.15, 0.075])
    button = wg.Button(ax2, 'Bayesian estimation')
 
    # 実測濃度情報の入力
    # - 実測血中濃度[μg/mL]
    # - 採血タイミング ('trough', 'c1', 'c1.5',.., または時間[hr]で入力)
    # - n回投与後
    cobs_info = [
        (10, 'trough', 3),
        (25, 'c1', 4)
    ]
 
    obs_times, c_obs = create_cobs_info_list(pop_model, cobs_info)
    print(f'採血時間:{obs_times}, 実測血中濃度:{c_obs}')
 
    # 個体間変動（CV→ωに変換）
    omega_CL = np.sqrt(np.log(1 + 0.369 ** 2))
    omega_Vc = np.sqrt(np.log(1 + 0.333 ** 2))
    omega_Vp = np.sqrt(np.log(1 + 0.728 ** 2))
    
    # 残差変動（CV=σ）
    sigma = 0.181
    
    variability = [omega_CL, omega_Vc, omega_Vp, sigma]
    initial_eta = [0, 0, 0]
 
    def button_click(event):
        multiplier = 3.0  # ±3ω
        bounds = calc_bounds([omega_CL, omega_Vc, omega_Vp], multiplier)
        result = minimize(objective,
                          initial_eta,
                          bounds=bounds,
                          method='Powell',
                          args=(params, variability, dosage, obs_times, c_obs)
                          )
        eta_cl_est, eta_vc_est, eta_vp_est = result.x
        minimized_ofv = result.fun
        
        # 最適化されたパラメータを使って血中濃度のグラフを描画
        CL_est = CL_pop * np.exp(eta_cl_est)
        Vc_est = Vc_pop * np.exp(eta_vc_est)
        Vp_est = Vp_pop * np.exp(eta_vp_est)
        
        print(f"Estimated eta_cl: {eta_cl_est}")
        print(f"Estimated eta_vc: {eta_vc_est}")
        print(f"Estimated eta_vp: {eta_vp_est}")
        print(f"Estimated CL: {round_v(CL_est, '0.001')} L/h")
        print(f"Estimated Vc: {round_v(Vc_est, '0.001')} L")
        print(f"Estimated Vp: {round_v(Vp_est, '0.001')} L")
        print(result)
 
        est_params = [CL_est, Vc_est, Vp_est, Q_fixed]
        est_model = TwoCompInfModel(*est_params, **dosage)
 
        for obs_time, obs in zip(obs_times, c_obs):
            ax.scatter(obs_time, obs, color='black')
        ax.text(0.99, 0.99, 'OFV=' + round_v(minimized_ofv),
                transform=ax.transAxes, ha='right', va='top')
 
        color = 'green'
        plot_conc(est_model, ax, color)
        disp_trough_values(est_model, ax, -7, color)
        disp_pk_result(est_model, ax, -0.1, '●ベイズ推定', color)
        
        # 目標AUCのための推奨1日投与量を計算
        new_tdd = target_auc * est_model.CL
        single_dose, tau, n_per_day = calc_recommended_dose(new_tdd, ccr)
        sim_count = n_per_day * 3
        new_dosage = {'new_dosage': [single_dose, 1, tau, sim_count]}
        
        new_model = TwoCompInfModel(*est_params, **new_dosage)
        disp_recommended_dose(new_model, ax, single_dose, tau, target_auc)
        fig.canvas.draw()
 
    button.on_clicked(button_click)
    plt.show()

*1:記号∝は比例関係を表します

*2:自分調べ

*3:環境によって最適化の収束値がわずかに異なる場合があります

*4:Python上のグラフを拡大すると下図は完全に実測点に一致していることが確認できます

*5:より正確に言うと、投与後の経過時間によって