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

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

トップ > 薬物動態 > 抗菌薬でモンテカルロシミュレーション①－乱数を用いたPKパラメータの生成－

抗菌薬でモンテカルロシミュレーション①－乱数を用いたPKパラメータの生成－

薬物動態

はじめに

前回、前々回とPK/PDパラメータを求めるコードを自前のクラスに実装したので、今度はこれらのパラメータを使うモンテカルロシミュレーションに挑戦したいと思います。
ネットの情報や本で調べた範囲で書いた記事で、実際にソフト等を使って結果を検証したわけではありません。私自身、理解が曖昧なまま書いている部分もありますので、その点ご留意ください。
今回の記事では、手順や方法はこちらの文献（薬理学的観点からみたPK-PD理論とブレイクポイント）を参考にさせてもらいました。詳しくはリンク先をご参照ください。

モンテカルロシミュレーションとは

モンテカルロ法（別名：モンテカルロシミュレーション）とは、『数値計算手法の一つで、乱数を用いた試行を繰り返すことにより近似解を求める手法。確率論的な事象についての推定値を得る場合を特に「モンテカルロシミュレーション」と呼ぶ。』（IT用語辞典 e-Words）
抗菌薬においては、患者のPKパラメータや原因菌のMICの分布を乱数で発生させることにより、血中濃度推移からPK/PDパラメータを得ることで有効率や治療効果を予測することになります。


目的と手順の確認

今回の目的は抗菌薬のメロペネムにおいて、一定の基準以上の有効率が得られるMICの最大値をPK/PDブレイクポイントとして設定することです。
ブレイクポイントとは抗菌薬の有効性の判定基準となるMIC値のことで、米国のCLSI、日本化学療法学会のものなど様々あるようです。今回はPK/PD理論を用いたブレイクポイントということになります。
では手順の概要を見てみます。

  1. メロペネムの母集団PKパラメータの平均と分散から10,000例のPKパラメータのセットを発生させる。
  2. 1回の用量、1日の投与回数、点滴時間を設定する。
  3. 10,000例のメロペネムの薬物血中濃度推移を得る。
  4. MICを0.063μg/mLから128μg/mLまで順次変更させ、各MICにおける%T>MICを求める。
  5. この%T>MICのうち、目標値(40%)を満たした例数からPK/PD目標達成確率を算出する。

他、細かいところとしてはPK/PDブレイクポイントは血中濃度ベースの値とする、想定する患者としてCCrが（高い・低い）、体重が（大きい・小さい）の組み合わせを用いる、等ありますが、この辺は順にコードを書いていく中で見ていきます。

乱数を用いてPKパラメータを発生させる

メロペネムの母集団PKパラメータとして、以下の表のデータを用います。

fixed effects parameter（θ） interindividual variability（ω）
CL 0.0905×CCr＋2.03 41.1%
Vc 0.199×Wt 39.8%
Q 4.02 32.8%
Vp 4.55 29.9%
CL：全身クリアランス（L/h）
Vc：中心コンパートメントの分布容積（L）
Q：コンパートメント間のクリアランス（L/h）
Vp：抹消コンパートメントの分布容積（L）
CCr：クレアチニンクリアランス（mL/min）
Wt：体重（kg）

ここで、fixed effectsは固定効果と訳され、母集団の平均値を指すようです。また、 これらのパラメータには変動部分があり、それをinterindividual variability（個体間変動）と呼ぶとのこと。この変動（誤差）を表す分布として、今回は以下の式に従う対数正規分布を仮定します。

P_i = \overline{P} × \exp(η_i)

P_i：個人のPKパラメータ
\overline{P}：母集団PKパラメータの平均値
η_i：平均0, 分散ω^2の正規分布に従う誤差

つまり\ln{(Pi)}は平均\ln{(\overline{P})}、分散ω^2の正規分布に従う、ということです。（対数を取ったものが正規分布になる）

ではここまでをコードで書きます。

import numpy as np
 
sample_size = 10000
CCr = 100      # クレアチニンクリアランス(mL/min)
Wt = 40        # 体重(kg)
 
CL = 0.0905 * CCr + 2.03
Vc = 0.199 * Wt
 
np.random.seed(1234)   # seedを設定し発生する乱数を固定する.
 
X = np.random.randn(4, sample_size)     # 引数で指定した形状の標準正規乱数を生成
 
X_cl = CL * np.exp(X[0, :] * 0.411)      # CL(L/h)
X_vc = Vc * np.exp(X[1, :] * 0.398)      # Vc(L)
X_q = 4.02 * np.exp(X[2, :] * 0.328)     # Q(L/h)
X_vp = 4.55 * np.exp(X[3, :] * 0.299)    # Vp(L)


乱数の生成はnumpyのrandomモジュールを使用します。Python標準にも乱数生成モジュールはありますが、発生させた大量のデータを後でまとめて処理するのにnumpyは高速でコードも簡潔に書くことができ、都合が良いからです。
また、今どきのメジャーな言語同様、乱数生成アルゴリズムに周期の長いメルセンヌ・ツイスタが使われており、乱数の質としても問題ありません。

まずサンプルサイズ、CCr、Wtを設定し、次にnp.random.seed(1234)で疑似乱数のジェネレータを初期化しています。乱数の出力を同じにしたい場合はこのようにseedを設定（数字は何でもかまいません）する必要があります。
np.ramdom.randn()は標準正規分布に従う乱数を返します。ここでは引数を指定して（４行×サンプルサイズ列）のnumpy配列として取得しています。
その後、X_cl～X_vpの行のnp.exp()内で、発生させた平均０、分散１の標準正規分布にωをかけて平均０、分散ω^2の正規分布に変換しています。
ここで少し補足説明します。対数正規分布LN(μ, σ^2)（このμ、 σ^2はこの分布自体の値ではなく、正規分布となるln(X)の平均と分散。紛らわしい）に従う確率変数Xに対し、平均E(X)、分散V(X)は以下の式で与えられます。

E(X) = e^{μ＋σ^2/2}
V(X) = (e^{σ^2}-1)e^{2μ＋σ^2}

またこのＸの変動係数の２乗CV^2をこれらから求めるとe^{σ^2} - 1となります。
ここで、CV\ll1のときe^xのマクローリン展開よりe^{σ^2} ≒ 1 + σ^2と近似でき、CV ≒ σとなります。私は最初なぜ変動係数のωを標準偏差として扱っているかわからなかったので、ここに記しておきます。
np.randomモジュールには対数正規分布からサンプリングするlognormal()という関数があります。これを使って、

 
X_cl = np.random.lognormal(np.log(CL), 0.411, sample_size)

としても同じです。引数は前から平均、標準偏差、サイズですが、平均と標準偏差は基となる正規分布の値です。こちらの方がシンプルですが、上の方で数式を書いたので、コードはそれに沿ったやり方で書きました。

では発生させたPKパラメータのうち、ＣＬの分布を見てみます。
グラフ描画部分のコードです。

from scipy.stats import norm

# PKパラメータCLの分布を確認
plt.figure()
 
# X_clの分布(対数正規分布)
plt.subplot(1, 2, 1)    # (行数, 列数, 何番目のプロットか)
# (データ配列, bins:ビン(棒)の数, density:Trueにすると面積の合計が1になるよう正規化, alpha:透明度, ec: 棒の縁の色(edgecolor)
plt.hist(X_cl, bins=50, density=True, alpha=0.5, color='green', ec='black')
plt.xlabel('X_cl')
 
# ln(X_cl)の分布(対数値が正規分布になる)
plt.subplot(1, 2, 2)
x = np.arange(0, 5, 0.01)
y = norm.pdf(x, np.log(CL), 0.411)    
plt.plot(x, y, color='red', linewidth='2', alpha=0.8, label='μ=ln(CL), σ=0.411')
plt.hist(np.log(X_cl), bins=20, density=True, alpha=0.5, color='blue', ec='black')
plt.xlabel('ln(X_cl)')
plt.legend()
 
# 座標軸の一覧取得
axs = plt.gcf().get_axes()
for ax in axs:
    # current axsを設定
    plt.sca(ax)
    # 軸ラベルを表示
    plt.ylabel('density')
    # グリッドの表示
    plt.grid(which='both', linestyle='dashed')
 
print('ln(X_cl)の平均: ', np.mean(np.log(X_cl)))
print('ln(X_cl)の標準偏差: ', np.std(np.log(X_cl)))
print('X_clの変動係数: ', np.std(X_cl) / np.mean(X_cl))
plt.show()


f:id:enokisaute:20200307144808p:plain

ln(X_cl)の平均:  2.411769656398536
ln(X_cl)の標準偏差:  0.40901894119910176
X_clの変動係数:  0.42030981923847915

左側の緑のヒストグラムが発生させたCLの分布、右側がその対数を取ったものです。右側の赤いラインは平均μ=\ln{(\overline{CL})}、標準偏差σ＝0.411の正規分布のグラフです。対数値が正規分布になっていることが見て取れます。X_clの変動係数ωはある程度の大きさがあるせいか、ln(X_cl)の標準偏差とは少し差がありますね。

長くなってきたので今回はこのへんで。続きます。
今回のコードを載せておきます。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
 
# 抗菌薬でモンテカルロシミュレーション1

sample_size = 10000
CCr = 100       # クレアチニンクリアランス
Wt = 40         # 体重
 
CL = 0.0905 * CCr + 2.03
Vc = 0.199 * Wt
 
# seedを設定し発生する乱数を固定する.
np.random.seed(1234)
 
# 引数で指定した形状(4×sample_size)の標準正規乱数を生成
X = np.random.randn(4, sample_size)
 
# 各PKパラメータの分布は対数正規分布を仮定する
X_cl = CL * np.exp(X[0, :] * 0.411)       # CL(L/h)
X_vc = Vc * np.exp(X[1, :] * 0.398)       # Vc(L)
X_q = 4.02 * np.exp(X[2, :] * 0.328)      # Q(L/h)
X_vp = 4.55 * np.exp(X[3, :] * 0.299)     # Vp(L)
 
# PKパラメータ(CL)の分布を確認
plt.figure()
# X_clの分布(対数正規分布)
plt.subplot(1, 2, 1)    # (行数, 列数, 何番目のプロットか)
# (データ配列, bins:ビン(棒)の数, density:Trueにすると面積の合計が1になるよう正規化, alpha: 透明度, ec: 棒の縁の色(edgecolor)
plt.hist(X_cl, bins=50, density=True, alpha=0.5, color='green', ec='black')
plt.xlabel('X_cl')
 
# ln(X_cl)の分布(対数値が正規分布になる)
plt.subplot(1, 2, 2)
x = np.arange(0, 5, 0.01)
y = norm.pdf(x, np.log(CL), 0.411)  # 正規分布の確率密度関数
plt.plot(x, y, color='red', linewidth='2', alpha=0.8, label='μ=ln(CL), σ=0.411')
plt.hist(np.log(X_cl), bins=20, density=True, alpha=0.5, color='blue', ec='black')
plt.xlabel('ln(X_cl)')
plt.legend()
 
# 座標軸の一覧取得
axs = plt.gcf().get_axes()
for ax in axs:
    # current axsを設定
    plt.sca(ax)
    # 軸ラベルを表示
    plt.ylabel('density')
    # グリッドの表示
    plt.grid(which='both', linestyle='dashed')
 
print('ln(X_cl)の平均: ', np.mean(np.log(X_cl)))
print('ln(X_cl)の標準偏差: ', np.std(np.log(X_cl)))
print('X_clの変動係数: ', np.std(X_cl) / np.mean(X_cl))
plt.show()


参考

1. 薬理学的観点からみた PK-PD 理論とブレイクポイント
2.モンテカルロシミュレーションに対応したMicrosoft ® Office Excel による抗菌薬の PK / PDシミュレーションソフトの開発