抗菌薬でモンテカルロシミュレーション③ーPK/PDブレイクポイントを求めるー

一定の基準以上の有効率が得られるMICの最大値をPK/PDブレイクポイントとして求めます。
内容的にはこちらの続きとなります。
www.yakupro.info

目標達成確率の算出
目標達成確率-MICのグラフを作成する
PK/PDブレイクポイントを確認する
プログラムの処理時間について
参考

目標達成確率の算出

前回までで各MICにおける10,000例の%T>MICをdata_matrixに取得することができました。次は、メロペネムの最大殺菌効果発現に必要な40%を目標値として、その目標値以上の例数をカウントします。この例数の割合をもって『PK/PDパラメータ目標達成確率(%)』とし、その値を算出します。

# key: 投与方法(dosage)
# value: (%T>MICが40以上となった数 / sample_size * 100, グラフ表示マーカー)
result = {}
 
for dosage, value in recipe.items():
    #######################
    # A, Bを求める処理(省略)
    #######################
    for si in range(sample_size):
    bc = BloodConcTwoComp2(A[si], alpha[si], B[si], beta[si], T, intervals)
        for mi, mic in enumerate(mic_list):
            data_matrix[mi, si] = bc.get_time_above_mic(mic)
 
    result[dosage] = (np.sum(data_matrix >= 40, axis=1) / sample_size * 100, value[-1])

一番外側のfor文の外で辞書型のresultを宣言しておきます。投与方法（dosage）をkeyとして、目標達成確率(%)、グラフ表示マーカーのタプルをvalueに取ります。また、np.sumのところで40以上の%T>MICをカウントしています。
data_matrixは２次元配列なので、axis=1と指定することで行（MIC）ごとの処理となります。具体的な例で見てみます。

>>> data = np.array([[55, 42, 40, 37], [47, 43, 39, 32]])
>>> data
array([[55, 42, 40, 37],
       [47, 43, 39, 32]])
>>> data >= 40    # ndarrayをスカラと比較するとbool値の配列を返してくる
array([[ True,  True,  True, False],
       [ True,  True, False, False]])
>>> np.sum(data >= 40, axis=1)    # True(=1)の数を行ごとにカウント
array([3, 2])

目標達成確率-MICのグラフを作成する

次はグラフの作成部分です。

    # resultをプロットする
    fig, ax = plt.subplots()
    # 余白の調整(左と下を少し広げる). デフォルト値はleft=0.125, bottom=0.1
    fig.subplots_adjust(left=0.15, bottom=0.15)
 
    for dosage, (prob, marker) in result.items():
        plt.plot(mic_list, prob, label=dosage, marker=marker, markersize=6)
 
    # 達成確率80%に水平ラインを引く
    plt.hlines(y=80, xmin=0, xmax=mic_list[-1], alpha=0.5, color='red', linestyles='dashed')
    plt.xlim(mic_list[0], )
    plt.xscale('log', basex=2)
    plt.xticks(mic_list, rotation=30, labels=mic_list)
    plt.ylim(0, )
    plt.yticks(np.arange(0, 110, 10))
    plt.ylabel('Probability of 40%T>MIC\n  attainment(%)')
    plt.xlabel('MIC(μg/mL)')
    plt.grid(linestyle='dashed', linewidth=0.5)
    # テキスト(患者タイプ, CCr, Wt)を右上に表示
    plt.text(0.68, 0.93, pt_type, size=12,  horizontalalignment='left',verticalalignment='top', transform=ax.transAxes)
    # グラフラベル(dosage)を左下に表示
    plt.legend(loc='lower left', frameon=False)
    plt.show()

特に難しいところはなく、resultをfor文で回して投与方法（dosage）毎に各MICにおける%T>MIC（４０以上）値をプロットしていきます。横軸（MIC）はplt.xscale()の引数に'log'、対数の底（basex）2を指定して対数スケールにしています。

４種類の患者タイプでプログラムを走らせたところ、以下のグラフが得られました。

f:id:enokisaute:20200307151441p:plain
f:id:enokisaute:20200307151605p:plain

PK/PDブレイクポイントを確認する

最後は、目標達成確率80%を達成できる最大のMIC値を『個別化 PK-PD ブレイクポイント』として表にします。といっても、コードを実行するときの患者タイプ１種類のPK/PDブレイクポイントをDataFrameで表示するだけです。

    import pandas as pd
    # key: dosage, value: 80%を達成できる最大のMIC値
    table = {}
 
    for dosage, (prob, marker) in result.items():
        plt.plot(mic_list, prob, label=dosage, marker=marker, markersize=6)
        table[dosage] = str(mic_list[np.argmax(np.where(prob >= 80))])
 
    # 辞書からDataFrameを作成. keyが行であるため、orientには'index'を渡す
    df = pd.DataFrame.from_dict(table, orient='index')
    df.columns = [pt_type[:9]]
    print('PK/PD breakpoint(μg/mL)')
    print(df)

先ほどのグラフ作成のforループの中で辞書を作成しています。key, valueはコメントの通りです。
np.where()で条件を指定し、得られた配列の中で最大のインデックスをnp.argmax()で取得します。ここで、np.argmax()の引数は「prob>=80」では駄目です。上でもやりましたが、これは0か1のbool値の配列を返すため、これではうまくいきません。最初ちょっとつまづいたのでメモ。
あとはDataFrame（２次元配列のようなもの）を辞書から作成しています。
Patient Aの場合では、このコードによって次のような出力が得られます。

PK/PD breakpoint(μg/mL)
                  Patient A
0.25g twice           0.125
0.25g 3 times           0.5
0.5g twice             0.25
0.5g 3 times            1.0
1g twice                0.5
1g 3 times(daily)       2.0

Patient A～Dの結果はこのようになりました。

メロペネム投与方法 (0.5時間点滴)	Patient A	Patient B	Patient C	Patient D
0.25g twice	0.125	0.25	0.5	1.0
0.25g 3 times	0.5	1.0	1.0	2.0
0.5g 3 twice	0.25	0.5	1.0	2.0
0.5g 3 times	1.0	2.0	2.0	4.0
1g 3 twice	0.5	1.0	2.0	4.0
1g 3 times	2.0	4.0	4.0	8.0

私のプログラムでも参考文献「薬理学的観点からみた PK-PD 理論とブレイクポイント」の表と同じ結果となりました。しかし、前回触れましたが私のプログラムの場合は蛋白結合率は考慮していない値となっています。
また、上の４つのグラフを見ると、参考文献のものよりもほんの僅かですが、値が全体的に低いような気がします。前回記事の「%T>MICの確認」や「PK/PDパラメータのCmax/MIC, AUC/MICを求める」の「TDM解析ソフトとの比較」でもそうでしたが、私のプログラムでは参考文献やソフトの値に比べて値が少し低く出ていたので、今回の結果にもそれが表れているのかもしれません。

プログラムの処理時間について

このようなシミュレーションプログラムというのは、手軽にいろいろな条件を試してみることができ、少しずつ条件を変えながら、何度も実行することになります。そのため、プログラムの処理時間はできるだけ短い方が好ましいでしょう。
処理時間の大部分を占めると思われる、一番外側のループの処理時間（sample_size=10000）を計測したところ、私のパソコン（Core i5 8600K）ではだいたい78秒（パソコンのスペックによっては数分かかるかも）でした。処理時間を計る部分のコードは次のように処理の終了時刻から開始時刻を引くというやり方で求めました。

    import time
    # ループの処理時間を計測
    start_time = time.time()
 
    for dosage, value in recipe.items():

        # サンプルサイズのループ処理(省略)
    
    end_time = time.time()
    processing_time = end_time - start_time
    print('processing_time(sec): ', processing_time)

次回はコードの改善によって、この処理時間がどれくらい短くなるかを見てみようと思います。
www.yakupro.info

最後にコード全文を載せておきます。

import numpy as np
import matplotlib.pyplot as plt
from collections import OrderedDict
import pandas as pd
import time
 
# 抗菌薬でモンテカルロシミュレーション3
 
class BloodConcTwoComp2:
    STEP = 0.01 # 経過時間の刻み幅
    MAGNI = 7  # 半減期の何倍で定常状態とするか
 
    def __init__(self, A, alpha, B, beta, T, intervals):
        self.A = A
        self.alpha = alpha
        self.B = B
        self.beta = beta
        self.T = T
        self.intervals = intervals
        self.idx = 0  # 現在の投与間隔のインデックス
        self.n = 1  # 何回目の投与かを表す
        self.ex_c = [0 for _ in range(int(intervals[0] / self.STEP))]  # 延長分濃度. 0で初期化
        t_half = np.log(2) / beta
        # 投与間隔の倍数のうち、半減期のMAGNI倍以上の最小値を定常状態到達時間とする
        self.css_t = (t_half * self.MAGNI // intervals[0] + 1) * intervals[0]
        # 定常状態までの繰り返し投与回数. 最低3回
        cycle_to_css = 3 if self.css_t // intervals[0] < 3 else int(self.css_t / intervals[0])
        c = []
        for _ in range(cycle_to_css):
            c.extend(self.calc_section_conc())
        self.last_interval = c[-int(self.intervals[0] / self.STEP):]
 
    # 点滴開始～終了直前
    def calc_during_div(self):
        td = np.linspace(0, self.T, num=int(self.T / self.STEP), endpoint=False)
        return self.A * (1 - np.exp(-self.alpha * td)) + self.B * (1 - np.exp(-self.beta * td))
 
    # 点滴終了時～次回投与直前
    def calc_after_div(self, begin, end):
        section = np.linspace(begin, end, num=int((end - begin) / self.STEP), endpoint=False)
        return self.A * (np.exp(-self.alpha * (section - self.T)) - np.exp(-self.alpha * section)) \
               + self.B * (np.exp(-self.beta * (section - self.T)) - np.exp(-self.beta * section))
 
    # 点滴開始から次回投与までの血中濃度を返す.
    # 今回投与分と前回投与時までの延長分(今回投与間隔分)の濃度の合計を足し合わせて返す.
    def calc_section_conc(self):
        c1 = self.calc_during_div()
        c2 = self.calc_after_div(self.T, self.intervals[self.idx])
        c = np.concatenate([c1, c2], axis=0)  # c1とc2を最初の軸(0次元目)で連結
        sum_c = c + self.calc_extended_conc(self.n, self.ex_c)  # 単回投与＋延長分
        # 次回投与時のために以下の変数を更新
        self.n += 1
        self.idx = (self.idx + 1) % len(self.intervals)
        return sum_c.tolist()
 
    # 1回目~前回までの全ての投与回の延長(残存)分の血中濃度の合計を返す.
    # n: 投与回数, ex_c: 延長分の濃度の合計, begin: 求めるべき延長分濃度の計算開始位置
    def calc_extended_conc(self, n, ex_c, begin=0):
        if n == 1:
            return ex_c
        else:
            prev_idx = (n - 2) % len(self.intervals)  # 1つ前の投与間隔のインデックス
            prev_interval = self.intervals[prev_idx]
            begin += prev_interval
            end = begin + self.intervals[self.idx]
            ex_c += self.calc_after_div(begin, end)
            return self.calc_extended_conc(n - 1, ex_c, begin)
 
    # %T>MICを計算する
    def get_time_above_mic(self, mic):
        above_mic = len([c for c in self.last_interval if c > mic])
        tau = len(self.last_interval)
        return above_mic / tau * 100
 
 
if __name__ == '__main__':
    sample_size = 10000
 
    CCr = 100       # クレアチニンクリアランス(mL/min)
    Wt = 40         # 体重(kg)
    # グラフに貼るテキスト(患者タイプ, CCr, BWを表示)の作成
    pt_type = 'Patient A:\n'+'  CCr=' + str(CCr) + 'mL/min,\n  BW=' + str(Wt) + 'kg'
 
    CL = 0.0905 * CCr + 2.03
    Vc = 0.199 * Wt
    T = 0.5         # 点滴時間
 
    # 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)
 
    # 血中濃度推移に対して試すMIC(μg/mL)のリスト
    mic_list = np.array([0.063, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128])
 
    # 順序が保持される辞書を使用
    recipe = OrderedDict()
    # key: 投与方法
    # value: (投与量(mg), 1日の投与回数, グラフ表示マーカー(o:丸, ^:三角, s:四角))
    recipe = {'0.25g twice': (250, 2, 'o'),
              '0.25g 3 times': (250, 3, 'o'),
              '0.5g twice': (500, 2, '^'),
              '0.5g 3 times': (500, 3, '^'),
              '1g twice': (1000, 2, 's'),
              '1g 3 times(daily)': (1000, 3, 's')}
 
    # 発生させたパラメータから血中濃度計算に必要なalpha, betaを求める
    ke = X_cl / X_vc    # CL / Vc
    k12 = X_q / X_vc    # Q / Vc
    k21 = X_q / X_vp    # Q / Vp
    s = k12 + k21 + ke  # alpha + beta
    p = k21 * ke        # alpha * beta
    alpha = (s + np.sqrt(s ** 2 - 4 * p)) / 2
    beta = s - alpha
 
    # key: 投与方法(dosage)
    # value: (%T>MICが40以上となった数 / sample_size * 100, グラフ表示マーカー)
    result = {}
 
    # ループの処理時間を計測
    start_time = time.time()
 
    for dosage, value in recipe.items():
        print(dosage)
        x0 = value[0]         # 投与量
        k0 = x0 / T           # 点滴速度(mg/hr) = 投与量 / 点滴時間
        intervals = (int(24 / value[1]), )  # 投与間隔(hr) = 24 / 投与回数
 
        # MICリストの要素数×サンプルサイズの行列. ここに%T>MICを格納していく
        data_matrix = np.zeros((len(mic_list), sample_size))
 
        # 発生させたパラメータから血中濃度計算に必要なA, Bを求める
        A = (k0 * (alpha - k21)) / (X_vc * alpha * (alpha - beta))
        B = (k0 * (k21 - beta)) / (X_vc * beta * (alpha - beta))
 
        for si in range(sample_size):
            bc = BloodConcTwoComp2(A[si], alpha[si], B[si], beta[si], T, intervals)
            for mi, mic in enumerate(mic_list):
                data_matrix[mi, si] = bc.get_time_above_mic(mic)
 
        result[dosage] = (np.sum(data_matrix >= 40, axis=1) / sample_size * 100, value[-1])
 
    end_time = time.time()
    processing_time = end_time - start_time
    print('processing_time(sec): ', processing_time)
 
    # resultをプロットする
    fig, ax = plt.subplots()
    # 余白の調整(左と下を少し広げる). デフォルト値はleft=0.125, bottom=0.1
    fig.subplots_adjust(left=0.15, bottom=0.15)
 
    # key: dosage, value: 80%を達成できる最大のMIC値
    table = {}
 
    for dosage, (prob, marker) in result.items():
        plt.plot(mic_list, prob, label=dosage, marker=marker, markersize=6)
        table[dosage] = str(mic_list[np.argmax(np.where(prob >= 80))])
 
    # 辞書からDataFrameを作成. keyが行であるため、orientには'index'を渡す
    df = pd.DataFrame.from_dict(table, orient='index')
    df.columns = [pt_type[:9]]
    print('PK/PD breakpoint(μg/mL)')
    print(df)
 
    # 達成確率80%に水平ラインを引く
    plt.hlines(y=80, xmin=0, xmax=mic_list[-1], alpha=0.5, color='red', linestyles='dashed')
    plt.xlim(mic_list[0], )
    plt.xscale('log', basex=2)
    plt.xticks(mic_list, rotation=30, labels=mic_list)
    plt.ylim(0, )
    plt.yticks(np.arange(0, 110, 10))
    plt.ylabel('Probability of 40%T>MIC\n  attainment(%)')
    plt.xlabel('MIC(μg/mL)')
    plt.grid(linestyle='dashed', linewidth=0.5)
    # テキスト(患者タイプ, CCr, Wt)を右上に表示
    plt.text(0.68, 0.93, pt_type, size=12,  horizontalalignment='left',verticalalignment='top', transform=ax.transAxes)
    # グラフラベル(dosage)を左下に表示
    plt.legend(loc='lower left', frameon=False)
    plt.show()