医薬品には併用してはならない組合せがあります。 このような併用してはならない組合せを併用禁忌・併用注意と呼び、 医療用医薬品添付文書への記載が義務付けられています[1]。 例えばインフルエンザの治療薬であるタミフルは、血栓症の治療薬であるワルファリンを併用注意として記載しています[2]。
処方されている医薬品が少なければ、併用禁忌・併用注意の回避は容易です。 しかし、日本では75歳以上の高齢者人口のうち四割は5種以上の医薬品を処方されており[3]、 併用禁忌・併用注意を避けるための医師・薬剤師の労力が大きくなっています。
このペーパーでは、併用禁忌・併用注意のリストから SBMにより併用禁忌・併用注意を避ける医薬品組合せの最大独立集合を求める一連の手順を説明します。 まず、仮想的に生成した併用禁忌・併用注意のリストを読み込みます。 次に、併用禁忌・併用注意を避けるよう、SBMへの入力となるQUBO行列を構築します[4]。 最後に、SBMを用いてQUBO行列を最小化する解ベクトルを求め、結果を評価します。
また、線形計画法ソルバーで解いた場合とSBMで解いた場合について大きい問題の処理時間を比較します。
import re
import numpy as np
import requests
import io
医薬品の種類数、併用不可リストの項目数を定義します。 併用不可リストの項目は医薬品のID二つからなっており、該当するIDの医薬品二つを同時に処方するのを禁止しています。 併用不可リストのIDを乱数で生成する際、併用不可リストのID二つが異なること、併用不可リストの中に同じ組合せのIDが登場しないこととします。
想定として、 医薬品の種類数にある医薬品はすべて、特定の患者にとって飲んだ方が良い医薬品だとします。 併用不可リストに記載されている医薬品の組合せには、 医療用医薬品添付文書に記載されている併用禁忌・併用注意の項目だけではなく、 同じ症状に対する医薬品二種などの好ましくない組合せも含むものとします。
np.random.seed(1234)
# 医薬品の種類数
N = 100
# 併用不可リストの項目数
M = 3000
Edge = np.zeros((N, N), dtype=int)
for m in range(M):
while True:
i = np.random.randint(N)
j = np.random.randint(N)
# 併用不可リストのID二つが異なること
# 併用不可リストの中に同じ組合せのIDが登場しないこと
if i != j and Edge[i, j] == 0:
break
Edge[i, j] += 1
print("医薬品の種類数", N)
print("併用禁忌・併用注意の項目数", M)
医薬品の種類数 100 併用禁忌・併用注意の項目数 3000
SBMへの入力となるQUBO行列を構築します。 併用禁忌・併用注意を避けたうえでなるべく多くの医薬品を有効にするという問題であり、 最大独立集合の一種として扱うことができます。 最大独立集合は多項式時間アルゴリズムが存在しておらず、 通常の逐次計算よりもSBMを使った方が高速に解ける可能性があります。
具体的には解ベクトルを以下のようになるようモデル化します。
QUBO行列の各要素を以下のように設定します。
非対角項の2という数字は、 すでに処方することを決めている医薬品$i$に対して 併用禁忌・併用注意となる医薬品$j$を追加するとき、 対角項の$-1$よりも非対角項の$+2$の方が大きくなるよう決めています。
# DIGITを大きくすると、ファイルサイズが大きくなる代わりに、丸め誤差が小さくなる
DIGIT = 8
def triangular(M):
if isinstance(M, np.ndarray):
return np.tril(M) + np.triu(M, k=+1).T
else:
return M
def writeMM(M):
MM = triangular(M)
N = MM.shape[0]
with io.StringIO() as f:
f.write('%%MatrixMarket matrix coordinate real symmetric\n')
f.write(f'{N} {N} {np.sum(MM!=0)}\n')
for i in range(N):
for j in range(N):
if MM[i, j] != 0:
f.write(f'{i+1} {j+1} {round(MM[i,j], DIGIT)}\n')
return f.getvalue()
解ベクトルを求め、評価します。 併用禁忌・併用注意を避けて医薬品を処方できたことが分かります。
%%time
Qubo = 2 * Edge.copy() - np.eye(N)
headers = {'content-type': 'application/octet-stream'}
params = {
"loops": 0,
"timeout": 0.1,
}
response = requests.post(f"http://aa.bb.cc.dd:8000/solver/autoising",
data=writeMM(Qubo),
headers=headers,
params=params)
X = np.array(response.json()["result"], dtype=bool)
print("選択できた薬品の数", sum(X))
print("併用禁忌・併用注意の該当数", Edge @ X @ X)
選択できた薬品の数 9 併用禁忌・併用注意の該当数 0 CPU times: total: 15.6 ms Wall time: 893 ms
同じ問題を線形計画法ソルバーでも解き、SBMで解いた場合と処理時間を比較します。 結果は以下のようになり、SBMの方が高速となりました。
| SBM v1.6 | Python MIP 1.13.0 / CBC | |
|---|---|---|
| 最適解までの処理時間 | 0.9 秒 | 26.8 秒 |
| SBM v1.6 (client) | SBM v1.6 (server) | Python MIP 1.13.0 / CBC | |
|---|---|---|---|
| CPU | Core-i7 8700 (3.20 GHz, 6 Core) | AWS EC2 p3.2xlarge (vCPU 8) | Core-i7 8700 (3.20 GHz, 6 Core) |
| GPU | on board | AWS EC2 p3.2xlarge (vCPU 8) | on board |
| RAM | 16 GB | 61 GB + 16 GB | 16 GB |
%%time
import mip
model = mip.Model()
x = model.add_var_tensor((N, ), 'x', var_type=mip.BINARY)
model.objective = mip.maximize(mip.xsum(x))
for (i, j) in zip(*np.where(Edge)):
model += (x[i] + x[j] <= 1)
if model.optimize() == mip.OptimizationStatus.OPTIMAL:
X = [x[n].x for n in range(N)]
print("選択できた薬品の数", sum(X))
print("併用禁忌・併用注意の該当数", Edge @ X @ X)
選択できた薬品の数 9.0 併用禁忌・併用注意の該当数 0.0 CPU times: total: 26.5 s Wall time: 26.8 s
