非線形薬物の血中濃度解析③ 可視化

非線形薬物の血中濃度解析を行うためのプログラム作成の続きです。作成の流れは、

1. データ入力

2. Lineweaver-Burkプロット（直線回帰）

3. 任意の投与量から血中濃度を計算

4. 可視化（グラフ化）

という順番で行っています。前回の投稿で1～2は作成済で、以下のコードが出来ました。

import numpy as np

num = int(input('採血ポイント数：'))
y = []
x = []
dose_list = []
conc_list = []

for i in range(num):
  dose = float(input('投与量：'))
  conc = float(input('血中濃度:'))
  dose_list.append(dose)
  conc_list.append(conc)
  y.append(1/dose)
  x.append(1/conc)

slope, intercept = np.polyfit(x, y, 1)
Vmax = 1/intercept
Km = slope * Vmax

今回はこの続きからです。

3. 任意の投与量から血中濃度を計算

KmとVmaxの定数を求めることができたので、投与量が分かれば血中濃度を予測することが出来ます。

適当な投与量を設定して1つ1つ計算していっても良いですが、せっかくなので全体像が分かるように網羅的に計算していきます。ループ処理を行えば、同じ計算の繰り返しも簡単にできそうですね。ループさせるためには、投与量のリストを作成する必要があるので、range()を使用します。range()を使用することで、指定した範囲内の等差数列を作成できます。

今回は、0～600の範囲で、5刻みの投与量リストを作成しdose_rangeという変数に格納します。

dose_range = range(0, 600, 5)

print(list(dose_range))

＞[0, 5, 10, 15, 20, 25, 30, 35, …, 585, 590, 595]
投与量のリストが完成しました。あとは、リストの要素1つ1つに対して繰り返し濃度計算を行い、計算結果をresultsというリストに格納します。

※注意※
計算にはKmとVmaxを使用するので、この記事の一番最初に示したコードを実行してから（サンプルデータは前回と同じ3点を使用します）下記コードを実行してください。

# 血中濃度を計算する関数を定義
def conc_calc(Km, Vmax, dose):
    return Km*dose/(Vmax-dose)

# 各投与量に対する血中濃度を計算し、リストに格納
results = []
for i in dose_range:
  conc = conc_calc(Km, Vmax, i)
  results.append(conc)

これで血中濃度の計算も終了しました。

4. 可視化（グラフ化）

データ作成と整形

投与量リスト（dose_range）とそれに対する血中濃度のリスト（results）が完成しているので、DataFrameを作成しグラフ化する準備をします。
pandasをimportし、{列名1：リスト1, 列名2：リスト2, …}というように辞書型でデータを渡してDataFrameを作成します。

import pandas as pd

df = pd.DataFrame({
    'dose(mg/d)': dose_range,
    'conc(mg/L)': results
    })

今回は非線形薬物を扱っているため、投与量が大きくなると血中濃度は際限なく大きくなってしまい治療域付近のデータが確認しにくくなってしまう可能性があります。例えば、ボリコナゾールは4 mg/L以上で肝障害が起こりやすくなるため、表示する血中濃度範囲は0～10 mg/Lぐらいで十分です。確認したい範囲に合わせてquery()を使用しデータを絞っておきます。

df = df.query("`conc(mg/L)` <= 10 and `conc(mg/L)` >= 0")

df.head()

df.tail()

血中濃度域を0～10 mg/Lに絞ることが出来ました。

上記のDataFrameとは別に、実測値のデータを集めたDataFrameも作成しておきます。曲線上に実測値をプロットするために使用します。作成方法は同じで、データ入力の時に作成したdose_listとconc_listを利用します。

actualValues = pd.DataFrame({
              'dose(mg/d)': dose_list,
              'conc(mg/L)': conc_list
              })

Altair

今回はAltairで可視化したいと思います。散布図（scatter plot）を使用します。

import altair as alt

# 曲線部分
scatter = alt.Chart(df).mark_line().encode(
            x='dose(mg/d)',
            y='conc(mg/L)',
            )

実測値のプロットには、alt.Chart( 1 ).mark_point( 2 ).encode( 3 )を使用します。

1. データソースとなるDataFrame（actualValues）を指定

2. プロットの設定（塗りつぶしやサイズ等）

3. 軸の設定

# 実測値プロット
points = alt.Chart(actualValues).mark_point(filled=False, size=50).encode(
        x=alt.X('dose(mg/d):Q', scale=alt.Scale(domain=[0,450])), 
        y=alt.Y('conc(mg/L):Q', scale=alt.Scale(domain=[0,10])),
        )

scatter + points

曲線と実測値のプロットが入ったグラフを描くことが出来ました。

おまけ

実務的には、「最大量としてどこまで増量できるか？」や「治療域をグラフ上に明示したい」などグラフに補助線を入れたくなります。また、投与量などのデータもテキストとしてまとめて出力されるとレポート作成時に便利です。
ボリコナゾールの解析を想定して、減量のボーダーラインとなる血中濃度4 mg/Lとそれに対する投与量を強調するようにコードを追記すると以下のようになります。

import numpy as np
import altair as alt
import pandas as pd

num = int(input('採血ポイント数：'))
y = []
x = []
dose_list = []
conc_list = []

for i in range(num):
  dose = float(input('投与量：'))
  conc = float(input('血中濃度:'))
  dose_list.append(dose)
  conc_list.append(conc)
  y.append(1/dose)
  x.append(1/conc)

slope, intercept = np.polyfit(x, y, 1)
Vmax = 1/intercept
Km = slope * Vmax

def conc_calc(Km, Vmax, dose):
    return Km*dose/(Vmax-dose)

#採血ポイント情報の表示
print('------------------------')
for i in range(num):
  print(f'データ{i+1}:投与量 {dose_list[i]}(mg/d)、血中濃度 {conc_list[i]}(mg/L)')

dose_range = range(0, 600, 5)
results = []
for i in dose_range:
  conc = conc_calc(Km, Vmax, i)
  results.append(conc)

#最大量を計算（conc=4）
limit_dose = Vmax*4/(Km+4)
print('------------------------')
print(f'Km = {Km:.2f}(mg/L)  \n'
      + f'Vmax = {Vmax:.2f}(mg/d)  \n'
      + f'最大投与量 = {limit_dose:.1f}(mg/d)'
      )

#プロットを作成
df = pd.DataFrame({
    'dose(mg/d)': dose_range,
    'conc(mg/L)': results
    })
#queryで 0<conc(mg/L)<10 の範囲のみへ
df = df.query("`conc(mg/L)` <= 10 and `conc(mg/L)` >= 0")

#実測ポイントプロット用のデータフレーム
actualValues = pd.DataFrame({
              'dose(mg/d)': dose_list,
              'conc(mg/L)': conc_list
              })
#曲線
scatter = alt.Chart(df).mark_line().encode(
            x='dose(mg/d)',
            y='conc(mg/L)',
            )
#補助線
rule = (alt.Chart().mark_rule(strokeDash=[5, 5], size=1, color='red').encode(y=alt.datum(4)))
#補助線2（最大投与量）
rule2 = (alt.Chart().mark_rule(strokeDash=[5, 5], size=1, color='red').encode(x=alt.datum(limit_dose)))
#実測ポイントをプロット
points = alt.Chart(actualValues).mark_point(filled=False, size=50).encode(
              x=alt.X('dose(mg/d):Q',
                      scale=alt.Scale(domain=[0,450])),
              y=alt.Y('conc(mg/L):Q',
                      scale=alt.Scale(domain=[0,10])),
               )

scatter + points + rule + rule2

---

計算するには少し面倒な非線形薬物の血中濃度解析について、Pythonを用いて計算ツールを作成しました。これがあれば投与量と血中濃度を入力するだけで自動的に計算と可視化が行われます。1つ1つのコードを理解すれば、自分用にアレンジすることも出来ます。是非ともオリジナルツール作りに挑戦してみてください！