今回は２-コンパートメント点滴静注モデルの繰り返し投与のプログラムを書いてみたいと思います。
2-コンパートメントモデルの繰り返し投与
1-コンパートメント点滴静注のようにｎ回目の血中濃度を表す一般化された式を作れば簡単ですが、あまり面白みがないのでこの方法はやりません。
前記事の「バンコマイシンの血中濃度データに2-コンパートメントモデルを当てはめる」の(2)の式
を使い、1-コンパートメント点滴静注モデルでやったやり方（前回投与の延長分を今回投与分の血中濃度と足し合わせる方法）ができるかもと思ったのですが、どうやらこの式からだけでは前回投与の延長分の血中濃度を求めることができず、下図の赤線部分のようになってしまいました。
そこで、別の実装方法を考えてみます。
www.yakupro.info
で見たように、２-コンパートメントモデルといっても投与終了後、充分に時間が経過していれば血中濃度は次のような近似式で表すことができるのでした。
前回投与の延長分の血中濃度は、投与終了後から充分に時間が経過しているとみなせば、この式を使って計算できそうです。ということで、この式を使った実装を書いてみたいと思います。
投与終了何時間後から近似するか
基本的な考え方は１-コンパートメントモデルの「1-コンパートメント点滴静注ー不均等の間隔・投与量に対応したプログラムを作るー」でやったときと同じです。
一次速度で減少していく血中濃度をで計算していましたが、この
は点滴終了時の血中濃度を使いました。今回はこの
を点滴終了してから充分に時間が経ってからの血中濃度に置き換えることになります。
プログラムを実装する上では、”充分に時間が経ってから”とは具体的にどれくらいか、を決めなくてはいけませんが、ここでは、取りあえず”点滴終了してから10時間後”としてみます。それだと、じゃあ１日３回投与（投与間隔8時間）とかだとどうなるの？ってまたなってくるのですが、今はあまり細かいことは考えずにこれで話を進めていきます。
近似を使って血中濃度を求めるクラスを書く
上で説明したように、延長分の血中濃度の計算部分は以前やったものと同じなので、大きく変わりません。まずはコンストラクタから。
class BloodConcTwoComp1: STEP = 0.01 # 経過時間の刻み幅 def __init__(self, A, alpha, B, beta, T, intervals): self.A = A self.alpha = alpha self.B = B self.beta = beta self.idx = 0 self.T = T self.intervals = intervals # タプルで受け取る. self.c_peak = 0 self.peak_t = 10 # 点滴終了してこの時間以降は近似式で計算する. (投与間隔 - 点滴時間)以下で設定する. self.extended_c = [0 for _ in range(int(intervals[0] / self.STEP))] # 0で初期化
インスタンス生成時にパラメータとしてA、α、B、β、点滴時間、投与間隔のタプルを受け取ります。今回はインスタンスを生成する前にCCr（クレアチニンクリアランス）と体重からA、α、B、βを計算*1しています。本来はCCrと体重をクラスに渡し、クラス内でこれらのパラメータを計算するべきだと思いますが、今回はクラスが複雑にならないようこのような実装としました。
self.c_peakが以前の実装のに当たります。変数の名前をc_peakとしていますが、これはいわゆる分布が終了した時点の「ピーク濃度」のことではありませんのでご注意を。
その次の行のself.peak_t = 10についてです。今回は決め打ちで10（時間）と入れていますが、点滴終了してこの時間が経過した時点の濃度をc_peakとし、これ以降の濃度の計算をでするということになります。
次に、単回投与の血中濃度を計算するメソッドですが、点滴中と点滴終了後の計算式は以下のものをそのまま使っています。
＜点滴中＞
＜点滴終了後＞
# 点滴開始～終了直前 def __calc_during_div(self): td = np.linspace(0, self.T, num=int(self.T / self.STEP), endpoint=False) c = self.A * (1 - np.exp(-self.alpha * td)) + self.B * (1 - np.exp(-self.beta * td)) return c.tolist() # 点滴終了時～次回投与直前 def __calc_after_div(self): ta = np.linspace(self.T, self.intervals[self.idx], num=int((self.intervals[self.idx] - self.T) / self.STEP), endpoint=False) c = self.A * (np.exp(-self.alpha * (ta - self.T)) - np.exp(-self.alpha * ta)) \ + self.B * (np.exp(-self.beta * (ta - self.T)) - np.exp(-self.beta * ta)) return c.tolist()
延長分の血中濃度と単回投与の血中濃度を足し合わせるメソッドは次のようになります。ロジックとしては「1-コンパートメント点滴静注ー不均等の間隔・投与量に対応したプログラムを作るー」でやったときとまったく同じですが、点滴終了してpeak_t（１０時間）後の濃度をc_peakとして保存している部分が異なります。
# 点滴開始から次回投与までのcを計算する. # 前回投与時の延長分と今回投与分の血中濃度を合計して返す. def calc_section_conc(self): c1 = self.__calc_during_div() c2 = self.__calc_after_div() c1.extend(c2) sum_c = [x + y for (x, y) in zip(c1, self.extended_c)] # c_peakは点滴終了してpeak_t時間後とする if self.T + self.peak_t == self.intervals[self.idx]: peak_idx = int((self.T + self.peak_t) / self.STEP) - 1 else: peak_idx = int((self.T + self.peak_t) / self.STEP) self.c_peak = sum_c[peak_idx] # 次回投与のために延長分を計算しておく self.extended_c = self.__calc_extended_conc() return sum_c
また、次回投与間隔分の血中濃度計算メソッド（__calc_extended_conc()）も以下の式以外は同じです。
# c_peakからの減少は近似式で計算
extended_c = self.c_peak * np.exp(-self.beta * ex_t)
グラフの形を確認してみる
ではこれでグラフを描いてみます。
赤いラインが延長分の血中濃度で、これに単回投与の血中濃度を足してグラフを描画しています。パッと見うまく描けているように見えますが、投与間隔を6時間、peak_tを5時間に短くして描画してみると次のようになります。
赤いラインの延長分の血中濃度がだいぶ不自然な形になっています。これは1-コンパートメントの式で近似するには５時間では短かったためで、先の投与間隔12時間の場合と比べると、より誤差が大きくなっています。また、A、B、α、βのパラメータによってはもっと不自然になる場合もあると考えられます。
というわけで、今回のプログラムは繰り返し投与の描画には使えなさそうだということがわかりました。
次回は別の実装を考えてみたいと思います。
www.yakupro.info
最後に、今回の全コードを載せておきます。
import numpy as np import matplotlib.pyplot as plt import math # 2コンパートメントモデル繰り返し点滴静注の薬物血中濃度推移 # 1次速度による近似バージョン class BloodConcTwoComp1: STEP = 0.01 # 経過時間の刻み幅 def __init__(self, A, alpha, B, beta, T, intervals): self.A = A self.alpha = alpha self.B = B self.beta = beta self.idx = 0 self.T = T self.intervals = intervals # タプルで受け取る. self.c_peak = 0 self.peak_t = 10 # 点滴終了してこの時間以降は近似式で計算する. (投与間隔 - 点滴時間)以下で設定する. self.extended_c = [0 for _ in range(int(intervals[0] / self.STEP))] # 0で初期化 # 延長分の血中濃度のグラフ表示用. 区間begin～endを表す self.b = intervals[0] self.e = intervals[0] * 2 # 点滴開始～終了直前 def __calc_during_div(self): td = np.linspace(0, self.T, num=int(self.T / self.STEP), endpoint=False) c = self.A * (1 - np.exp(-self.alpha * td)) + self.B * (1 - np.exp(-self.beta * td)) return c.tolist() # 点滴終了時～次回投与直前 def __calc_after_div(self): ta = np.linspace(self.T, self.intervals[self.idx], num=int((self.intervals[self.idx] - self.T) / self.STEP), endpoint=False) c = self.A * (np.exp(-self.alpha * (ta - self.T)) - np.exp(-self.alpha * ta)) \ + self.B * (np.exp(-self.beta * (ta - self.T)) - np.exp(-self.beta * ta)) return c.tolist() # 点滴開始から次回投与までのcを計算する. # 前回投与時の延長分と今回投与分の血中濃度を合計して返す. def calc_section_conc(self): c1 = self.__calc_during_div() c2 = self.__calc_after_div() c1.extend(c2) sum_c = [x + y for (x, y) in zip(c1, self.extended_c)] # c_peakは点滴終了してpeak_t時間後とする if self.T + self.peak_t == self.intervals[self.idx]: peak_idx = int((self.T + self.peak_t) / self.STEP) - 1 else: peak_idx = int((self.T + self.peak_t) / self.STEP) self.c_peak = sum_c[peak_idx] # 次回投与のために延長分を計算しておく self.extended_c = self.__calc_extended_conc() return sum_c # 次回投与間隔分のcを計算する. def __calc_extended_conc(self): n_idx = (self.idx + 1) % len(self.intervals) # 次回投与間隔のインデックス current_interval = self.intervals[self.idx] # 現在の投与間隔 next_interval = self.intervals[n_idx] # 次回の投与間隔 begin = current_interval - self.T - self.peak_t end = begin + next_interval ex_t = np.linspace(begin, end, num=int((end - begin) / self.STEP), endpoint=False) # c_peakからの減少は近似式で計算 extended_c = self.c_peak * np.exp(-self.beta * ex_t) self.idx = n_idx ############################################ # 延長分の血中濃度を赤いラインで表示 t = np.linspace(self.b, self.e, num=int((self.e - self.b) / self.STEP), endpoint=False) plt.plot(t, extended_c, 'r-') self.b += intervals[0] self.e += intervals[0] ############################################ return extended_c if __name__ == '__main__': CCr = 50 # クレアチニンクリアランス(mL/min) Wt = 50 # 体重(kg) x0 = 500 # 投与量(mg) T = 1 # 点滴時間(hr) intervals = (12, ) # 投与間隔 cycle = 12 # intervalsを何サイクル繰り返すか graphlabel = 'VCM500mg×12時間毎' # バンコマイシンTDMソフト「MEEK」の薬物動態パラメータを使用 if CCr >= 85: CL = 3.83 else: CL = 0.537 * CCr * 60 / 1000 + 0.32 # L/h Vc = 2.32 * 0.206 * Wt # 中央コンパートメントの分布容積(L) Vt = 60.5578 # 抹消コンパートメントの分布容積(L) Q = 8.81 # コンパートメント間における血流速度(L/h) k12 = Q / Vc # 中央→抹消コンパートメントへの消失速度定数(hr^-1) k21 = Q / Vt # 抹消→中央コンパートメントへの消失速度定数(hr^-1) ke = CL / Vc X = k12 + k21 + ke # alpha + beta Y = k21 * ke # alpha * beta alpha = (X + math.sqrt(X ** 2 - 4 * Y)) / 2 beta = X - alpha A = (x0 / T * (k21 - alpha)) / (Vc * alpha * (beta - alpha)) B = (x0 / T * (k21 - beta)) / (Vc * (alpha - beta) * beta) print(A, B) print(alpha, beta) # パラメータを渡してインスタンスを生成 bc = BloodConcTwoComp1(A, alpha, B, beta, T, intervals) conc1 = [] bc = {'blue': (bc, conc1, cycle, graphlabel)} for color, value in bc.items(): _bc, conc, cycle, g_label = value for _ in range(len(_bc.intervals) * cycle): conc.extend(_bc.calc_section_conc()) period = sum(_bc.intervals) * cycle # グラフの描画区間 time = np.linspace(0, period, num=int(period / BloodConcTwoComp1.STEP), endpoint=False) plt.plot(time, conc, color=color, alpha=0.8, label=g_label) plt.ylabel('血中濃度[μg/mL]') plt.xlabel('時間[hr]') plt.xticks(np.arange(0, time[-1] + 24, 24)) plt.xlim(0,) plt.ylim(0,) plt.legend() plt.grid(linestyle='dashed') plt.show()
*1:薬物動態パラメータはバンコマイシンTDM解析ソフト「MEEK」で使われているものと同じものを使用しました
