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

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

経口繰り返し投与のグラフを描く

n回繰り返し投与後の血中濃度

投与間隔τ(タウ)でn回繰り返し投与後のt時間経過時の血中濃度C_nは次の式のようになります。

\large{C_n = \frac{F \cdot D \cdot ka}{Vd(ka - ke)}\{e^{-k_e \cdot t}(\frac{1-e^{-n \cdot k_e \cdot τ}}{1-e^{-k_e \cdot τ}})-e^{-k_a \cdot t}(\frac{1-e^{-n \cdot k_a \cdot τ}}{1-e^{-k_a \cdot τ}})\}}\ \ \ \ \ (1)


以前は静注モデル「1-コンパートメント点滴静注繰り返し投与のグラフ描画」でやりましたが、そのときと同じように考えてこの式を導くことができます。

経口投与の定常状態における最高・最低血中濃度とAUC

定常状態の血中濃度C_{ss}についてですが、n=∞としたときはe^{nkτ}≒0となり、(1)式は以下のようになります。

\large{C_{ss} = \frac{F \cdot D \cdot ka}{Vd(ka - ke)}\{e^{-k_e \cdot t}(\frac{1}{1-e^{-k_e \cdot τ}})-e^{-k_a \cdot t}(\frac{1}{1-e^{-k_a \cdot τ}})\}}\ \ \ \ \ (2)


この(2)式において、\large{t=τ}のときが定常状態における最低血中濃度となります。

次に、Tmaxについて考えてみます。これも考え方は単回投与のときに「経口投与1-コンパートメントモデルのグラフを描く」でやったのと同じです。(2)式を時間tで微分して、その値が0となる時間で求められます。


\large{T_{\rm ss\_max}=\frac{ln\{ka(1-e^{-keτ})/ke(1-e^{-kaτ})\}}{ka - ke}}


また、Cmaxはこの求めたTmaxを(2)式に代入すれば求められます。

次に、定常状態における平均血中濃度(Css_ave)です。1-コンパートメント点滴静注モデルでもみましたが、これは定常状態における1投与間隔のAUCを投与間隔で割ったものとなります。

\large{C_{\rm ss\_ave} =\frac{1}{τ}・\int_0^τ C_{ss}dt=\frac{1}{τ} \cdot \frac{F \cdot D}{k_e \cdot Vd}}


Sympyによる微分と定積分

上でTmaxとAUCを求めるのに(2)式の微分と定積分の計算結果をさらっと書きましたが、ここではSympyを使ってこれらの式を求めてみます。
前回も少し触れましたが、Sympyは代数計算(数式処理)を行うためのライブラリで、方程式を解いたり、式の展開や因数分解、微分積分などを計算することができます。Anacondaをインストールした人は最初から入っているので、そのまま下のように書いて使うことができます。

import sympy
 
# 変数を定義. まとめてすることもできる
t = sympy.Symbol('t')
ka = sympy.Symbol('ka')
F = sympy.Symbol('F')
D = sympy.Symbol('D')
Vd = sympy.Symbol('Vd')
ke = sympy.Symbol('ke')
tau = sympy.Symbol('tau')
 
# 上の変数を使って定常状態の血中濃度式(2)を定義する
f = ka * F * D / Vd / (ka - ke) * \
        (sympy.exp(-ke * t) / (1-sympy.exp(-ke * tau))
        - sympy.exp(-ka * t) / (1-sympy.exp(-ka * tau)))
 
# 微分. 第二引数に対象の変数を指定する
d_f = sympy.diff(f, t)
 
# tについて解く
print(sympy.solve(d_f, t))
# [-log(ke*(exp(ka*tau) - 1)*exp(-tau*(ka - ke))/(ka*(exp(ke*tau) - 1)))/(ka - ke)]
 
# 定積分. (関数, (積分する変数, 始点, 終点))
s_f = sympy.integrate(f, (t, 0, tau))
 
# 数式を簡略化
print(sympy.simplify(s_f))
# D*F/(Vd*ke)

まず変数と数式を定義しています。対象の変数を指定してsympy.diff()で微分をして、sympy.solve()でtについて解いています。結果は#に書いた式となりましたが、 整理すると上に書いたTmaxの式となります。
次に定積分ですが、sympy.integrate()で行います。s_fのままだと長い複雑な式ですが、数式を簡略化するsympy.simplify()を使うとD*F/(Vd*ke)となります。ただし、このsimplify()は自動で式を簡略化してくれるのですが、思うような結果にならないこともあります。

経口繰り返し投与のクラス

ではクラスを作成します。n回目投与時の血中濃度、Css_max、Css_min、Css_ave、Tss_max、AUCを計算するメソッドを持たせます。

class PoOneComp2:
    STEP = 0.01     # 経過時間の刻み幅
 
    def __init__(self, F, D, vd, ka, ke, interval):
        self.F = F
        self.D = D
        self.vd = vd
        self.ka = ka
        self.ke = ke
        self.interval = interval
        self.A = F * D * ka / vd / (ka - ke)
        self.t = np.linspace(0, interval, num=int(interval / PoOneComp2.STEP), endpoint=False)
 
    # n回目の1投与区間の血中濃度
    def get_section_conc(self, n):
        c = self.A * np.exp(-self.ke * self.t) * (1 - np.exp(-n * self.ke * self.interval)) \
            / (1 - np.exp(-self.ke * self.interval)) \
            - self.A * np.exp(-self.ka * self.t) * (1 - np.exp(-n * self.ka * self.interval)) \
            / (1 - np.exp(-self.ka * self.interval))
        return c.tolist()
 
    # 定常状態における最低血中濃度
    def get_cmin(self):
        cmin = self.A * np.exp(-self.ke * self.interval) / (1 - np.exp(-self.ke * self.interval)) \
            - self.A * np.exp(-self.ka * self.interval) / (1 - np.exp(-self.ka * self.interval))
        return cmin
 
    # 定常状態における最高血中濃度
    def get_cmax(self):
        tmax = self.get_tmax()
        cmax = self.A * np.exp(-self.ke * tmax) / (1 - np.exp(-self.ke * self.interval)) \
               - self.A * np.exp(-self.ka * tmax) / (1 - np.exp(-self.ka * self.interval))
        return cmax
 
    # 定常状態における最高血中濃度到達時間
    def get_tmax(self):
        tmax = np.log(self.ka / self.ke * (1 - np.exp(-self.ke * self.interval))
                      / (1 - np.exp(-self.ka * self.interval))) / (self.ka - self.ke)
        return tmax
 
    # AUC(投与間隔)
    def get_auc(self):
        return self.F * self.D / self.ke / self.vd
 
    # 定常状態における平均血中濃度
    def get_cssave(self):
        return self.get_auc() / self.interval

n回目投与時の血中濃度は1つの数式で書き表せるため、静注モデルのときと比べるとコードの複雑さは幾分マシだと思います。数式中の変数tについては、コンストラクタで0~\large{τ}までの時間をnp.linspace()を使って生成したものを使っています。
その他のメソッドについても、上に書いた数式をそのままタイプしただけです。

血中濃度推移の描画

今回各パラメータはテオフィリンのインタビューフォームの値を用いています。描画する際はこれらの値でインスタンスを生成し、投与日数から描画区間の投与回数を求めます。次にループを回してget_section_conc(i)のiに投与回数を1、2、3、・・・と渡し、その結果のリストを繋ぎ合わせていくことで血中濃度を生成しています。
このあたりのコードは1-コンパートメント点滴静注繰り返し投与ー最高・最低血中濃度とAUCーのグラフ描画部分とほぼ同じです。

f:id:enokisaute:20200314215743p:plain

color:  blue
投与量: 200, 1日投与回数: 2
Css_min:  5.181729261402002
Css_max:  7.873290748913346
Tss_max:  3.932301987630879
AUC(投与間隔):  82.30452674897118
Css_ave:  6.858710562414266
------------------------------
color:  red
投与量: 400, 1日投与回数: 1
Css_min:  2.7798919395863186
Css_max:  10.447306275126877
Tss_max:  5.242428650135141
AUC(投与間隔):  164.60905349794237
Css_ave:  6.858710562414266
------------------------------



今回のコードです。

import numpy as np
import matplotlib.pyplot as plt
 
# 経口1-コンパートメントモデル繰り返し投与
 
class PoOneComp2:
    STEP = 0.01     # 経過時間の刻み幅
 
    def __init__(self, F, D, vd, ka, ke, interval):
        self.F = F
        self.D = D
        self.vd = vd
        self.ka = ka
        self.ke = ke
        self.interval = interval
        self.A = F * D * ka / vd / (ka - ke)
        self.t = np.linspace(0, interval, num=int(interval / PoOneComp2.STEP), endpoint=False)
 
    # n回目の1投与区間の血中濃度
    def get_section_conc(self, n):
        c = self.A * np.exp(-self.ke * self.t) * (1 - np.exp(-n * self.ke * self.interval)) \
            / (1 - np.exp(-self.ke * self.interval)) \
            - self.A * np.exp(-self.ka * self.t) * (1 - np.exp(-n * self.ka * self.interval)) \
            / (1 - np.exp(-self.ka * self.interval))
        return c.tolist()
 
    # 定常状態における最低血中濃度
    def get_cmin(self):
        cmin = self.A * np.exp(-self.ke * self.interval) / (1 - np.exp(-self.ke * self.interval)) \
            - self.A * np.exp(-self.ka * self.interval) / (1 - np.exp(-self.ka * self.interval))
        return cmin
 
    # 定常状態における最高血中濃度
    def get_cmax(self):
        tmax = self.get_tmax()
        cmax = self.A * np.exp(-self.ke * tmax) / (1 - np.exp(-self.ke * self.interval)) \
               - self.A * np.exp(-self.ka * tmax) / (1 - np.exp(-self.ka * self.interval))
        return cmax
 
    # 定常状態における最高血中濃度到達時間
    def get_tmax(self):
        tmax = np.log(self.ka / self.ke * (1 - np.exp(-self.ke * self.interval))
                      / (1 - np.exp(-self.ka * self.interval))) / (self.ka - self.ke)
        return tmax
 
    # AUC(投与間隔)
    def get_auc(self):
        return self.F * self.D / self.ke / self.vd
 
    # 定常状態における平均血中濃度
    def get_cssave(self):
        return self.get_auc() / self.interval
 
 
if __name__ == '__main__':
    # 薬物動態パラメータはテオフィリンのインタビューフォームより
    F = 1
    vd = 0.45 * 60
    ka = 0.29
    ke = 0.09
    days = 4  # 何日投与するか
 
    # 1回の投与量(mg), 投与間隔(hr)を変えて投与する.
    D_1, interval_1 = 200, 12
    D_2, interval_2 = 400, 24
 
    # インスタンス生成
    po1 = PoOneComp2(F, D_1, vd, ka, ke, interval_1)
    po2 = PoOneComp2(F, D_2, vd, ka, ke, interval_2)
 
    # 描画区間の投与回数
    n1 = int(days * 24 / interval_1)
    n2 = int(days * 24 / interval_2)
 
    conc1, conc2 = [], []
 
    # key: 色, value: (インスタンス, 投与回数, 血中濃度)
    po = {'blue': (po1, n1, conc1), 'red': (po2, n2, conc2)}
 
    period = 24 * days
    time = np.linspace(0, period, int(period / PoOneComp2.STEP), endpoint=False)
 
    for color, value in po.items():
        _po, n, conc = value
        # 各クラス毎に1~n回目の投与までの血中濃度を加えていく
        for i in range(1, n + 1):
            conc.extend(_po.get_section_conc(i))
        # 1回投与量, 1日投与回数
        dose = int(_po.D)
        num = int(24 / _po.interval)
        plt.plot(time, conc, c=color, alpha=0.8, label=str(dose) + 'mg×' + str(num) + '回/日')
        print('color: ', color)
        print('投与量: {}, 1日投与回数: {}'.format(dose, num))
        print('Css_min: ', _po.get_cmin())
        print('Css_max: ', _po.get_cmax())
        print('Tss_max: ', _po.get_tmax())
        print('AUC(投与間隔): ', _po.get_auc())
        print('Css_ave: ', _po.get_cssave())
        print('-' * 30)
 
    plt.title('テオフィリンの血中濃度推移(イメージ)')
    # 軸ラベルを表示
    plt.ylabel('血中濃度[μg/mL]')
    plt.xlabel('時間[hr]')
    plt.xlim(0, )
    # 目盛りを24時間毎にする
    plt.xticks(np.arange(0, time[-1] + 24, 24))
    plt.ylim(0, 14)
    plt.grid(linestyle='dashed')
    plt.legend()
    plt.show()
 


参考文献

・徹底解説 薬物動態の数学ー微積分と対数、非線形ー