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

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

トップ > 薬物動態 > トラフ値からkeを推定する

トラフ値からkeを推定する

薬物動態

今回の記事は参考文献『４ステップ　臨床力UPエクササイズ４　TDM領域』の『Step3 薬物治療の個別化（薬物投与設計）の実践』の項目を基に書いています。
この項目の概要としては、喘息症状が悪化した患者のスクリーニングを通して薬物クリアランスの変動を考慮した投与設計を行っていく、という内容で、その中に喫煙時と非喫煙時の実測のトラフ値からkeを推定する、という部分があります。
参考文献の中では、Excelでkeを変化させながら定常状態におけるトラフ値が実測値に近づくようフィッティングしていますが、ここではPythonでコードを書いてやってみたいと思います。


喫煙がテオフィリンの血中濃度に与える影響

テオフィリンは主にCYP1A2で代謝されますが、喫煙によりこの酵素が誘導され、クリアランスが上昇することが知られています。そのため、喫煙時にはテオフィリン血中濃度が低下し、非喫煙時よりも高い用量が必要になることがあります。

定常状態におけるトラフ値

前回記事でも書きましたが、最低血中濃度（トラフ値）は定常状態での血中濃度\large{C_{\rm ss}}を表す式において、\large{t=τ}のときとなるため以下の式になります。

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

この\large{C_{\rm ss\_min}}がある特定の値のときのkeを求める、ということになります。

Sympyで方程式を代数的に解く

上で示した式において、ke以外のパラメータの値（ F, D, Vd, ka, τ）はすべて既知の値であるとします。ここで、非喫煙時のトラフ値が9μg/mL、喫煙時のトラフ値が
3μg/mLだったとして、それぞれのトラフ値におけるkeを求めます。

import sympy
 
# 変数を定義
ke = sympy.Symbol('ke')
 
F = 1                 # バイオアベイラビリティ
D = 200               # 投与量(mg)
vd = 0.46 * 63.6      # 分布容積(L/kg) * 理想体重(kg)
ka = 0.36             # 吸収速度定数(hr^-1)
tau = 12              # 投与間隔(hr)
trough = 9            # 実測トラフ値(μg/mL)
 

# 数式fが見にくくなるのでAとして先にまとめておく
A = F * D * ka / vd / (ka - ke)
 
# 数式を定義
f = A * sympy.exp(-ke * tau) / (1 - sympy.exp(-ke * tau)) \
    - A * sympy.exp(-ka * tau) / (1 - sympy.exp(-ka * tau)) - trough
 
# 方程式を解く
print(sympy.solve(f, ke))

各パラメータは参考文献のもの（これもインタビューフォームの値です）をそのまま使用しています。分布容積については、この問題で想定する患者は100kgですが、実測体重で薬物量を計算すると過量となるためBMIが22となる理想体重で計算しています。
トラフ値はまずは非喫煙時として9を設定しました。また、sympy.solve()では()内の式が＝0とした場合の解が得られるため、数式ｆではtroughを移項した形にしています。
このようにして、今までと同様にSympyのsolve()を使って解こうとしましたが、次のようなエラーメッセージが表示されうまくいきません。

No algorithms are implemented to solve equation

「方程式を解くためのアルゴリズムが実装されていない」と出ています。どうやら代数計算では解けないようなので、数値計算で解を求めることにしました。
そこで、調べたところsympyで非線形方程式を数値的に解く関数としてnsolve()というものがあり、これを使おうと思ったのですが、関数に渡す初期値がシビア（正解に近い値でないと解が得られなかった）だったため、結局今までも何度かお世話になっているscipy.optimizeを使用することにしました。

Css_minの式を可視化する

先にsympyで解きたい方程式ｆを可視化してグラフの概形を確認しておきます。
sympyではplottingモジュールを使うことで、定義した数式をそのまま渡すだけでプロットすることができます。上で書いたコードの下に、次のように書きます。
上で示したコードのprint(sympy.solve(f, ke))はエラーが出るので消しておいて下さい。

from sympy.plotting import plot
 
# keが[0.01, 0.15]の範囲でプロット
plot(f, (ke, 0.01, 0.15))

keの範囲はインタビューフォームの値から当たりをつけて上記で設定してみました。
f:id:enokisaute:20200315155345p:plain
これを見るとだいたい0.05を過ぎたあたりに解がありそうです。ちなみに、このグラフだけみると単純な形をしているように見えますが、keの範囲を広げてmatplotlibでプロットすると次のような形をしていました。
f:id:enokisaute:20200317215154p:plain

Scipyで方程式の数値解を求める

scipy.optimize モジュールのnewton()関数を使います。公式ドキュメントによれば、Newton-Raphson法(または割線法)を使用して実関数または複素関数のゼロを求めるとのこと。ではやってみます。

from scipy import optimize
import numpy as np
 
F = 1                 # バイオアベイラビリティ
D = 200               # 投与量(mg)
vd = 0.46 * 63.6      # 分布容積(L/kg) * 理想体重(kg)
ka = 0.36             # 吸収速度定数(hr^-1)
tau = 12              # 投与間隔(hr)
 

def func(ke, param, trough):
    _F = param[0]
    _D = param[1]
    _vd = param[2]
    _ka = param[3]
    _tau = param[4]
 
    _A = _F * _D * _ka / _vd / (_ka - ke)
 
    return _A * np.exp(-ke * _tau) / (1 - np.exp(-ke * _tau)) \
        - _A * np.exp(-_ka * _tau) / (1 - np.exp(-_ka * _tau)) - trough
 
param = [F, D, vd, ka, tau]
 
trough1 = 9               # 非喫煙時の実測トラフ値(μg/mL)
# optimize.newton(根を求めたい関数, 初期推定値, 追加の引数)
ke1 = optimize.newton(func, 0.1, args=(param, trough1))
print('ke1:', ke1)        # 0.05255675488882937
 
trough2 = 3               # 喫煙時の実測トラフ値(μg/mL)
ke2 = optimize.newton(func, 0.1, args=(param, trough2))
print('ke2:', ke2)        # 0.12123491149221788

コメントに書きましたがnewton()の引数は前から順に根を求めたい関数, 初期推定値, 追加の引数となっています（他にも導関数や許容誤差、反復の最大数などを指定できる）。
今度はうまく解が求まりました。上でみたグラフの値とも一致しているようにみえます。

求めたkeを使って繰り返し投与を描画

この推定したkeを使って繰り返し投与をプロットしてみました。描画部分のプログラムは前回記事で作成したものをほぼそのまま使っています。
f:id:enokisaute:20200315173951p:plain

color:  blue
Css_min:  8.999999999999266
------------------------------
color:  red
Css_min:  3.0000000000000018
------------------------------

定常状態でのCss_minも9と3μg/mLとなっています。
参考文献ではこの推定されたkeについて、個別化されたパラメータを使用していないため正確ではないとしつつも、喫煙によりクリアランスが２倍以上上昇したと考え、アミノフィリンの投与設計、維持投与量の考察へと話が続いていきます。

参考

  • 4ステップ 臨床力UPエクササイズ4TDM領域


最後に、今回のコードを載せておきます。

import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt
 
# 喫煙時と非喫煙時のトラフ値からkeを推定する
 
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             # バイオアベイラビリティ
    D = 200           # 投与量(mg)
    vd = 0.46 * 63.6  # 分布容積(L/kg) * 理想体重(kg)
    ka = 0.36         # 吸収速度定数(hr^-1)
    tau = 12          # 投与間隔(hr)
    days = 5          # 何日投与するか
 
    # optimize.newton()でf(ke)=0を解く
    def func(ke, param, trough):
        _F = param[0]
        _D = param[1]
        _vd = param[2]
        _ka = param[3]
        _tau = param[4]
 
        _A = _F * _D * _ka / _vd / (_ka - ke)
 
        return _A * np.exp(-ke * _tau) / (1 - np.exp(-ke * _tau)) \
            - _A * np.exp(-_ka * _tau) / (1 - np.exp(-_ka * _tau)) - trough
 
 
    param = [F, D, vd, ka, tau]
 
    trough1 = 9  # 非喫煙時の実測トラフ値(μg/mL)
    # newton(根を求めたい関数, 初期推定値, 追加の引数)
    ke1 = optimize.newton(func, 0.1, args=(param, trough1))
    print('ke1:', ke1)
 
    trough2 = 3  # 喫煙時の実測トラフ値(μg/mL)
    ke2 = optimize.newton(func, 0.1, args=(param, trough2))
    print('ke2:', ke2)
 
    # インスタンス生成
    po1 = PoOneComp2(F, D, vd, ka, ke1, tau)   # 非喫煙時(blue)
    po2 = PoOneComp2(F, D, vd, ka, ke2, tau)   # 喫煙時(red)
 
    # 描画区間の投与回数
    n = int(days * 24 / tau)
 
    conc1, conc2 = [], []
 
    # key: 色, value: (インスタンス, 投与回数, 血中濃度, 実測トラフ値, 推定したke)
    po = {'blue': (po1, n, conc1, trough1, ke1),
          'red': (po2, n, conc2, trough2, ke2)}
 
    period = 24 * days
    time = np.linspace(0, period, int(period / PoOneComp2.STEP), endpoint=False)
 
    for color, value in po.items():
        _po, n, conc, tgh, ke = 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='$\\rm C_{trough}$='+str(tgh) + 'μg/mLでフィッティング\nke='+str(round(ke,4)))
        print('color: ', color)
        print('Css_min: ', _po.get_cmin())
        print('-' * 30)
 
    plt.title('実測トラフ値からのke推定(イメージ図)')
    # 軸ラベルを表示
    plt.ylabel('血中濃度[μg/mL]')
    plt.xlabel('時間[hr]')
    plt.xlim(0, )
    # 目盛りを24時間毎にする
    plt.xticks(np.arange(0, time[-1] + 24, 24))
    plt.ylim(0, 16)
    plt.grid(linestyle='dashed')
    plt.legend()
    plt.show()