経口投与のフリップ・フロップをグラフで見る

Flip-Flop(フリップ・フロップ)

前回記事の最後に、消失相の傾きはkaとkeの大小関係によってどちらが反映されるかは変わるということを書きました。
通常の ke<kaという条件においては消失が律速となりkeが反映されますが、ke>kaのときはkaが反映されます。薬物は体内に吸収されるよりも早くは消失しないため、吸収される速度が律速となるためです。
これをFlip-Flop(フリップ・フロップ)といいます。
今回はこれをグラフで確認してみようということで、kaとkeの値をスライダーを使って変化させることのできるグラフを作ります。

各条件で用いる数式

今までもいくつかmatplotlibでインタラクティブなグラフを作りましたが、同じようにmatplotlib.widgetsを使います。変化するkeとkaの値の取得、血中濃度の値を再計算して更新、再描画するなどのコードを書く必要があります。血中濃度式は、

=Aとおくと、

（吸収過程）

（消失過程）

（経口投与モデルの血中濃度式）

でした。各条件（ke<ka、 ke>ka、 ke=ka）における血中濃度式をみてみます。

ke<kaの場合

これは前回やったケースです。いわゆる普通の薬で、消化管からの吸収が速やかな場合は消失相（β相）の傾きにはkeが反映されます。この場合は上の３つの式をそのままプロットします。

ke>kaの場合

吸収されるまでに時間のかかる徐放性製剤などの場合に起こることがあります。
上で書いたようにこの場合はkaが反映されます。上の式でAをとおいていましたが、この条件（ke>ka）の場合はAがマイナスになってしまいます。最終的にCは同じように求めることができますが、各指数項もプロットしたいので、

を次のように変形した形でプロットします。

（前の項の分母のka, keを入れ替えて、カッコ内の項の順番を入れ替えた）

ke=kaの場合

この場合はAの分母が０となり、上の式は定義されなくなってしまいます。
この場合の血中濃度を表す式は

という式になるとのこと。ここでは各指数項の値はNone（なし）とし、k=ka=keとしてCのみプロットすることとします。

ではこれらをコードに書いていきます。

class PoOneComp:
    def __init__(self, F, D, vd, t):
        self.F = F
        self.D = D
        self.vd = vd
        self.t = t
 

    def absorption_process(self, a, ka):
            return a * np.exp(-ka * self.t)
 
    def elimination_process(self, a, ke):
            return a * np.exp(-ke * self.t)
 
    # ka, keの大小関係に応じて消失過程, 吸収過程, 血中濃度を計算する
    def calc_conc(self, ka, ke):
        if ka == ke:
            k = ka
            c_ab = None
            c_eli = None
            c = (k * self.F * self.D / self.vd) * np.exp(-k * self.t) * self.t
        elif ka < ke:
            # flip-flop
            A = (ka * self.F * self.D) / (self.vd * (ke - ka))
            c_ab = self.absorption_process(A, ka)
            c_eli = self.elimination_process(A, ke)
            # cを求める式は前後が入れ替わる
            c = c_ab - c_eli
        else:
            A = (ka * self.F * self.D) / (self.vd * (ka - ke))
            c_ab = self.absorption_process(A, ka)
            c_eli = self.elimination_process(A, ke)
            c = c_eli - c_ab
        return c_ab, c_eli, c

変数ka, keはスクリプトファイルのあちこちで参照・変更されることになるためクラスを作成し、変化するka, keの値に応じて結果を返すメソッドを実装していきます。
今回のプログラムでは変更されることのないF, D, vd, t(経過時間)をコンストラクタの引数としています。
コメントにある通り、メソッドcalc_concはka, keを引数に取り、その大小関係に応じて消失過程, 吸収過程, 血中濃度の３つの値を返します。各条件における数式は、上で書いたとおりのものです。

TmaxとCmax

またグラフの変化に応じてTmaxを図中に点線で表示させるようにします。そのためにCmaxの値も必要です。「経口投与１-コンパートメントモデルのグラフを描く」でこれらを求める式を示しましたが、問題はke=kaの場合です。この場合でも、上の項「ke=kaの場合」でみた血中濃度Cを表す式を微分して０とおきます。
ここは代数計算（数式処理）を行うライブラリのSympyにやってもらいました。

>>> import sympy            # モジュールのインポート
>>> t = sympy.Symbol('t')   # 変数を定義していく
>>> F = sympy.Symbol('F')
>>> D = sympy.Symbol('D')
>>> Vd = sympy.Symbol('Vd')
>>> k = sympy.Symbol('k')
>>> expr = sympy.diff(F * D * k / Vd * sympy.exp(-k * t) * t, t)    # 微分
>>> expr
-D*F*k**2*t*exp(-k*t)/Vd + D*F*k*exp(-k*t)/Vd
>>> sympy.solve(expr, t)    # 方程式を解く
[1/k]

Anacondaをインストールした方は最初から入っているため、そのまま上のように書いて使うことができます。
t = sympy.Symbol('t')と書くことでtを変数として定義しています。F, D, Vd, kについても同様です。次に、sympy.diff()を使って微分を行い、結果をexprに代入しています。複数の変数を含む今回のような式では、第二引数に対象の変数（ここではｔ）を指定します。
最後に、exprが０となるｔを求めるのにsympy.solve()を使っています。ここでも、第二引数に変数（ｔ）を指定することで、ｔについて解くことができます。

結果、ka = keのときのTmaxは1/kとなり、メソッドは次のように書きました。

    # Cmax到達時間
    def t_max(self, ka, ke):
        if ka == ke:
            return 1 / ka
        return np.log(ka / ke) / (ka - ke)

Cmaxはここで求めたTmaxを元の（微分する前の）式に入れることで計算します。

グラフと文字の更新処理

このへんのコードは以前書いたものとほぼ同じですのでそちらも参考にしてください。
薬物動態パラメータを変化させた場合の血中濃度推移

ここではコールバック関数（マウスによるスライダーの操作が発生したときに行う処理）について見てみます。

def ke_update(slider_val):
    ke = slider_val
    ka = ka_slider.val
    c_ab, c_eli, c = po.calc_conc(ka, ke)
    graph1.set_ydata(c_ab)
    graph2.set_ydata(c_eli)
    graph3.set_ydata(c)
    tmax = po.t_max(ka, ke)
    cmax = po.c_max(ka, ke)
    graph4.set_data([tmax, tmax], [0, cmax])
    # textを更新
    tmax = round(po.t_max(ka, ke), 1)
    tmax_text.set_text('Tmax(hr): ' + str(tmax))
    state_text.set_text(PoOneComp.get_state(ka, ke))
    fig.canvas.draw_idle()

これはクラスのメソッドではありません。スライダーの現在の値をke, kaとして取得し、メソッドcalc_conc()に渡して消失過程、吸収過程、血中濃度とTmax、Cmaxを再計算→グラフの値を更新→再描画という流れになっています。
また、グラフ中に現在のkeとkaの大小関係を表示するための文字列を返すget_stateというスタティックメソッドを作っています。

    # kaとkeの大小を不等号で示した文字列を返す
    @ staticmethod
    def get_state(ka, ke):
        s = '=' if ka == ke else '>' if ka > ke else '<'
        return 'ka ' + s + ' ke'

三項演算子をネストにしているので少しわかりにくいかもしれませんが、ka, keの大小関係（またはイコール）に応じて'='、'>'、'<'をsに代入しています。
一番左の最初の条件式ka == keがTrueならば'='を返し、Falseならば'>'を返すのではなく、次の条件式であるka > keを評価します。
ちなみに、このメソッドは単純にkaとkeの関係を文字列として返すだけなので、クラスに入れる必要はありません。別に関数で良かったのですが、このクラスに関係する機能の一部ということでこのようにしました。

フリップフロップを確認する

では実際にスライダーを動かしてみていきます。
ka=0.5, ke=0.15からスタートです。この状態からkaを少しずつ減らしていきます。

kaを0.3まで減少させました。最初のグラフよりもCmaxは低くなり、Tmaxは3.4から4.6に延長しているのが見て取れますが、Cmax以降の濃度推移の形はまだ最初のものと似ています。ここからkaをさらに0.07まで変化させると下のグラフのようになりました。フリップフロップの状態です。

Cmaxはさらに低くなり、この濃度が長い時間持続しています。このように、分布容積とクリアランスが同じ（=keが同じ）であっても、徐放性製剤にするなどの工夫によってkaを小さくすると、濃度推移のグラフの形は大きく変化するということがわかります。そして、上で書いたようにこのフリップフロップの状態では緑の濃度曲線の末端部分はkaを反映することになります。

作ってから思ったんですが、keのスライダーはいらなかったかも。

参考

• 薬物動態解析入門　はじめての薬物速度論
• 4ステップ　臨床力UPエクササイズ４　TDM領域
• Lecture 7.3: Flip-flop kinetics(You Tube)
• 徹底解説 薬物動態の数学―微積分と対数、非線形

今回のコードです。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
import sympy
 
# kaを変化させてフリップフロップを確認する
 
class PoOneComp:
    def __init__(self, F, D, vd, t):
        self.F = F
        self.D = D
        self.vd = vd
        self.t = t
 
    def absorption_process(self, a, ka):
            return a * np.exp(-ka * self.t)
 
    def elimination_process(self, a, ke):
            return a * np.exp(-ke * self.t)
 
    # 各ケースに応じて消失過程, 吸収過程, 血中濃度を計算する
    def calc_conc(self, ka, ke):
        if ka == ke:
            k = ka
            c_ab = None
            c_eli = None
            c = (k * self.F * self.D / self.vd) * np.exp(-k * self.t) * self.t
        elif ka < ke:
            # flip-flop
            A = (ka * self.F * self.D) / (self.vd * (ke - ka))
            c_ab = self.absorption_process(A, ka)
            c_eli = self.elimination_process(A, ke)
            # cを求める式は前後が入れ替わる
            c = c_ab - c_eli
        else:
            A = (ka * self.F * self.D) / (self.vd * (ka - ke))
            c_ab = self.absorption_process(A, ka)
            c_eli = self.elimination_process(A, ke)
            c = c_eli - c_ab
        return c_ab, c_eli, c
 
    # Cmax到達時間
    def t_max(self, ka, ke):
        if ka == ke:
            return 1 / ka
        return np.log(ka / ke) / (ka - ke)
 
    # 最高血中濃度
    def c_max(self, ka, ke):
        if ka == ke:
            k = ka
            tmax = float(self.t_max(ka, ke))
            return (k * self.F * self.D / self.vd) * np.exp(-k * tmax) * tmax
        return self.F * self.D / self.vd * (ka / ke) ** (ke / (ke - ka))
 
    # kaとkeの大小を不等号で示した文字列を返す
    @ staticmethod
    def get_state(ka, ke):
        s = '=' if ka == ke else '>' if ka > ke else '<'
        return 'ka ' + s + ' ke'
 
 
def ke_update(slider_val):
    ke = slider_val
    ka = ka_slider.val
    c_ab, c_eli, c = po.calc_conc(ka, ke)
    graph1.set_ydata(c_ab)
    graph2.set_ydata(c_eli)
    graph3.set_ydata(c)
    tmax = po.t_max(ka, ke)
    cmax = po.c_max(ka, ke)
    graph4.set_data([tmax, tmax], [0, cmax])
    tmax = round(po.t_max(ka, ke), 1)
    tmax_text.set_text('Tmax: ' + str(tmax))
    state_text.set_text(PoOneComp.get_state(ka, ke))
    fig.canvas.draw_idle()
 
 
def ka_update(slider_val):
    ka = slider_val
    ke = ke_slider.val
    c_ab, c_eli, c = po.calc_conc(ka, ke)
    graph1.set_ydata(c_ab)
    graph2.set_ydata(c_eli)
    graph3.set_ydata(c)
    tmax = po.t_max(ka, ke)
    cmax = po.c_max(ka, ke)
    graph4.set_data([tmax, tmax], [0, cmax])
    tmax = round(po.t_max(ka, ke), 1)
    tmax_text.set_text('Tmax: ' + str(tmax))
    state_text.set_text(PoOneComp.get_state(ka, ke))
    fig.canvas.draw_idle()
 
 
F = 1 
D = 100
vd = 10
Ka = 0.5
Ke = 0.15
 
A = (Ka * F * D) / (vd * (Ka - Ke))
time = np.arange(0, 50, 0.1)
po = PoOneComp(F, D, vd, time)
 
c_ab = po.absorption_process(A, Ka)
c_eli = po.elimination_process(A, Ke)
t_max = po.t_max(Ka, Ke)
c_max = po.c_max(Ka, Ke)
 
fig, ax = plt.subplots()
plt.subplots_adjust(bottom=0.3, top=0.9)
graph1, = plt.plot(time, c_ab, color='blue', alpha=0.8, label='c2:吸収過程')
graph2, = plt.plot(time, c_eli, color='red', alpha=0.8, label='c1:消失過程')
graph3, = plt.plot(time, c_eli - c_ab, color='green', alpha=0.8, label='c:血中濃度', linewidth=2)
# tmaxに補助線を入れる
graph4, = plt.plot([t_max, t_max], [0, c_max], color='crimson', linestyle='dotted')
# 軸ラベルを表示
plt.ylabel('conc[μg/mL]')
plt.xlabel('time[hr]')
plt.xlim(0, 40)
plt.ylim(0, 20)
plt.legend()
 
# fig内でのaxes(グラフ)座標を取得
ax_pos = ax.get_position()
# axes([left, under, width, height])
ke_slider_pos = plt.axes([ax_pos.x0, ax_pos.y0 - 0.15, ax_pos.width, 0.03])
ka_slider_pos = plt.axes([ax_pos.x0, ax_pos.y0 - 0.23, ax_pos.width, 0.03])
 
ke_slider = Slider(ke_slider_pos, 'ke', 0.01, 1.0, valinit=Ke, valstep=0.01, valfmt="%0.3f")
ka_slider = Slider(ka_slider_pos, 'ka', 0.01, 0.5, valinit=Ka, valstep=0.01, valfmt="%0.3f")
ke_slider.on_changed(ke_update)
ka_slider.on_changed(ka_update)
 
# 図に文字列を挿入(横位置, 縦位置, 文字列). textオブジェクトへの参照を持っておく
state_text = fig.text(ax_pos.width - 0.05, ax_pos.y0 + 0.36, PoOneComp.get_state(Ka, Ke), fontsize=18)
tmax_text = fig.text(ax_pos.width - 0.06, ax_pos.y0 + 0.26, "Tmax: " + str(round(po.t_max(Ka, Ke), 1)),
                     fontsize=15)
plt.show()