JP2015181853A - 薬物動態パラメータの推定方法及び薬物動態パラメータの推定プログラム - Google Patents

Abstract

【課題】推定する薬物動態パラメータの数にかかわらず、当該薬物動態パラメータを発散の可能性を低く抑えつつ推定できる方法を提供する。【解決手段】患者の血液中の薬物の濃度から薬物動態パラメータxiを推定する方法であり、電子計算機のCPUにより実行される。CPUは、データ入力画面３１における入力データなどに基づいて、薬物動態パラメータxiを推定する方法として最尤法またはベイズ推定を選択する。CPUは、最尤法またはベイズ推定の基本式である薬物動態パラメータxiのベクトルx=(x1,x2,・・・,xn)Tの関数s(x)に対して２階偏微分を実行し、当該２階偏微分によって得られた式から同一の薬物動態パラメータxiによる２階偏導関数の項を削除した近似式に修正マルカート法を適用して、関数s(x)が最小となるベクトルxを求めることで、薬物動態パラメータxiを推定する。【選択図】図３

Description

本発明は、患者から採取された血液における薬物の血中濃度から薬物動態パラメータを推定する方法、及び当該方法を実行するための薬物動態パラメータの推定プラグラムに関する。

個々の患者に適した薬物の投与設計を行う治療薬物モニタリング(Therapeutic Drug Monitoring,TDM)が知られている。TDMでは、患者の血液における薬物の血中濃度を測定して解析した結果に基づいて、最適な薬用量、投与法が設定される。また、TDMでは、測定された血中濃度から薬物動態パラメータが推定される。

上述したようなTDMを実施するためのシステムの一例として、ベイズ推定を用いて予測血中濃度を計算する薬物投与解析システムが、特許文献１に開示されている。また、薬物動態パラメータの推定には、数値計算法の一つであるニュートン法やガウス・ニュートン法が用いられることが多い。

特開平７−１２４１２５号公報

ベイズ推定は、一人の患者に対する血中濃度の測定回数が少なく、測定回数が推定するべき薬物動態パラメータの数以下である場合に、過去のデータから得られた母集団パラメータの値を使用して薬物動態パラメータを推定する数値計算法である。しかし、一人の患者に対する血中濃度の測定回数が多いときには、母集団パラメータを用いずとも薬物動態パラメータを推定することができるため、ベイズ推定以外の数値計算法を用いる方がよい場合がある。

また、ニュートン法は、方程式f(x)=0におけるxの値の近似値を反復することによって求める手法である。しかし、ニュートン法では、方程式に最初に代入する初期値x⁽⁰⁾によっては、xの値が収束せずに発散してしまう。そして、ニュートン法では、xの値を収束させるための初期値x⁽⁰⁾の収束域が狭いため、xの値を発散させない適切な初期値x⁽⁰⁾を選択することが困難であるという問題がある。

また、ニュートン法よりも収束性を改善させるための手法として、ガウス・ニュートン法がある。しかし、ガウス・ニュートン法は、２乗和の形の関数に適した手法であるため、２乗和以外の項がある関数に適用した場合、収束性の改善は期待できない。

本発明は、上記事情に鑑みてなされたものであり、その目的は、推定する薬物動態パラメータの数にかかわらず、当該薬物動態パラメータを、発散の可能性を低く抑えつつ推定することができる方法、及び当該方法を実行するための薬物動態パラメータの推定プラグラムを提供することにある。

(1) 本発明に係る薬物動態パラメータの推定方法は、患者から採取された血液における薬物の血中濃度から薬物動態パラメータx_i(i=1,2,・・・,n)を推定する方法であって、薬物動態パラメータx_iを推定する方法として最尤法またはベイズ推定の何れかが選択される方法選択ステップと、上記方法選択ステップにおいて選択された方法における基本式である上記薬物動態パラメータx_iのベクトルx=(x₁,x₂,・・・,x_n)^Tの関数s(x)に対して２階偏微分を実行し、当該２階偏微分によって得られた式から同一の薬物動態パラメータx_iによる２階偏導関数の項を削除した近似式に修正マルカート法を適用して上記関数s(x)を最小化し、最小化した際の上記ベクトルxを求めることによって、上記薬物動態パラメータx_iを推定する推定ステップと、を含む。

本方法は、方法選択ステップを含む。これにより、血中濃度の測定回数が推定すべき薬物動態パラメータx_iの個数よりも多い場合、過去のデータから得られた母集団パラメータを用いることなく、最尤法によって薬物動態パラメータx_iを推定することができる。一方、血中濃度の測定回数が推定すべき薬物動態パラメータx_iの個数以下の場合、母集団パラメータを用いて、ベイズ推定によって薬物動態パラメータx_iを推定することができる。

また、本方法では、修正マルカート法が用いられている。このため、ニュートン法よりも収束性を改善することができる。また、修正マルカート法は、ガウス・ニュートン法とは異なり、２乗和以外の項がある関数に適用した場合でも、収束性を改善することができる。

また、本願の出願人は、最尤法及びベイズ推定の基本式を２階偏微分して得られた式から２階偏導関数の項を削除した近似式を修正マルカート法を適用することによって、関数s(x)が発散する可能性を低くすることができることを発見した。

(2) 本発明に係る薬物動態パラメータの推定方法において、最尤法における上記関数s(x)は、時刻t_jにおける血中濃度測定値c_jが正規分布N(c(t_j),σ_c ²)に従う場合、以下の（式１）であり、上記血中濃度測定値c_jが正規分布N(c(t_j),{ω_c×c(t_j)}²)に従う場合、以下の（式２）であり、上記血中濃度測定値c_jの対数値log(c_j)が正規分布N(log{c(t_j)},(σ-_c ²))(但し、"-"は文字σのアッパーラインの意味である。)に従う場合、以下の（式３）である。

以上において、
m：血中濃度測定回数
t_j：j回目の測定時刻(j=1,2,・・・,m)
c_j：時刻t_jにおける血中濃度測定値
c(t_j)=c(t_j,x)：時刻t_jにおける血中濃度理論値
σ-_c(但し、"-"は文字σのアッパーラインの意味である。)：jによらず同じ値をとる場合のlog(c_j)の標準偏差
ω_c：σ_j／c(t_j,x)で与えられる変動係数
である。上記方法選択ステップにおいて最尤法が選択された場合に、当該３つの式の何れかが選択される式選択ステップを更に含む。

従来、最尤法による薬物動態パラメータの推定において、関数s(x)として（式１）が用いられていた。しかしながら、（式１）を用いた場合、推定した薬物動態パラメータの精度の点に問題があり、精度の向上が求められていた。この問題の解決のために、関数s(x)として（式２）を用いることが考えられる。しかしながら、（式２）を用いた場合、xの値の収束性が良くないという問題がある。つまり、推定過程においてxの値の収束が止まってしまうことがある。xの値の収束が止まってしまった場合、収束したxの値と実測値とが乖離してしまう。

そこで、本願発明者は上記（１）のような薬物パラメータの推定方法を考え出した。当該推定方法によると、（式２）を用いた場合でも、xの値の収束性を良好にすることができる。その結果、本方法のように（式１）及び（式２）を選択的に使用することができ、しかも何れの式を使用した場合であっても、薬物動態パラメータを正確に推定することができる。

また、本方法によれば、血中濃度測定値が正規分布に従う場合（（式１）及び（式２））と対数正規分布に従う場合（式３）とに応じて式を使い分けることができる。また、血中濃度測定値が正規分布に従う場合には更に、血中濃度測定値の誤差の分散が一定の場合（式１）と、血中濃度測定値の誤差の分散が血中濃度に比例する場合（式２）とに応じて式を使い分けることができる。

(3) 本発明に係る薬物動態パラメータの推定方法において、ベイズ推定における上記関数s(x)は、時刻t_jにおける血中濃度測定値c_jが正規分布N(c(t_j),σ_c ²)に従い且つ薬物動態パラメータx_iが正規分布N(μ_j,σ_i ²)に従う場合、以下の（式４）であり、上記血中濃度測定値c_jが正規分布N(c(t_j),{ω_c×c(t_j)}²)に従い且つ薬物動態パラメータx_iが正規分布N(μ_j,(ω_i×μ_i)²)に従う場合、以下の（式５）であり、上記血中濃度測定値c_jの対数値log(c_j)が正規分布N(log{c(t_j)},(σ-_c ²))(但し、"-"は文字σのアッパーラインの意味である。)に従い且つ薬物動態パラメータx_iの対数値log(x_i)が正規分布N(μ-_i,σ-_i ²)(但し、"-"は文字μ,σそれぞれのアッパーラインの意味である。)に従う場合、以下の（式６）である。

以上において、
m：血中濃度測定回数
n：薬物動態パラメータ数
t_j：j回目の測定時刻(j=1,2,・・・,m)
c_j：時刻t_jにおける血中濃度測定値
c(t_j)=c(t_j,x)：時刻t_jにおける血中濃度理論値
σ_c：jによらず同じ値をとる場合のc_jの標準偏差
σ-_c(但し、"-"は文字σのアッパーラインの意味である。)：jによらず同じ値をとる場合のlog(c_j)の標準偏差
μ_i：x_i(i=1,2,・・・,n)の平均
μ-_i(但し、"-"は文字μのアッパーラインの意味である。)：log(x_i)の平均
σ_i：x_iの標準偏差
σ-_i(但し、"-"は文字σのアッパーラインの意味である。)：log(x_i)の標準偏差
ω_c：σ_c／c(t_j,x)で与えられる変動係数
ω_i：σ_i／μ_iで与えられる変動係数
ｘ-_i(但し、"-"は文字xのアッパーラインの意味である。)：log(x_i)
である。上記方法選択ステップにおいてベイズ推定が選択された場合に、当該３つの式の何れかが選択される式選択ステップを更に含む。

従来、ベイズ推定による薬物動態パラメータの推定において、関数s(x)として（式４）が用いられていた。しかしながら、（式４）を用いた場合、推定した薬物動態パラメータの精度の点に問題があり、精度の向上が求められていた。この問題の解決のために、関数s(x)として（式５）を用いることが考えられる。しかしながら、（式５）を用いた場合、xの値の収束性が良くないという問題がある。つまり、推定過程においてxの値の収束が止まってしまうことがある。xの値の収束が止まってしまった場合、収束したxの値と実測値とが乖離してしまう。

そこで、本願発明者は上記（１）のような薬物パラメータの推定方法を考え出した。当該推定方法によると、（式４）を用いた場合でも、xの値の収束性を良好にすることができる。その結果、本方法のように（式４）及び（式５）を選択的に使用することができ、しかも何れの式を使用した場合であっても、薬物動態パラメータを正確に推定することができる。

また、本方法によれば、血中濃度の測定において混入した誤差及び薬物動態パラメータが正規分布に従う場合（（式４）及び（式５））と対数正規分布に従う場合（式６）とに応じて式を使い分けることができる。また、上記誤差が正規分布に従う場合には更に、血中濃度測定値の誤差の分散が一定の場合（式４）と、血中濃度測定値の誤差の分散が血中濃度に比例する場合（式５）とに応じて式を使い分けることができる。

(4) 本発明に係る薬物動態パラメータの推定方法において、上記推定ステップは、上記ベクトルxの初期値x⁽⁰⁾=(x₁ ⁽⁰⁾,x₂ ⁽⁰⁾,・・・,x_n ⁽⁰⁾)^Tを準備し、最小値が更新されたか否かを表す指標であるrの初期値をr=0と準備する初期値準備ステップと、上記式選択ステップにおいて選択された式から、上記ベクトルxのi（iは０以上の整数）番目の近似値であるx⁽ⁱ⁾=(x₁ ⁽ⁱ⁾,x₂ ⁽ⁱ⁾,・・・,x_n ⁽ⁱ⁾)^T、及びΔx=(Δx₁,Δx₂,・・・,Δx_n）^Tを含む行列による式を導出し、当該行列による式において上記Δxに対するn元連立一次方程式を解き、x'=x⁽ⁱ⁾+Δxとおく計算ステップと、上記計算ステップにおいて上記n元連立一次方程式を解いた結果、s(x')<s(x⁽ⁱ⁾)であることを条件として、rをmin(0,r-1)とし且つiに１を加えてx⁽ⁱ⁾=x'とし、s(x')>=s(x⁽ⁱ⁾)であることを条件として、rに１を加える変数変更ステップと、上記変数変更ステップの実行後に、i>=Mまたはr>=N（但し、M,Nは予め設定された正の整数定数）である第１条件を具備することを条件として、s(x)が収束しないと判定し、|s(x⁽ⁱ⁺¹⁾)-s(x⁽ⁱ⁾)|<ε且つ|Δx|<δ且つr=0（但し、ε,δは予め設定された実数正数）である第２条件を具備することを条件として、s(x)が収束すると判定する判定ステップと、上記判定ステップにおいて、s(x)が収束すると判定されたことを条件として、そのときのx⁽ⁱ⁾の各成分を薬物動態パラメータx_iとする決定ステップと、を含む。上記計算ステップは、上記初期値準備ステップの実行が完了したこと、または、上記判定ステップにおいて上記第１条件及び上記第２条件を具備しないことを条件として実行される。

本方法のように、行列による式を用いることによって、計算が容易となる。

(5) 本発明に係る薬物動態パラメータの推定方法において、上記推定ステップは、上記初期値準備ステップにおいて、上記ベクトルxの初期値x⁽⁰⁾を複数個準備する。上記決定ステップは、複数個の初期値x⁽⁰⁾の各々について上記計算ステップ、上記変数変更ステップ、及び上記判定ステップが実行された結果、s(x)が最も小さく収束したときのx⁽ⁱ⁾の各成分を薬物動態パラメータx_iと決定する。

薬物動態パラメータx_iの推定過程において、xが発散してしまうと、薬物動態パラメータx_iを推定することができない。また、xが収束した場合でも、当該収束がs(x)を最小にする解以外の局所解に収束した場合には、適切な薬物動態パラメータx_iを推定することができていない。そこで、本方法のように、初期値x⁽⁰⁾を複数個準備して、複数の関数s(x)の解を求め、当該解が最も小さく収束したときのx⁽ⁱ⁾の各成分を薬物動態パラメータx_iと決定することで、局所解ではない適切な薬物動態パラメータx_iを推定することができる。

(6) 本発明に係る薬物動態パラメータの推定方法では、上記初期値準備ステップにおいて、上記ベクトルxの初期値x⁽⁰⁾は、ボックス・ミューラー法によって作成される。

ボックス・ミューラー法を用いることにより、作成された初期値x⁽⁰⁾の各成分を正規分布に従う独立な正規乱数とすることができる。特に、初期値x⁽⁰⁾を複数個準備した場合、複数の解xを求めた際に、当該解の少なくとも一部が適切に収束する可能性を高くすることができる。

(7) 本発明は、患者から採取された血液における薬物の血中濃度から薬物動態パラメータx_i(i=1,2,・・・,n)を推定するためのプログラムであって、薬物動態パラメータx_iを推定する方法として最尤法またはベイズ推定の何れかが選択される方法選択ステップと、上記方法選択ステップにおいて選択された方法における基本式である上記薬物動態パラメータx_iのベクトルx=(x₁,x₂,・・・,x_n)^Tの関数s(x)に対して２階偏微分を実行し、当該２階偏微分によって得られた式から同一の薬物動態パラメータx_iによる２階偏導関数の項を削除した近似式に修正マルカート法を適用して上記関数s(x)を最小化し、最小化した際の上記ベクトルxを求めることによって、上記薬物動態パラメータx_iを推定する推定ステップと、を実行させる薬物動態パラメータの推定プログラムとして捉えることもできる。

(8) 本発明に係る薬物動態パラメータの推定プログラムにおいて、最尤法における上記関数s(x)は、時刻t_jにおける血中濃度測定値c_jが正規分布N(c(t_j),σ_c ²)に従う場合、以下の（式１）であり、上記血中濃度測定値c_jが正規分布N(c(t_j),{ω_c×c(t_j)}²)に従う場合、以下の（式２）であり、上記血中濃度測定値c_jの対数値log(c_j)が正規分布N(log{c(t_j)},(σ-_c ²))(但し、"-"は文字σのアッパーラインの意味である。)に従う場合、以下の（式３）である。

以上において、
m：血中濃度測定回数
t_j：j回目の測定時刻(j=1,2,・・・,m)
c_j：時刻t_jにおける血中濃度測定値
c(t_j)=c(t_j,x)：時刻t_jにおける血中濃度理論値
σ-_c(但し、"-"は文字σのアッパーラインの意味である。)：jによらず同じ値をとる場合のlog(c_j)の標準偏差
ω_c：σ_j／c(t_j,x)で与えられる変動係数
である。上記方法選択ステップにおいて最尤法が選択された場合に、当該３つの式の何れかが選択される式選択ステップを更に実行させる。

(9) 本発明に係る薬物動態パラメータの推定プログラムにおいて、ベイズ推定における上記関数s(x)は、時刻t_jにおける血中濃度測定値c_jが正規分布N(c(t_j),σ_c ²)に従い且つ薬物動態パラメータx_iが正規分布N(μ_j,σ_i ²)に従う場合、以下の（式４）であり、上記血中濃度測定値c_jが正規分布N(c(t_j),{ω_c×c(t_j)}²)に従い且つ薬物動態パラメータx_iが正規分布N(μ_j,(ω_i×μ_i)²)に従う場合、以下の（式５）であり、上記血中濃度測定値c_jの対数値log(c_j)が正規分布N(log{c(t_j)},(σ-_c ²))(但し、"-"は文字σのアッパーラインの意味である。)に従い且つ薬物動態パラメータx_iの対数値log(x_i)が正規分布N(μ-_i,σ-_i ²)(但し、"-"は文字μ,σそれぞれのアッパーラインの意味である。)に従う場合、以下の（式６）である。

以上において、
m：血中濃度測定回数
n：薬物動態パラメータ数
t_j：j回目の測定時刻(j=1,2,・・・,m)
c_j：時刻t_jにおける血中濃度測定値
c(t_j)=c(t_j,x)：時刻t_jにおける血中濃度理論値
σ_c：jによらず同じ値をとる場合のc_jの標準偏差
σ-_c(但し、"-"は文字σのアッパーラインの意味である。)：jによらず同じ値をとる場合のlog(c_j)の標準偏差
μ_i：x_i(i=1,2,・・・,n)の平均
μ-_i(但し、"-"は文字μのアッパーラインの意味である。)：log(x_i)の平均
σ_i：x_iの標準偏差
σ-_i(但し、"-"は文字σのアッパーラインの意味である。)：log(x_i)の標準偏差
ω_c：σ_c／c(t_j,x)で与えられる変動係数
ω_i：σ_i／μ_iで与えられる変動係数
ｘ-_i(但し、"-"は文字xのアッパーラインの意味である。)：log(x_i)
である。上記方法選択ステップにおいてベイズ推定が選択された場合に、当該３つの式の何れかが選択される式選択ステップを更に実行させる。

(10) 本発明に係る薬物動態パラメータの推定プログラムにおいて、上記推定ステップは、上記ベクトルxの初期値x⁽⁰⁾=(x₁ ⁽⁰⁾,x₂ ⁽⁰⁾,・・・,x_n ⁽⁰⁾)^Tを準備し、最小値が更新されたか否かを表す指標であるrの初期値をr=0と準備する初期値準備ステップと、上記式選択ステップにおいて選択された式から、上記ベクトルxのi（iは０以上の整数）番目の近似値であるx⁽ⁱ⁾=(x₁ ⁽ⁱ⁾,x₂ ⁽ⁱ⁾,・・・,x_n ⁽ⁱ⁾)^T、及びΔx=(Δx₁,Δx₂,・・・,Δx_n）^Tを含む行列による式を導出し、当該行列による式において上記Δxに対するn元連立一次方程式を解き、x'=x⁽ⁱ⁾+Δxとおく計算ステップと、上記計算ステップにおいて上記n元連立一次方程式を解いた結果、s(x')<s(x⁽ⁱ⁾)であることを条件として、rをmin(0,r-1)とし且つiに１を加えてx⁽ⁱ⁾=x'とし、s(x')>=s(x⁽ⁱ⁾)であることを条件として、rに１を加える変数変更ステップと、上記変数変更ステップの実行後に、i>=Mまたはr>=N（但し、M,Nは予め設定された正の整数定数）である第１条件を具備することを条件として、s(x)が収束しないと判定し、|s(x⁽ⁱ⁺¹⁾)-s(x⁽ⁱ⁾)|<ε且つ|Δx|<δ且つr=0（但し、ε,δは予め設定された実数正数）である第２条件を具備することを条件として、s(x)が収束すると判定する判定ステップと、上記判定ステップにおいて、s(x)が収束すると判定されたことを条件として、そのときのx⁽ⁱ⁾の各成分を薬物動態パラメータx_iとする決定ステップと、を実行させる。上記計算ステップは、上記初期値準備ステップの実行が完了したこと、または、上記判定ステップにおいて上記第１条件及び上記第２条件を具備しないことを条件として実行される。

(11) 本発明に係る薬物動態パラメータの推定プログラムにおいて、上記推定ステップは、上記初期値準備ステップにおいて、上記ベクトルxの初期値x⁽⁰⁾を複数個準備する。上記決定ステップは、複数個の初期値x⁽⁰⁾の各々について上記計算ステップ、上記変数変更ステップ、及び上記判定ステップが実行された結果、s(x)が最も小さく収束したときのx⁽ⁱ⁾の各成分を薬物動態パラメータx_iと決定する。

(12) 本発明に係る薬物動態パラメータの推定プログラムでは、上記初期値準備ステップにおいて、上記ベクトルxの初期値x⁽⁰⁾は、ボックス・ミューラー法によって作成される。

本発明に係る薬物動態パラメータの推定方法及び薬物動態パラメータの推定プログラムによると、推定する薬物動態パラメータx_iの数にかかわらず、当該薬物動態パラメータx_iを、発散の可能性を低く抑えつつ推定することができる。

図１は、電子計算機１０のブロック図である。 図２は、薬物動態パラメータを推定するためのシステムの処理について説明するためのフローチャートである。 図３は、表示部１６に表示されたデータ入力画面３１を示す図である。 図４は、表示部１６に表示されたモデル選択画面３５を示す図である。 図５は、表示部１６に表示されたモデル選択画面３６を示す図である。 図６は、表示部１６に表示された血中濃度推移画面４１を示す図である。 図７は、推定ステップの処理について説明するためのフローチャートである。 図８は、変形例１における表示部１６に表示されたデータ入力画面３１を示す図である。 図９は、変形例３における推定ステップの処理について説明するためのフローチャートである。

以下、適宜図面を参照して本発明の実施形態について説明する。なお、以下に説明される実施形態は本発明の一例にすぎず、本発明の要旨を変更しない範囲で、本発明の実施形態を適宜変更できることは言うまでもない。

［薬物動態パラメータを推定するためのシステムの構成］
本発明に係る薬物動態パラメータの推定方法によって推定される薬物動態パラメータは、コンパートメントモデルに追随して決まる。ここで、コンパートメントモデルは、薬物動態に関する体内のモデルである。コンパートメントモデルは、対象となる薬剤によって決まる。

コンパートメントモデルとしては、例えば点滴静脈注射のコンパートメントモデル、経口・筋肉注射のコンパートメントモデル、及び点滴のコンパートメントモデルなどがある。各コンパートメントモデルは、定義された複数の薬物動態パラメータと、当該複数の薬物動態パラメータ間で成立する少なくとも１つの関係式よりなる。つまり、コンパートメントモデルが決定されると、使用する未知の薬物動態パラメータのパラメータ数が決まる。なお、コンパートメントモデルは、上述した３種類に限らず、他にも存在する。また、同じ種類のコンパートメントモデルであっても、式の変換により別の薬物動態パラメータが定義されることもある。

以下、点滴のコンパートメントモデルの１つについて、具体的に説明する。この点滴のコンパートメントモデルの１つ（例示コンパートメントモデルと称する。）では、以下の７つの薬物動態パラメータが定義される。つまり、CL:全身クリアランス(単位はリットル／時間(L/hr))、V₁:体循環コンパートメントの分布容積(L)、V₂:末梢コンパートメントの分布容積(L)、Q:体循環コンパートメントと末梢コンパートメントの間のクリアランス(L/hr)、k₁₀:消失速度定数(1/hr)、k₁₂:体循環から末梢コンパートメントへの薬物移行を表す消失速度定数(1/hr)、k₂₁:末梢から体循環コンパートメントへの薬物移行を表す消失速度定数(1/hr)である。また、例示コンパートメントモデルの薬物動態パラメータでは、以下の関係式が成り立つ。つまり、CL=k₁₀V₁,Q=k₁₂V₁,Q=k₂₁V₂である。

この場合、７つのパラメータに対して関係式が３つであるため、未知のパラメータ数は４つとなる。本発明に係る薬物動態パラメータの推定方法によって推定する薬物動態パラメータの選択については、原則として任意であるが、コンパートメントモデルが有する関係式によって選択可能な薬物動態パラメータが制限される場合がある。例えば、例示コンパートメントモデルの場合、CL=k₁₀V₁,関係式が存在するために、推定する薬物動態パラメータとして、(V₁,k₁₀,k₁₂,k₂₁)の４つを選択することはできるが、(V₁,k₁₀,CL,k₂₁)の４つを選択することはできない。

本発明に係る方法による薬物動態パラメータの推定は、図１に示されるパーソナルコンピュータなどの電子計算機１０で構成されたシステムが当該方法のプログラムを実行することによって実現される。つまり、上記プログラムは、本発明に係る薬物動態パラメータの推定プログラムに相当する。また、本実施形態において、上記プログラムは、後述するようにRAM１２に格納されているが、上記プログラムは、ＣＤ−ＲＯＭなどの記録媒体に格納されていてもよい。つまり、本発明に係る薬物動態パラメータの推定プログラムは、上記プログラムが記録されたＣＤ−ＲＯＭなどの記録媒体として捉えられてもよい。

電子計算機１０は、中央演算装置(CPU)１１とRAM１２と大容量記憶部１３とASIC１４とを有する計算機本体２０、入力部１５、及び表示部１６を備えている。また、計算機本体２０は、CPU１１、RAM１２、大容量記憶部１３、及びASIC１４を相互に接続する内部バス２１を備えている。

CPU１１は、RAM１２に格納されたプログラムに従って演算を実行する。

大容量記憶部１３は、大容量の情報を格納可能な記憶媒体であり、ハードディスク、光磁気ディスク、フラッシュメモリなどで構成されている。大容量記憶部１３には、本発明に係る方法によって薬物動態パラメータの推定を実行可能なプログラムや、血中濃度に関する測定データなどが格納されている。大容量記憶部１３に格納された情報は、電子計算機１０の電源がオフとなっても消去されない。

RAM１２は、大容量記憶部１３よりも小容量の情報を格納可能な記憶媒体である。RAM１２には、上記プログラム、及びCPU１１が上記プログラムを実行する際に用いるデータや信号などが、一時的に格納される。

ASIC１４は、入力部１５及び表示部１６と接続されている。ASIC１４には、入力部１５及び表示部１６とCPU１１との信号のやり取りを仲介する回路が組み込まれている。

入力部１５は、後述する患者に関する個人情報、患者に投与する医薬品情報、及び患者の血中濃度に関する測定データを入力する際に使用されるキーボード１７やマウス１８などで構成されている。

表示部１６は、液晶ディスプレイなどで構成されており、入力部１５によって入力された個人情報、医薬品情報、及び測定データや、後述する推定ステップにおいてCPU１１による演算によって推定された薬物動態パラメータに基づく血中濃度推移線などが表示される（図３〜図６参照）。

［薬物動態パラメータを推定するためのシステムの処理の概要］
CPU１１がRAM１２に格納された上記プログラムを実行することによって、薬物動態パラメータを推定する処理が図２のフローチャートに従って実行される。以下に、当該処理の概要が説明される。

電子計算機１０のユーザがマウス１８を操作することによって、表示部１６に表示された上記プログラムを実行するためのアイコン（不図示）がダブルクリックされると、上記プログラムが大容量記憶部１３からRAM１２にコピーされて、CPU１１がRAM１２にコピーされた上記プログラムを実行する。これにより、表示部１６には、図３に示されるようなデータ入力画面３１が表示される。データ入力画面３１には、患者に関する個人情報が入力される個人情報入力領域８１、患者に投与する医薬品情報が入力される医薬品情報入力領域８２、患者の血中濃度に関する測定データが入力される測定データ入力領域８３、及び「データ入力完了」ボタン３４が設けられている。

ユーザは、データ入力画面３１を参照しつつ入力部１５（キーボード１７及びマウス１８）を操作することによって、上述した各領域８１、８２、８３に、個人情報、医薬品情報、及び測定データを入力する（図２のステップＳ１０）。

図３に示されるように、個人情報入力領域８１に入力される個人情報は、患者の性別、年齢、身長、体重である。また、医薬品情報入力領域８２に入力される医薬品情報は、患者に投与される医薬品の品名、剤形などである。なお、個人情報入力領域８１及び医薬品情報入力領域８２は、上述した情報以外の情報が入力されるように構成されていてもよい。

また、測定データ入力領域８３に入力される測定データは、患者の血中濃度を測定した時刻と、当該時刻における血中濃度とで構成される。例えば、血中濃度が１５分間隔且つ２４時間に亘って測定された場合、測定時刻と当該測定時刻に対応する血中濃度とのデータ数は、それぞれ９６個である。図３では、測定データは、キーボード１７の操作による測定時刻と血中濃度との入力と、マウス１８の操作による「次データの入力」ボタン３２のクリックとが交互に実行されることによって、測定時刻と血中濃度とが１組ずつ入力される。そして、「測定データ入力完了」ボタン３３がクリックされることによって、今まで入力された複数組の測定時刻と血中濃度とが測定データとして確定される。なお、測定データは、上述のように１組ずつ入力されるのではなく、大容量記憶部１３に予めデータベースとして格納されていてもよい。

ユーザは、個人情報、医薬品情報、及び測定データを入力を完了すると、マウス１８を操作することによって、「データ入力完了」ボタン３４をクリックする。これにより、後述する方法選択ステップが実行され（図２のステップＳ２０）、薬物動態パラメータを推定する方法として、最尤法またはベイズ推定が選択される。最尤法が選択された場合、表示部１６に、図４に示されるようなモデル選択画面３５が表示される。一方、ベイズ推定が選択された場合、表示部１６に、図５に示されるようなモデル選択画面３６が表示される。

モデル選択画面３５には、最尤法において選択可能な３つのモデル（正規分布モデル１、正規分布モデル２、及び対数正規分布モデル）が表示された最尤法モデル選択領域８４、「推定実行」ボタン３７、及び「戻る」ボタン３９が設けられている。また、モデル選択画面３６には、ベイズ推定において選択可能な３つのモデル（正規分布モデル１、正規分布モデル２、及び対数正規分布モデル）が表示されたベイズ推定モデル選択領域８５、「推定実行」ボタン３８、及び「戻る」ボタン４０が設けられている。

ユーザは、マウス１８の操作によって、各モデル選択領域８４、８５に表示された３つのモデルの何れかを選択し、「推定実行」ボタン３７、３８をクリックする。これにより、後述する式選択ステップが実行される（図２のステップＳ３０）。なお、モデル選択画面３５、３６の「戻る」ボタン３９、４０がクリックされると、表示部１６に表示される画面がモデル選択画面３５、３６からデータ入力画面３１へと遷移する。

ステップＳ３０に続けて、ステップＳ３０において選択された式に基づいて、後述する推定ステップが実行される（図２のステップＳ４０）。推定ステップが実行されることにより、CPU１１が後述する演算を実行する。これにより、薬物動態パラメータが推定される。

ステップＳ４０が実行されて薬物動態パラメータが推定されると、図６に示される血中濃度推移画面４１が表示される。血中濃度推移画面４１には、描画領域８６及び「戻る」ボタン４４が設けられている。

描画領域８６は、横軸が時刻（日時）であり縦軸が血中濃度(mEq/L)であるグラフを描画するための領域である。描画領域８６には、２本の血中濃度推移線４２、４３が描画される（図２のステップＳ５０）。

図６中に細線で示された血中濃度推移線４２は、現行の薬物の投与スケジュールに基づく血中濃度推移線である。現行の薬物の投与スケジュールに基づく血中濃度推移線は、ステップＳ１０において入力された測定データや、過去の多数の患者に対する測定データなどに基づいて描画されるグラフである。

図６中に太線で示された血中濃度推移線４３は、ステップＳ１０において入力された個人情報、医薬品情報、及び測定データや、ステップＳ４０において推定された薬物動態パラメータに基づいて描画されるグラフであり、現行の薬物の投与スケジュールに基づく血中濃度推移線よりも、医薬品を投与される患者に適した投与スケジュールに基づく血中濃度推移線である。

なお、各種情報、各種データ、及びステップＳ４０において推定された薬物動態パラメータなどから血中濃度推移線４２、４３を導き出す過程は、ここでは、その詳細な説明を省略する。

血中濃度推移画面４１の「戻る」ボタン４４がクリックされると、表示部１６に表示される画面が血中濃度推移画面４１から直近に表示されていたモデル選択画面３５、３６へと遷移する。

［方法選択ステップ］
以下、図２のステップＳ２０において実行される方法選択ステップが説明される。方法選択ステップは、CPU１１によって実行される。CPU１１は、未知の薬物動態パラメータの個数と、測定回数とを比較する。

未知の薬物動態パラメータの個数は、対象となった薬剤によって決まったコンパートメントモデルが有する複数の薬物動態パラメータのうち、当該コンパートメントモデルが有する関係式によって制限されていない範囲で、推定すべき薬物動態パラメータとして選択された薬物動態パラメータの個数である。また、測定回数は、図２のステップＳ１０において入力された測定データの組数である。

CPU１１は、測定回数が未知の薬物動態パラメータの個数よりも多い場合、最尤法を選択し、表示部１６にモデル選択画面３５（図４参照）を表示させる。このような場合に最尤法が選択される理由は、このような場合においては方程式の数を変数の個数以上とすることができるために、連立方程式によって未知の薬物動態パラメータを求めることができるからである。よって、測定回数が未知の薬物動態パラメータの個数よりも多い場合、連立方程式を解くことによって変数の値、つまり薬物動態パラメータを求めることができる。

一方、CPU１１は、測定回数が未知の薬物動態パラメータの個数以下であり、薬物動態パラメータの母集団パラメータがあり、且つ当該母集団パラメータが互いに独立である場合、ベイズ推定を選択し、表示部１６にモデル選択画面３６（図５参照）を表示させる。ここで、薬物動態パラメータの母集団パラメータは、過去における多数の患者の血中濃度の測定から求められた薬物動態パラメータのデータ群であり、大容量記憶部１３に格納されている。

なお、測定回数が未知の薬物動態パラメータの個数以下であるが、薬物動態パラメータの母集団パラメータがない場合、或いは当該母集団パラメータがあるが当該母集団パラメータに含まれる薬物動態パラメータが互いに独立でない場合、CPU１１は、最尤法及びベイズ推定の何れも選択せずにその旨のメッセージを表示部１６に表示させてもよいし、ベイズ推定ではなく最尤法を選択するために更に測定データを入力して測定回数を増やすように要求するメッセージを表示部１６に表示させてもよい。

［式選択ステップ］
以下、図２のステップＳ３０において実行される式選択ステップが説明される。式選択ステップは、ユーザがモデルを選択するステップと、CPU１１がユーザによって選択されたモデルに応じて、後述する推定ステップ（図２のステップＳ４０）において用いる基本式に決定するステップとを含む。ここで、基本式は、薬物動態パラメータx_iのベクトルx=(x₁,x₂,・・・,x_n)^Tの関数s(x)で表される。

ユーザがモデルを選択するステップは、上述したように、ユーザがモデル選択領域８４、８５（図４及び図５参照）に表示された３つのモデルの何れかを選択して「推定実行」ボタン３７、３８をクリックするステップである。

ＣＰＵ１１は、方法選択ステップ（図２のステップＳ２０）で最尤法が選択された場合において、正規分布モデル１が選択された場合、以下の（式１）を基本式とし、正規分布モデル２が選択された場合、以下の（式２）を基本式とし、対数正規分布モデルが選択された場合、以下の（式３）を基本式とする。

ここで、m：血中濃度測定回数、t_j：j回目の測定時刻(j=1,2,・・・,m)
c_j：時刻t_jにおける血中濃度測定値、c(t_j)=c(t_j,x)：時刻t_jにおける血中濃度理論値、σ-_c(但し、"-"は文字σのアッパーラインの意味である。)：jによらず同じ値をとる場合のlog(c_j)の標準偏差、ω_c：σ_j／c(t_j,x)で与えられる変動係数である。

なお、（式１）は、時刻t_jにおける血中濃度測定値c_jが正規分布N(c(t_j),σ_c ²)に従う場合、つまり誤差の分散が一定の場合の最尤法の基本式であり、（式２）は、血中濃度測定値c_jが正規分布N(c(t_j),{ω_c×c(t_j)}²)に従う場合、つまり誤差の分散が血中濃度に比例する場合の最尤法の基本式であり、（式３）は、上記血中濃度測定値c_jの対数値log(c_j)が正規分布N(log{c(t_j)},(σ-_c ²))(但し、"-"は文字σのアッパーラインの意味である。)に従う場合、つまり誤差の分散が一定の場合の最尤法の基本式である。

また、CPU１１は、方法選択ステップ（図２のステップＳ２０）でベイズ推定が選択された場合において、正規分布モデル１が選択された場合、以下の（式４）を基本式とし、正規分布モデル２が選択された場合、以下の（式５）を基本式とし、対数正規分布モデルが選択された場合、以下の（式６）を基本式とする。

ここで、m：血中濃度測定回数、n：薬物動態パラメータ数、t_j：j回目の測定時刻(j=1,2,・・・,m)、c_j：時刻t_jにおける血中濃度測定値、c(t_j)=c(t_j,x)：時刻t_jにおける血中濃度理論値、σ_c：jによらず同じ値をとる場合のc_jの標準偏差、σ-_c(但し、"-"は文字σのアッパーラインの意味である。)：jによらず同じ値をとる場合のlog(c_j)の標準偏差、μ_i：x_i(i=1,2,・・・,n)の平均、μ-_i(但し、"-"は文字μのアッパーラインの意味である。)：log(x_i)の平均、σ_i：x_iの標準偏差、σ-_i(但し、"-"は文字σのアッパーラインの意味である。)：log(x_i)の標準偏差、ω_c：σ_c／c(t_j,x)で与えられる変動係数、ω_i：σ_i／μ_iで与えられる変動係数、ｘ-_i(但し、"-"は文字xのアッパーラインの意味である。)：log(x_i)である。

なお、（式４）は、時刻t_jにおける血中濃度測定値c_jが正規分布N(c(t_j),σ_c ²)に従い且つ薬物動態パラメータx_iが正規分布N(μ_j,σ_i ²)に従う場合、つまり誤差の分散が一定の場合のベイズ推定の基本式であり、（式５）は、血中濃度測定値c_jが正規分布N(c(t_j),{ω_c×c(t_j)}²)に従い且つ薬物動態パラメータx_iが正規分布N(μ_j,(ω_i×μ_i)²)に従う場合、つまり誤差の分散が血中濃度に比例する場合のベイズ推定の基本式であり、（式６）は、血中濃度測定値c_jの対数値log(c_j)が正規分布N(log{c(t_j)},(σ-_c ²))(但し、"-"は文字σのアッパーラインの意味である。)に従い且つ薬物動態パラメータx_iの対数値log(x_i)が正規分布N(μ-_i,σ-_i ²)(但し、"-"は文字μ,σそれぞれのアッパーラインの意味である。)に従う場合のベイズ推定の基本式である。なお、（式６）は、血中濃度測定値c_jの変動が大きい場合に適用される。

［推定ステップ］
以下、図２のステップＳ４０において実行される推定ステップの詳細が、図７のフローチャートに基づいて説明される。推定ステップは、CPU１１が上記プログラムに従って以下に詳述する演算を実行することによって実現される。

推定ステップは、以下に詳述するように、方法選択ステップ（図２のステップＳ２０）において選択された方法（最尤法またはベイズ推定）における基本式（（式１）〜（式６））である薬物動態パラメータx_iのベクトルx=(x₁,x₂,・・・,x_n)^Tの関数s(x)に対して２階偏微分を実行し、当該２階偏微分によって得られた式から同一の薬物動態パラメータx_iによる２階偏導関数の項を削除した近似式に修正マルカート法を適用して関数s(x)を最小化し、最小化した際のベクトルxを求めることによって、薬物動態パラメータx_iを推定する。

上記のような、修正マルカート法による薬物動態パラメータx_iの推定は、図７のフローチャートに従った処理が実行されることによって実現される。つまり、推定ステップは、初期値準備ステップ（Ｓ２１０）、計算ステップ（Ｓ２２０）、変数変更ステップ（Ｓ２３０）、判定ステップ（Ｓ２４０）、及び決定ステップ（Ｓ２７０）を含む。

なお、図７のフローチャートにおいて、計算ステップ（Ｓ２２０）は、式選択ステップ（図２のステップＳ３０）において基本式とされた（式１）〜（式６）に応じて異なる行列による式が導出されるが、ステップＳ２２０以外においては、（式１）〜（式６）の何れが基本式とされても同処理である。よって、以下の説明では、式選択ステップ（図２のステップＳ３０）において（式１）が基本式とされたと仮定して図７のフローチャートの説明を行い、その後で、他の式（（式２）〜（式６））が基本式とされた場合の計算ステップ（Ｓ２２０）について説明する。

初期値準備ステップ（Ｓ２１０）において、CPU１１は、薬物動態パラメータx_iのベクトルxの初期値x⁽⁰⁾=(x₁ ⁽⁰⁾,x₂ ⁽⁰⁾,・・・,x_n ⁽⁰⁾)^Tを準備する。本実施形態では、以下に詳述するように、初期値x⁽⁰⁾の各成分(x₁ ⁽⁰⁾,x₂ ⁽⁰⁾,・・・,x_n ⁽⁰⁾)^T、つまり各薬物動態パラメータx_iは、ボックス・ミューラー法を用いて作成される。なお、初期値x⁽⁰⁾は、以下に詳述する方法以外の公知の方法によって準備されてもよい。

以下に、初期値x⁽⁰⁾を作成する手順について詳述する。最初に、独立な区間[0,1)の一様乱数を複数組用意する。例えば、統計数理研究所によって提供されている多種多様な物理乱数・算術乱数を用意する。次に、２組ずつの一様乱数からボックス・ミューラー法により、正規分布N(0,1)に従う正規乱数を作成する。なお、ボックス・ミューラー法は、u,vを独立な区間[0,1)の一様乱数として、以下の（式７）とする方法である。

このとき、x,yはN(0,1)）に従う独立な正規乱数となる。更に、薬物動態パラメータx_i(i=1,2,・・・,n)に対して、適当な平均値μ_xi及び変動係数ω_xiを考え、一様乱数列ｒ_iを用いて、以下の（式８）として初期値を作成する。

また、初期値準備ステップ（Ｓ２１０）において、CPU１１は、最小値が更新されたか否かを表す指標であるパラメータrの初期値を０とする。準備されたベクトルxの初期値x⁽⁰⁾及びパラメータrの初期値はRAM１２に格納される。

初期値準備ステップ（Ｓ２１０）の実行が完了すると、計算ステップ（Ｓ２２０）が実行される。計算ステップ（Ｓ２２０）において、CPU１１は、最初に、式選択ステップ（図２のステップＳ３０）において選択された基本式である薬物動態パラメータx_iのベクトルx=(x₁,x₂,・・・,x_n)^Tの関数s(x)、つまり（式１）に対して２階偏微分を実行する。以下に、（式１）に対して１階偏微分が実行された（式９）と２階偏微分が実行された（式１０）とが示される。

次に、CPU１１は、当該２階偏微分によって得られた式から同一の薬物動態パラメータx_iによる２階偏導関数の項を削除することによって、（式１１）に示される近似式を導出する。

次に、CPU１１は、近似式（式１１）を以下の（式１２）に適用する。ここで、（式１２）はニュートン法の式である。

すると、以下の（式１３）が導出される。そして、（式１３）を変形することにより（式１４）が導出される。

次に、CPU１１は、（式１４）から、ベクトルxのｉ（iは０以上の整数）番目の近似値であるx⁽ⁱ⁾=(x₁ ⁽ⁱ⁾,x₂ ⁽ⁱ⁾,・・・,x_n ⁽ⁱ⁾)^T、及びΔx=(Δx₁,Δx₂,・・・,Δx_n）^Tを含む行列による（式１５）を導出する。ここで、（式１５）においてA,Δx,及びb1は以下の（式１６）、（式１７）、（式１８）の通りである。また、（式１５）において、Iはm次の単位行列である。また、（式１４）、（式１５）において、λは正のパラメータで、2^r-1である。λを大きくすればするほど、Δxは最急降下方向に近づく

次に、CPU１１は、行列による（式１５）においてΔxに対するｎ元連立一次方程式を解き、x'=x⁽ⁱ⁾+Δxとおく。

変数変更ステップ（Ｓ２３０）において、CPU１１は、計算ステップ（Ｓ２２０）においてｎ元連立一次方程式を解いた結果、s(x')<s(x⁽ⁱ⁾)であることを条件として、パラメータrをmin(0,r-1)とし且つiに１を加えてx⁽ⁱ⁾=x'とする。一方、s(x')>=s(x⁽ⁱ⁾)であることを条件として、パラメータrに１を加える。つまり、変数変更ステップ（Ｓ２３０）では、パラメータr及び偏微分の階数iの値を調整する。

変数変更ステップ（Ｓ２３０）の実行後に、判定ステップ（Ｓ２４０）が実行される。判定ステップ（Ｓ２４０）において、CPU１１は、第１条件及び第２条件を具備するか否かを判断する。ここで、第１条件は、i>=Mまたはr>=Nであることである。このとき、M,Nは予め設定された正の整数定数である。また、第２条件は、|s(x⁽ⁱ⁺¹⁾)-s(x⁽ⁱ⁾)|<ε且つ|Δx|<δ且つr=0であることである。このとき、ε,δは予め設定された小さな実数正数である。

CPU１１は、第１条件を具備することを条件として、s(x)が収束しないと判定する（Ｓ２５０）。このとき、薬物動態パラメータx_iは推定されることなく、推定ステップは終了する。なお、CPU１１は、その旨のメッセージを表示部１６に表示させてもよい。

一方、CPU１１は、第２条件を具備することを条件として、s(x)が収束すると判定する（Ｓ２６０）。このとき、決定ステップ（Ｓ２７０）が実行される。決定ステップ（Ｓ２７０）において、CPU１１は、s(x)が収束すると判定された時点におけるベクトルx⁽ⁱ⁾の各成分を薬物動態パラメータx_iとする。つまり、決定ステップでは、s(x)の最小値を与える薬物動態パラメータx_iが選定される。

また、一方、CPU１１は、第１条件及び第２条件の何れも具備しないことを条件として、再び計算ステップ（Ｓ２２０）を実行する。この計算ステップ（Ｓ２２０）においては、以前に実行された計算ステップ（Ｓ２２０）及び変数変更ステップ（Ｓ２３０）によってパラメータrや偏微分の階数iが変更されている。そのため、次回の判定ステップ（Ｓ２４０）においては、前回と異なる判定結果となり得る。

以上のように、修正マルカート法では、初期値が設定された上で（Ｓ２１０）、判定ステップ（Ｓ２４０）においてs(x)が収束する（Ｓ２６０）或いは収束しない（Ｓ２５０）と判定されるまで、計算ステップ（Ｓ２２０）及び変数変更ステップ（Ｓ２３０）が繰り返し実行される。つまり、推定ステップでは、計算ステップ（Ｓ２２０）及び変数変更ステップ（Ｓ２３０）が反復されることによって、ｎ元連立非線形方程式が解かれる。

以下、式選択ステップ（図２のステップＳ３０）において（式２）が基本式とされた場合の計算ステップ（Ｓ２２０）について説明する。計算ステップ（Ｓ２２０）において、CPU１１は、最初に、式選択ステップ（図２のステップＳ３０）において選択された基本式である（式２）に対して２階偏微分を実行する。次に、CPU１１は、当該２階偏微分によって得られた式から同一の薬物動態パラメータx_iによる２階偏導関数の項を削除することによって、近似式を導出する。以下に、（式２）に対して１階偏微分が実行された（式１９）と２階偏微分が実行され且つ２階偏導関数の項を削除された（式２０）とが示される。

次に、CPU１１は、近似式（式２０）を（式１２）に適用する。これにより、以下の（式２１）が導出される。そして、（式２１）を変形することにより（式２２）が導出される。

次に、CPU１１は、（式２２）から、ベクトルxのi（iは０以上の整数）番目の近似値であるx⁽ⁱ⁾=(x₁ ⁽ⁱ⁾,x₂ ⁽ⁱ⁾,・・・,x_n ⁽ⁱ⁾)^T、及びΔx=(Δx₁,Δx₂,・・・,Δx_n）^Tを含む行列による（式２３）を導出する。ここで、（式２３）においてA,D,Δx,及びb1は以下の（式２４）、（式２５）、（式２６）、（式２７）の通りである。また、（式２３）において、Iはm次の単位行列とする。

次に、CPU１１は、行列による（式２３）においてΔxに対するｎ元連立一次方程式を解き、x'=x⁽ⁱ⁾+Δxとおく。

以下、式選択ステップ（図２のステップＳ３０）において（式３）が基本式とされた場合の計算ステップ（Ｓ２２０）について説明する。

（式３）が基本式とされた場合は、x,c_j,c(t_j,x)をそれぞれlog(x),log(c_j),log(c(t_j,x))とすることにより、（式１）が基本式とされた場合に対する計算方法をそのまま使用することが可能である。

以下、式選択ステップ（図２のステップＳ３０）において（式４）が基本式とされた場合の計算ステップ（Ｓ２２０）について説明する。計算ステップ（Ｓ２２０）において、CPU１１は、最初に、式選択ステップ（図２のステップＳ３０）において選択された基本式である（式４）に対して２階偏微分を実行する。次に、CPU１１は、当該２階偏微分によって得られた式から同一の薬物動態パラメータx_iによる２階偏導関数の項を削除することによって、近似式を導出する。以下に、（式４）に対して１階偏微分が実行された（式２８）と２階偏微分が実行され且つ２階偏導関数の項を削除された（式２９）とが示される。

次に、CPU１１は、近似式（式２９）を（式１２）に適用する。これにより、以下の（式３０）が導出される。そして、（式３０）を変形することにより（式３１）が導出される。

次に、CPU１１は、（式３１）から、ベクトルxのi（iは０以上の整数）番目の近似値であるx⁽ⁱ⁾=(x₁ ⁽ⁱ⁾,x₂ ⁽ⁱ⁾,・・・,x_n ⁽ⁱ⁾)^T、及びΔx=(Δx₁,Δx₂,・・・,Δx_n）^Tを含む行列による（式３２）を導出する。ここで、（式３２）においてA,B,Δx,b1,及びb3は以下の（式３３）、（式３４）、（式３５）、（式３６）、（式３７）の通りである。また、（式３２）において、Iはm次の単位行列とする。

次に、CPU１１は、行列による（式３２）においてΔxに対するｎ元連立一次方程式を解き、x'=x⁽ⁱ⁾+Δxとおく。

以下、式選択ステップ（図２のステップＳ３０）において（式５）が基本式とされた場合の計算ステップ（Ｓ２２０）について説明する。計算ステップ（Ｓ２２０）において、CPU１１は、最初に、式選択ステップ（図２のステップＳ３０）において選択された基本式である（式５）に対して２階偏微分を実行する。次に、CPU１１は、当該２階偏微分によって得られた式から同一の薬物動態パラメータx_iによる２階偏導関数の項を削除することによって、近似式を導出する。以下に、（式４）に対して１階偏微分が実行された（式３８）と２階偏微分が実行され且つ２階偏導関数の項を削除された（式３９）とが示される。

次に、CPU１１は、近似式（式３９）を（式１２）に適用する。これにより、以下の（式４０）が導出される。そして、（式４０）を変形することにより（式４１）が導出される。

次に、CPU１１は、（式４１）から、ベクトルxのi（iは０以上の整数）番目の近似値であるx⁽ⁱ⁾=(x₁ ⁽ⁱ⁾,x₂ ⁽ⁱ⁾,・・・,x_n ⁽ⁱ⁾)^T、及びΔx=(Δx₁,Δx₂,・・・,Δx_n）^Tを含む行列による（式４２）を導出する。ここで、（式３２）においてA,D,B,Δx,b2,及びb4は以下の（式４３）、（式４４）、（式４５）、（式４６）、（式４７）、（式４８）の通りである。また、（式４２）において、Iはm次の単位行列とする。

次に、CPU１１は、行列による（式４２）においてΔxに対するｎ元連立一次方程式を解き、x'=x⁽ⁱ⁾+Δxとおく。

以下、式選択ステップ（図２のステップＳ３０）において（式６）が基本式とされた場合の計算ステップ（Ｓ２２０）について説明する。

（式６）が基本式とされた場合は、x,c_j,c(t_j,x),σ_c,x_i,μ_i,σ_iをそれぞれlog(x),log(c_j),log{c(t_j,x)},log(σ_c),log(x_i),log(μ_i),log(σ_i)とすることにより、（式４）が基本式とされた場合に対する計算方法をそのまま使用することが可能である。

［実施形態の作用効果］
本実施形態は、方法選択ステップ（Ｓ２０）を含む。これにより、血中濃度の測定回数が薬物動態パラメータx_iの個数以上の場合、過去のデータから得られた母集団パラメータを用いることなく、最尤法によって薬物動態パラメータx_iを推定することができる。一方、血中濃度の測定回数が薬物動態パラメータx_iの個数未満の場合、母集団パラメータを用いて、ベイズ推定によって薬物動態パラメータx_iを推定することができる。

また、本実施形態では、修正マルカート法が用いられている。このため、ニュートン法よりも収束性を改善することができる。また、修正マルカート法は、ガウス・ニュートン法とは異なり、２乗和以外の項がある関数に適用した場合でも、収束性を改善することができる。

また、本願の出願人は、最尤法及びベイズ推定の基本式を２階偏微分して得られた式から２階偏導関数の項を削除した近似式を用いることによって、関数s(x)が発散する可能性を低くすることができることを発見した。

以上より、本実施形態によると、推定する薬物動態パラメータx_iの数にかかわらず、当該薬物動態パラメータx_iを、発散の可能性を低く抑えつつ推定することができる。

また、本実施形態によれば、最尤法が選択された場合において、血中濃度測定値が正規分布に従う場合（（式１）及び（式２））と対数正規分布に従う場合（式３）とに応じて式を使い分けることができる。また、血中濃度測定値が正規分布に従う場合には更に、血中濃度測定値の誤差の分散が一定の場合（式１）と、血中濃度測定値の誤差の分散が血中濃度に比例する場合（式２）とに応じて式を使い分けることができる。

また、本実施形態によれば、ベイズ推定が選択された場合において、血中濃度の測定において混入した誤差が正規分布に従う場合（（式４）及び（式５））と対数正規分布に従う場合（式６）とに応じて式を使い分けることができる。また、上記誤差が正規分布に従う場合には更に、血中濃度測定値の誤差の分散が一定の場合（式４）と、血中濃度測定値の誤差の分散が血中濃度に比例する場合（式５）とに応じて式を使い分けることができる。

また、本実施形態のように、行列による式を用いることによって、計算が容易となる。

また、本実施形態のようにボックス・ミューラー法を用いることにより、作成された初期値x⁽⁰⁾の各成分を正規分布に従う独立な正規乱数とすることができる。

［変形例１］
上述の実施形態では、CPU１１が、未知の薬物動態パラメータの数と測定回数とを比較し、その比較の結果によって、最尤法とベイズ推定とを選択していた。しかし、最尤法とベイズ推定との選択は、ユーザによって行われてもよい。例えば、図８に示されるように、データ入力画面３１に、最尤法における正規分布モデル１、正規分布モデル２、及び対数正規分布モデル、並びに、ベイズ推定における正規分布モデル１、正規分布モデル２、及び対数正規分布モデルを選択するためのモデル選択領域４５を設けてもよい。これにより、ユーザは、データ入力画面３１において、推定方法及び当該推定方法におけるモデルを選択することができる。この場合、測定回数が未知の薬物動態パラメータの個数よりも多い場合であっても、ベイズ推定を選択することができる。

［変形例２］
上述の実施形態では、各推定方法（最尤法及びベイズ推定）におけるモデル選択は、ユーザがモデル選択画面３５において選択することによって行われていた。しかし、上記モデル選択は、CPU１１によって自動的に行われてもよい。例えば、CPU１１は、データ入力画面３１においてユーザに入力された測定データや、大容量記憶部１３に格納された母集団パラメータなどに基づいて、最適と判断されるモデルを選択してもよい。具体的には、CPU１１は、最尤法が選択されている場合において測定データが概ね正規分布に従い且つ誤差が略一定の場合に、正規分布モデル１を選択する。

［変形例３］
上述の実施形態では、初期値準備ステップ（図７のステップＳ２１０）において、CPU１１は、１つの初期値x⁽⁰⁾を準備したが、複数の初期値x⁽⁰⁾を準備してもよい。この際、各初期値x⁽⁰⁾は、上述の実施形態と同様に、ボックス・ミューラー法を用いて作成されてもよい。

変形例３では、推定ステップは、図９のフローチャートに従って実行される。つまり、複数の初期値x⁽⁰⁾が準備された場合、全ての初期値x⁽⁰⁾の各々について図７のステップＳ２２０〜Ｓ２７０が実行される。そして、CPU１１は、各初期値x⁽⁰⁾に基づく演算結果について、s(x)が収束するか否かを判断し、s(x)が収束する場合にはその時点におけるベクトルx⁽ⁱ⁾を求める。全ての初期値x⁽⁰⁾について、s(x)が収束するか否かが判断されると（Ｓ４１０：Ｙｅｓ）、CPU１１は、収束したベクトルx⁽ⁱ⁾の値を比較し、最も小さいベクトルx⁽ⁱ⁾の各成分を薬物動態パラメータｘ_ｉと決定する（Ｓ４２０）。つまり、CPU１１は、複数の初期値x⁽⁰⁾を準備した上で、初期値x⁽⁰⁾毎にs(x)の値を求め、最も小さい値に収束したs(x)の各成分を薬物動態パラメータx_iと決定する。

薬物動態パラメータx_iの推定過程において、関数s(x)が発散してしまうと、薬物動態パラメータx_iを推定することができない。また、関数s(x)が収束した場合でも、当該収束が本来収束して欲しい解以外の局所解に収束した場合には、適切な薬物動態パラメータx_iを推定することができていない。そこで、変形例３のように、初期値x⁽⁰⁾を複数個準備して、複数の関数s(x)の解を求め、当該解が最も小さく収束したときのx⁽ⁱ⁾の各成分を薬物動態パラメータx_iと決定することで、局所解ではない適切な薬物動態パラメータx_iを推定することができる。

また、初期値準備ステップにおいて、ボックス・ミューラー法を用いた場合、上述の実施形態と同様に、作成された初期値x⁽⁰⁾の各成分を正規分布に従う独立な正規乱数とすることができる。その結果、初期値x⁽⁰⁾を複数個準備して、複数の関数s(x)の解を求めた際に、当該解の少なくとも一部が適切に収束する可能性を高くすることができる。

１０・・・電子計算機
１１・・・CPU
３１・・・データ入力画面
３５・・・モデル選択画面
３６・・・モデル選択画面
４１・・・血中濃度推移画面

本発明は、患者から採取された血液における薬物の血中濃度から薬物動態パラメータを推定する方法、及び当該方法を実行するための薬物動態パラメータの推定プログラムに関する。

個々の患者に適した薬物の投与設計を行う治療薬物モニタリング(Therapeutic Drug Monitoring,TDM)が知られている。TDMでは、患者の血液における薬物の血中濃度を測定して解析した結果に基づいて、最適な薬用量、投与法が設定される。また、TDMでは、測定された血中濃度から薬物動態パラメータが推定される。

上述したようなTDMを実施するためのシステムの一例として、ベイズ推定を用いて予測血中濃度を計算する薬物投与解析システムが、特許文献１に開示されている。また、薬物動態パラメータの推定には、数値計算法の一つであるニュートン法やガウス・ニュートン法が用いられることが多い。

特開平７−１２４１２５号公報

ベイズ推定は、一人の患者に対する血中濃度の測定回数が少なく、測定回数が推定するべき薬物動態パラメータの数以下である場合に、過去のデータから得られた母集団パラメータの値を使用して薬物動態パラメータを推定する統計手法である。しかし、一人の患者に対する血中濃度の測定回数が多いときには、母集団パラメータを用いずとも薬物動態パラメータを推定することができるため、ベイズ推定以外の方法を用いる方がよい場合がある。

また、ニュートン法は、方程式f(x)=0におけるxの値の近似値を反復することによって求める手法である。しかし、ニュートン法では、方程式に最初に代入する初期値x⁽⁰⁾によっては、xの値が収束せずに発散してしまう。そして、ニュートン法では、xの値を収束させるための初期値x⁽⁰⁾の収束域が狭いため、xの値を発散させない適切な初期値x⁽⁰⁾を選択することが困難であるという問題がある。

また、f(x)が２乗和の形の関数の場合、ニュートン法よりも収束性を改善させるための手法として、ガウス・ニュートン法がある。しかし、ガウス・ニュートン法を、２乗和以外の項がある関数に適用した場合、収束性の改善は期待できない。

本発明は、上記事情に鑑みてなされたものであり、その目的は、推定する薬物動態パラメータの数にかかわらず、当該薬物動態パラメータを、発散の可能性を低く抑えつつ推定することができる方法、及び当該方法を実行するための薬物動態パラメータの推定プログラムを提供することにある。

(1) 本発明に係る薬物動態パラメータの推定方法は、患者から採取された血液における薬物の血中濃度から薬物動態パラメータx_i(i=1,2,・・・,n)を推定する方法であって、薬物動態パラメータx_iを推定する方法として最尤法またはベイズ推定の何れかが選択される方法選択ステップと、上記方法選択ステップにおいて選択された方法における基本式である上記薬物動態パラメータx_iのベクトルx=(x₁,x₂,・・・,x_n)^Tの関数s(x)に対して２階偏微分を実行し、当該２階偏微分によって得られた式から血中濃度に対する薬物動態パラメータx_iによる２階偏導関数の項を削除した近似式に修正マルカート法を適用して上記関数s(x)を最小化し、最小化した際の上記ベクトルxを求めることによって、上記薬物動態パラメータx_iを推定する推定ステップと、を含む。

本方法は、方法選択ステップを含む。これにより、血中濃度の測定回数が推定すべき薬物動態パラメータx_iの個数よりも多い場合、過去のデータから得られた母集団パラメータを用いることなく、最尤法によって薬物動態パラメータx_iを推定することができる。一方、血中濃度の測定回数が推定すべき薬物動態パラメータx_iの個数以下の場合、母集団パラメータを用いて、ベイズ推定によって薬物動態パラメータx_iを推定することができる。

また、本方法では、修正マルカート法が用いられている。このため、ニュートン法よりも収束性を改善することができる。また、修正マルカート法は、ガウス・ニュートン法とは異なり、２乗和以外の項がある関数に適用した場合でも、収束性を改善することができる。

また、本願の出願人は、最尤法及びベイズ推定の基本式を２階偏微分して得られた式から２階偏導関数の項を削除した近似式を修正マルカート法を適用することによって、関数s(x)が発散する可能性を低くすることができることを発見した。

(2) 本発明に係る薬物動態パラメータの推定方法において、最尤法における上記関数s(x)は、時刻t_jにおける血中濃度測定値c_jが正規分布N(c(t_j),σ_c ²)に従う場合、以下の（式１）であり、上記血中濃度測定値c_jが正規分布N(c(t_j),{ω_c×c(t_j)}²)に従う場合、以下の（式２）であり、上記血中濃度測定値c_jの対数値log(c_j)が正規分布N(log{c(t_j)},(σ-_c ²))(但し、"-"は文字σのアッパーラインの意味である。)に従う場合、以下の（式３）である。

以上において、
m：血中濃度測定回数
t_j：j回目の測定時刻(j=1,2,・・・,m)
c_j：時刻t_jにおける血中濃度測定値
c(t_j)=c(t_j,x)：時刻t_jにおける血中濃度理論値
σ-_c(但し、"-"は文字σのアッパーラインの意味である。)：jによらず同じ値をとる場合のlog(c_j)の標準偏差
ω_c：σ_j／c(t_j,x)で与えられる変動係数
である。上記方法選択ステップにおいて最尤法が選択された場合に、当該３つの式の何れかが選択される式選択ステップを更に含む。

従来、最尤法による薬物動態パラメータの推定において、関数s(x)として（式１）が用いられていた。しかしながら、（式１）を用いた場合、推定した薬物動態パラメータの精度の点に問題があり、精度の向上が求められていた。この問題の解決のために、関数s(x)として（式２）を用いることが考えられる。しかしながら、（式２）を用いた場合、xの値の収束性が良くないという問題がある。つまり、推定過程においてxの値が発散してしまうことがある。

そこで、本願発明者は上記（１）のような薬物パラメータの推定方法を考え出した。当該推定方法によると、（式２）を用いた場合でも、xの値の収束性を良好にすることができる。その結果、本方法のように（式１）及び（式２）を選択的に使用することができ、しかも何れの式を使用した場合であっても、薬物動態パラメータを正確に推定することができる。

また、本方法によれば、血中濃度測定値が正規分布に従う場合（（式１）及び（式２））と対数正規分布に従う場合（式３）とに応じて式を使い分けることができる。また、血中濃度測定値が正規分布に従う場合には更に、血中濃度測定値の誤差の分散が一定の場合（式１）と、血中濃度測定値の誤差の分散が血中濃度に比例する場合（式２）とに応じて式を使い分けることができる。

(3) 本発明に係る薬物動態パラメータの推定方法において、ベイズ推定における上記関数s(x)は、時刻t_jにおける血中濃度測定値c_jが正規分布N(c(t_j),σ_c ²)に従い且つ薬物動態パラメータx_iが正規分布N(μ_i,σ_i ²)に従う場合、以下の（式４）であり、上記血中濃度測定値c_jが正規分布N(c(t_j),{ω_c×c(t_j)}²)に従い且つ薬物動態パラメータx_iが正規分布N(μ_i,(ω_i×μ_i)²)に従う場合、以下の（式５）であり、上記血中濃度測定値c_jの対数値log(c_j)が正規分布N(log{c(t_j)},(σ-_c ²))(但し、"-"は文字σのアッパーラインの意味である。)に従い且つ薬物動態パラメータx_iの対数値log(x_i)が正規分布N(μ-_i,σ-_i ²)(但し、"-"は文字μ,σそれぞれのアッパーラインの意味である。)に従う場合、以下の（式６）である。

以上において、
m：血中濃度測定回数
n：薬物動態パラメータ数
t_j：j回目の測定時刻(j=1,2,・・・,m)
c_j：時刻t_jにおける血中濃度測定値
c(t_j)=c(t_j,x)：時刻t_jにおける血中濃度理論値
σ_c：jによらず同じ値をとる場合のc_jの標準偏差
σ-_c(但し、"-"は文字σのアッパーラインの意味である。)：jによらず同じ値をとる場合のlog(c_j)の標準偏差
μ_i：x_i(i=1,2,・・・,n)の平均
μ-_i(但し、"-"は文字μのアッパーラインの意味である。)：log(x_i)の平均
σ_i：x_iの標準偏差
σ-_i(但し、"-"は文字σのアッパーラインの意味である。)：log(x_i)の標準偏差
ω_c：σ_c／c(t_j,x)で与えられる変動係数
ω_i：σ_i／μ_iで与えられる変動係数
ｘ-_i(但し、"-"は文字xのアッパーラインの意味である。)：log(x_i)
である。上記方法選択ステップにおいてベイズ推定が選択された場合に、当該３つの式の何れかが選択される式選択ステップを更に含む。

従来、ベイズ推定による薬物動態パラメータの推定において、関数s(x)として（式４）が用いられていた。しかしながら、（式４）を用いた場合、推定した薬物動態パラメータの精度の点に問題があり、精度の向上が求められていた。この問題の解決のために、関数s(x)として（式５）を用いることが考えられる。しかしながら、（式５）を用いた場合、xの値の収束性が良くないという問題がある。つまり、推定過程においてxの値が発散してしまうことがある。

そこで、本願発明者は上記（１）のような薬物パラメータの推定方法を考え出した。当該推定方法によると、（式４）を用いた場合でも、xの値の収束性を良好にすることができる。その結果、本方法のように（式４）及び（式５）を選択的に使用することができ、しかも何れの式を使用した場合であっても、薬物動態パラメータを正確に推定することができる。

また、本方法によれば、血中濃度の測定において混入した誤差及び薬物動態パラメータが正規分布に従う場合（（式４）及び（式５））と対数正規分布に従う場合（式６）とに応じて式を使い分けることができる。また、上記誤差が正規分布に従う場合には更に、血中濃度測定値の誤差の分散が一定の場合（式４）と、血中濃度測定値の誤差の分散が血中濃度に比例する場合（式５）とに応じて式を使い分けることができる。

(4) 本発明に係る薬物動態パラメータの推定方法において、上記推定ステップは、上記ベクトルxの初期値x⁽⁰⁾=(x₁ ⁽⁰⁾,x₂ ⁽⁰⁾,・・・,x_n ⁽⁰⁾)^Tを準備し、最小値が更新されたか否かを表す指標であるrの初期値をr=0と準備する初期値準備ステップと、上記式選択ステップにおいて選択された式から、上記ベクトルxのi（iは０以上の整数）番目の近似値であるx⁽ⁱ⁾=(x₁ ⁽ⁱ⁾,x₂ ⁽ⁱ⁾,・・・,x_n ⁽ⁱ⁾)^T、及びΔx=(Δx₁,Δx₂,・・・,Δx_n）^Tを含む行列による式を導出し、当該行列による式において上記Δxに対するn元連立一次方程式を解き、x'=x⁽ⁱ⁾+Δxとおく計算ステップと、上記計算ステップにおいて上記n元連立一次方程式を解いた結果、s(x')<s(x⁽ⁱ⁾)であることを条件として、rをmin(0,r-1)とし且つiに１を加えてx⁽ⁱ⁾=x'とし、s(x')>=s(x⁽ⁱ⁾)であることを条件として、rに１を加える変数変更ステップと、上記変数変更ステップの実行後に、i>=Mまたはr>=N（但し、M,Nは予め設定された正の整数定数）である第１条件を具備することを条件として、s(x)が収束しないと判定し、|s(x⁽ⁱ⁺¹⁾)-s(x⁽ⁱ⁾)|<ε且つ|Δx|<δ且つr=0（但し、ε,δは予め設定された実数正数）である第２条件を具備することを条件として、s(x)が収束すると判定する判定ステップと、上記判定ステップにおいて、s(x)が収束すると判定されたことを条件として、そのときのx⁽ⁱ⁾の各成分を薬物動態パラメータx_iとする決定ステップと、を含む。上記計算ステップは、上記初期値準備ステップの実行が完了したこと、または、上記判定ステップにおいて上記第１条件及び上記第２条件を具備しないことを条件として実行される。

本方法のように、行列による式を用いることによって、計算が容易となる。

(5) 本発明に係る薬物動態パラメータの推定方法において、上記推定ステップは、上記初期値準備ステップにおいて、上記ベクトルxの初期値x⁽⁰⁾を複数個準備する。上記決定ステップは、複数個の初期値x⁽⁰⁾の各々について上記計算ステップ、上記変数変更ステップ、及び上記判定ステップが実行された結果、s(x)が最も小さく収束したときのx⁽ⁱ⁾の各成分を薬物動態パラメータx_iと決定する。

薬物動態パラメータx_iの推定過程において、xが発散してしまうと、薬物動態パラメータx_iを推定することができない。また、xが収束した場合でも、当該収束がs(x)を最小にする解以外の局所解に収束した場合には、適切な薬物動態パラメータx_iを推定することができていない。そこで、本方法のように、初期値x⁽⁰⁾を複数個準備して、複数の関数s(x)の解を求め、当該解が最も小さく収束したときのx⁽ⁱ⁾の各成分を薬物動態パラメータx_iと決定することで、局所解ではない適切な薬物動態パラメータx_iを推定することができる。

(6) 本発明に係る薬物動態パラメータの推定方法では、上記初期値準備ステップにおいて、上記ベクトルxの初期値x⁽⁰⁾は、ボックス・ミューラー法によって作成される。

ボックス・ミューラー法を用いることにより、作成された初期値x⁽⁰⁾の各成分を正規分布に従う独立な正規乱数とすることができる。特に、初期値x⁽⁰⁾を複数個準備した場合、複数の解xを求めた際に、当該解の少なくとも一部が適切に収束する可能性を高くすることができる。

(7) 本発明は、患者から採取された血液における薬物の血中濃度から薬物動態パラメータx_i(i=1,2,・・・,n)を推定するためのプログラムであって、薬物動態パラメータx_iを推定する方法として最尤法またはベイズ推定の何れかが選択される方法選択ステップと、上記方法選択ステップにおいて選択された方法における基本式である上記薬物動態パラメータx_iのベクトルx=(x₁,x₂,・・・,x_n)^Tの関数s(x)に対して２階偏微分を実行し、当該２階偏微分によって得られた式から血中濃度に対する薬物動態パラメータx_iによる２階偏導関数の項を削除した近似式に修正マルカート法を適用して上記関数s(x)を最小化し、最小化した際の上記ベクトルxを求めることによって、上記薬物動態パラメータx_iを推定する推定ステップと、を実行させる薬物動態パラメータの推定プログラムとして捉えることもできる。

(8) 本発明に係る薬物動態パラメータの推定プログラムにおいて、最尤法における上記関数s(x)は、時刻t_jにおける血中濃度測定値c_jが正規分布N(c(t_j),σ_c ²)に従う場合、以下の（式１）であり、上記血中濃度測定値c_jが正規分布N(c(t_j),{ω_c×c(t_j)}²)に従う場合、以下の（式２）であり、上記血中濃度測定値c_jの対数値log(c_j)が正規分布N(log{c(t_j)},(σ-_c ²))(但し、"-"は文字σのアッパーラインの意味である。)に従う場合、以下の（式３）である。

以上において、
m：血中濃度測定回数
t_j：j回目の測定時刻(j=1,2,・・・,m)
c_j：時刻t_jにおける血中濃度測定値
c(t_j)=c(t_j,x)：時刻t_jにおける血中濃度理論値
σ-_c(但し、"-"は文字σのアッパーラインの意味である。)：jによらず同じ値をとる場合のlog(c_j)の標準偏差
ω_c：σ_j／c(t_j,x)で与えられる変動係数
である。上記方法選択ステップにおいて最尤法が選択された場合に、当該３つの式の何れかが選択される式選択ステップを更に実行させる。

(9) 本発明に係る薬物動態パラメータの推定プログラムにおいて、ベイズ推定における上記関数s(x)は、時刻t_jにおける血中濃度測定値c_jが正規分布N(c(t_j),σ_c ²)に従い且つ薬物動態パラメータx_iが正規分布N(μ_i,σ_i ²)に従う場合、以下の（式４）であり、上記血中濃度測定値c_jが正規分布N(c(t_j),{ω_c×c(t_j)}²)に従い且つ薬物動態パラメータx_iが正規分布N(μ_i,(ω_i×μ_i)²)に従う場合、以下の（式５）であり、上記血中濃度測定値c_jの対数値log(c_j)が正規分布N(log{c(t_j)},(σ-_c ²))(但し、"-"は文字σのアッパーラインの意味である。)に従い且つ薬物動態パラメータx_iの対数値log(x_i)が正規分布N(μ-_i,σ-_i ²)(但し、"-"は文字μ,σそれぞれのアッパーラインの意味である。)に従う場合、以下の（式６）である。

以上において、
m：血中濃度測定回数
n：薬物動態パラメータ数
t_j：j回目の測定時刻(j=1,2,・・・,m)
c_j：時刻t_jにおける血中濃度測定値
c(t_j)=c(t_j,x)：時刻t_jにおける血中濃度理論値
σ_c：jによらず同じ値をとる場合のc_jの標準偏差
σ-_c(但し、"-"は文字σのアッパーラインの意味である。)：jによらず同じ値をとる場合のlog(c_j)の標準偏差
μ_i：x_i(i=1,2,・・・,n)の平均
μ-_i(但し、"-"は文字μのアッパーラインの意味である。)：log(x_i)の平均
σ_i：x_iの標準偏差
σ-_i(但し、"-"は文字σのアッパーラインの意味である。)：log(x_i)の標準偏差
ω_c：σ_c／c(t_j,x)で与えられる変動係数
ω_i：σ_i／μ_iで与えられる変動係数
ｘ-_i(但し、"-"は文字xのアッパーラインの意味である。)：log(x_i)
である。上記方法選択ステップにおいてベイズ推定が選択された場合に、当該３つの式の何れかが選択される式選択ステップを更に実行させる。

(10) 本発明に係る薬物動態パラメータの推定プログラムにおいて、上記推定ステップは、上記ベクトルxの初期値x⁽⁰⁾=(x₁ ⁽⁰⁾,x₂ ⁽⁰⁾,・・・,x_n ⁽⁰⁾)^Tを準備し、最小値が更新されたか否かを表す指標であるrの初期値をr=0と準備する初期値準備ステップと、上記式選択ステップにおいて選択された式から、上記ベクトルxのi（iは０以上の整数）番目の近似値であるx⁽ⁱ⁾=(x₁ ⁽ⁱ⁾,x₂ ⁽ⁱ⁾,・・・,x_n ⁽ⁱ⁾)^T、及びΔx=(Δx₁,Δx₂,・・・,Δx_n）^Tを含む行列による式を導出し、当該行列による式において上記Δxに対するn元連立一次方程式を解き、x'=x⁽ⁱ⁾+Δxとおく計算ステップと、上記計算ステップにおいて上記n元連立一次方程式を解いた結果、s(x')<s(x⁽ⁱ⁾)であることを条件として、rをmin(0,r-1)とし且つiに１を加えてx⁽ⁱ⁾=x'とし、s(x')>=s(x⁽ⁱ⁾)であることを条件として、rに１を加える変数変更ステップと、上記変数変更ステップの実行後に、i>=Mまたはr>=N（但し、M,Nは予め設定された正の整数定数）である第１条件を具備することを条件として、s(x)が収束しないと判定し、|s(x⁽ⁱ⁺¹⁾)-s(x⁽ⁱ⁾)|<ε且つ|Δx|<δ且つr=0（但し、ε,δは予め設定された実数正数）である第２条件を具備することを条件として、s(x)が収束すると判定する判定ステップと、上記判定ステップにおいて、s(x)が収束すると判定されたことを条件として、そのときのx⁽ⁱ⁾の各成分を薬物動態パラメータx_iとする決定ステップと、を実行させる。上記計算ステップは、上記初期値準備ステップの実行が完了したこと、または、上記判定ステップにおいて上記第１条件及び上記第２条件を具備しないことを条件として実行される。

(11) 本発明に係る薬物動態パラメータの推定プログラムにおいて、上記推定ステップは、上記初期値準備ステップにおいて、上記ベクトルxの初期値x⁽⁰⁾を複数個準備する。上記決定ステップは、複数個の初期値x⁽⁰⁾の各々について上記計算ステップ、上記変数変更ステップ、及び上記判定ステップが実行された結果、s(x)が最も小さく収束したときのx⁽ⁱ⁾の各成分を薬物動態パラメータx_iと決定する。

(12) 本発明に係る薬物動態パラメータの推定プログラムでは、上記初期値準備ステップにおいて、上記ベクトルxの初期値x⁽⁰⁾は、ボックス・ミューラー法によって作成される。

本発明に係る薬物動態パラメータの推定方法及び薬物動態パラメータの推定プログラムによると、推定する薬物動態パラメータx_iの数にかかわらず、当該薬物動態パラメータx_iを、発散の可能性を低く抑えつつ推定することができる。

図１は、電子計算機１０のブロック図である。 図２は、薬物動態パラメータを推定するためのシステムの処理について説明するためのフローチャートである。 図３は、表示部１６に表示されたデータ入力画面３１を示す図である。 図４は、表示部１６に表示されたモデル選択画面３５を示す図である。 図５は、表示部１６に表示されたモデル選択画面３６を示す図である。 図６は、表示部１６に表示された血中濃度推移画面４１を示す図である。 図７は、推定ステップの処理について説明するためのフローチャートである。 図８は、変形例１における表示部１６に表示されたデータ入力画面３１を示す図である。 図９は、変形例３における推定ステップの処理について説明するためのフローチャートである。

以下、適宜図面を参照して本発明の実施形態について説明する。なお、以下に説明される実施形態は本発明の一例にすぎず、本発明の要旨を変更しない範囲で、本発明の実施形態を適宜変更できることは言うまでもない。

［薬物動態パラメータを推定するためのシステムの構成］
本発明に係る薬物動態パラメータの推定方法によって推定される薬物動態パラメータは、コンパートメントモデルに追随して決まる。ここで、コンパートメントモデルは、薬物動態に関する体内のモデルである。コンパートメントモデルは、対象となる薬剤によって決まる。

コンパートメントモデルとしては、例えば点滴静脈注射のコンパートメントモデル、経口・筋肉注射のコンパートメントモデル、及び点滴のコンパートメントモデルなどがある。各コンパートメントモデルは、定義された複数の薬物動態パラメータと、当該複数の薬物動態パラメータ間で成立する少なくとも１つの関係式よりなる。つまり、コンパートメントモデルが決定されると、使用する未知の薬物動態パラメータのパラメータ数が決まる。なお、コンパートメントモデルは、上述した３種類に限らず、他にも存在する。また、同じ種類のコンパートメントモデルであっても、式の変換により別の薬物動態パラメータが定義されることもある。

以下、点滴のコンパートメントモデルの１つについて、具体的に説明する。この点滴のコンパートメントモデルの１つである２コンパートメントモデル（例示コンパートメントモデルと称する。）では、以下の７つの薬物動態パラメータが定義される。つまり、CL:全身クリアランス(単位はリットル／時間(L/hr))、V₁:体循環コンパートメントの分布容積(L)、V₂:末梢コンパートメントの分布容積(L)、Q:体循環コンパートメントと末梢コンパートメントの間のクリアランス(L/hr)、k₁₀:消失速度定数(1/hr)、k₁₂:体循環から末梢コンパートメントへの薬物移行を表す消失速度定数(1/hr)、k₂₁:末梢から体循環コンパートメントへの薬物移行を表す消失速度定数(1/hr)である。また、例示コンパートメントモデルの薬物動態パラメータでは、以下の関係式が成り立つ。つまり、CL=k₁₀V₁,Q=k₁₂V₁,Q=k₂₁V₂である。

この場合、７つのパラメータに対して関係式が３つであるため、未知のパラメータ数は４つとなる。本発明に係る薬物動態パラメータの推定方法によって推定する薬物動態パラメータの選択については、原則として任意であるが、コンパートメントモデルが有する関係式によって選択可能な薬物動態パラメータが制限される場合がある。例えば、例示コンパートメントモデルの場合、CL=k₁₀V₁,関係式が存在するために、推定する薬物動態パラメータとして、(V₁,k₁₀,k₁₂,k₂₁)の４つを選択することはできるが、(V₁,k₁₀,CL,k₂₁)の４つを選択することはできない。

本発明に係る方法による薬物動態パラメータの推定は、図１に示されるパーソナルコンピュータなどの電子計算機１０で構成されたシステムが当該方法のプログラムを実行することによって実現される。つまり、上記プログラムは、本発明に係る薬物動態パラメータの推定プログラムに相当する。また、本実施形態において、上記プログラムは、後述するようにRAM１２に格納されているが、上記プログラムは、ＣＤ−ＲＯＭなどの記録媒体に格納されていてもよい。つまり、本発明に係る薬物動態パラメータの推定プログラムは、上記プログラムが記録されたＣＤ−ＲＯＭなどの記録媒体として捉えられてもよい。

電子計算機１０は、中央演算装置(CPU)１１とRAM１２と大容量記憶部１３とASIC１４とを有する計算機本体２０、入力部１５、及び表示部１６を備えている。また、計算機本体２０は、CPU１１、RAM１２、大容量記憶部１３、及びASIC１４を相互に接続する内部バス２１を備えている。

CPU１１は、RAM１２に格納されたプログラムに従って演算を実行する。

大容量記憶部１３は、大容量の情報を格納可能な記憶媒体であり、ハードディスク、光磁気ディスク、フラッシュメモリなどで構成されている。大容量記憶部１３には、本発明に係る方法によって薬物動態パラメータの推定を実行可能なプログラムや、血中濃度に関する測定データなどが格納されている。大容量記憶部１３に格納された情報は、電子計算機１０の電源がオフとなっても消去されない。

RAM１２は、大容量記憶部１３よりも小容量の情報を格納可能な記憶媒体である。RAM１２には、上記プログラム、及びCPU１１が上記プログラムを実行する際に用いるデータや信号などが、一時的に格納される。

ASIC１４は、入力部１５及び表示部１６と接続されている。ASIC１４には、入力部１５及び表示部１６とCPU１１との信号のやり取りを仲介する回路が組み込まれている。

入力部１５は、後述する患者に関する個人情報、患者に投与する医薬品情報、及び患者の血中濃度に関する測定データを入力する際に使用されるキーボード１７やマウス１８などで構成されている。

表示部１６は、液晶ディスプレイなどで構成されており、入力部１５によって入力された個人情報、医薬品情報、及び測定データや、後述する推定ステップにおいてCPU１１による演算によって推定された薬物動態パラメータに基づく血中濃度推移線などが表示される（図３〜図６参照）。

［薬物動態パラメータを推定するためのシステムの処理の概要］
CPU１１がRAM１２に格納された上記プログラムを実行することによって、薬物動態パラメータを推定する処理が図２のフローチャートに従って実行される。以下に、当該処理の概要が説明される。

電子計算機１０のユーザがマウス１８を操作することによって、表示部１６に表示された上記プログラムを実行するためのアイコン（不図示）がダブルクリックされると、上記プログラムが大容量記憶部１３からRAM１２にコピーされて、CPU１１がRAM１２にコピーされた上記プログラムを実行する。これにより、表示部１６には、図３に示されるようなデータ入力画面３１が表示される。データ入力画面３１には、患者に関する個人情報が入力される個人情報入力領域８１、患者に投与する医薬品情報が入力される医薬品情報入力領域８２、患者の血中濃度に関する測定データが入力される測定データ入力領域８３、及び「データ入力完了」ボタン３４が設けられている。

ユーザは、データ入力画面３１を参照しつつ入力部１５（キーボード１７及びマウス１８）を操作することによって、上述した各領域８１、８２、８３に、個人情報、医薬品情報、及び測定データを入力する（図２のステップＳ１０）。

図３に示されるように、個人情報入力領域８１に入力される個人情報は、患者の性別、年齢、身長、体重である。また、医薬品情報入力領域８２に入力される医薬品情報は、患者に投与される医薬品の品名、剤形などである。なお、個人情報入力領域８１及び医薬品情報入力領域８２は、上述した情報以外の情報が入力されるように構成されていてもよい。

また、測定データ入力領域８３に入力される測定データは、患者の血中濃度を測定した時刻と、当該時刻における血中濃度とで構成される。例えば、血中濃度が１５分間隔且つ２４時間に亘って測定された場合、測定時刻と当該測定時刻に対応する血中濃度とのデータ数は、それぞれ９６個である。図３では、測定データは、キーボード１７の操作による測定時刻と血中濃度との入力と、マウス１８の操作による「次データの入力」ボタン３２のクリックとが交互に実行されることによって、測定時刻と血中濃度とが１組ずつ入力される。そして、「測定データ入力完了」ボタン３３がクリックされることによって、今まで入力された複数組の測定時刻と血中濃度とが測定データとして確定される。なお、測定データは、上述のように１組ずつ入力されるのではなく、大容量記憶部１３に予めデータベースとして格納されていてもよい。

ユーザは、個人情報、医薬品情報、及び測定データを入力を完了すると、マウス１８を操作することによって、「データ入力完了」ボタン３４をクリックする。これにより、後述する方法選択ステップが実行され（図２のステップＳ２０）、薬物動態パラメータを推定する方法として、最尤法またはベイズ推定が選択される。最尤法が選択された場合、表示部１６に、図４に示されるようなモデル選択画面３５が表示される。一方、ベイズ推定が選択された場合、表示部１６に、図５に示されるようなモデル選択画面３６が表示される。

モデル選択画面３５には、最尤法において選択可能な３つのモデル（正規分布モデル１、正規分布モデル２、及び対数正規分布モデル）が表示された最尤法モデル選択領域８４、「推定実行」ボタン３７、及び「戻る」ボタン３９が設けられている。また、モデル選択画面３６には、ベイズ推定において選択可能な３つのモデル（正規分布モデル１、正規分布モデル２、及び対数正規分布モデル）が表示されたベイズ推定モデル選択領域８５、「推定実行」ボタン３８、及び「戻る」ボタン４０が設けられている。

ユーザは、マウス１８の操作によって、各モデル選択領域８４、８５に表示された３つのモデルの何れかを選択し、「推定実行」ボタン３７、３８をクリックする。これにより、後述する式選択ステップが実行される（図２のステップＳ３０）。なお、モデル選択画面３５、３６の「戻る」ボタン３９、４０がクリックされると、表示部１６に表示される画面がモデル選択画面３５、３６からデータ入力画面３１へと遷移する。

ステップＳ３０に続けて、ステップＳ３０において選択された式に基づいて、後述する推定ステップが実行される（図２のステップＳ４０）。推定ステップが実行されることにより、CPU１１が後述する演算を実行する。これにより、薬物動態パラメータが推定される。

ステップＳ４０が実行されて薬物動態パラメータが推定されると、図６に示される血中濃度推移画面４１が表示される。血中濃度推移画面４１には、描画領域８６及び「戻る」ボタン４４が設けられている。

描画領域８６は、横軸が時刻（日時）であり縦軸が血中濃度(mEq/L)であるグラフを描画するための領域である。描画領域８６には、２本の血中濃度推移線４２、４３が描画される（図２のステップＳ５０）。

図６中に細線で示された血中濃度推移線４２は、現行の薬物の投与スケジュールに基づく血中濃度推移線である。現行の薬物の投与スケジュールに基づく血中濃度推移線は、ステップＳ１０において入力された測定データや、過去の多数の患者に対する測定データなどに基づいて描画されるグラフである。

図６中に太線で示された血中濃度推移線４３は、ステップＳ１０において入力された個人情報、医薬品情報、及び測定データや、ステップＳ４０において推定された薬物動態パラメータに基づいて描画されるグラフであり、現行の薬物の投与スケジュールに基づく血中濃度推移線よりも、医薬品を投与される患者に適した投与スケジュールに基づく血中濃度推移線である。

なお、各種情報、各種データ、及びステップＳ４０において推定された薬物動態パラメータなどから血中濃度推移線４２、４３を導き出す過程は、ここでは、その詳細な説明を省略する。

血中濃度推移画面４１の「戻る」ボタン４４がクリックされると、表示部１６に表示される画面が血中濃度推移画面４１から直近に表示されていたモデル選択画面３５、３６へと遷移する。

［方法選択ステップ］
以下、図２のステップＳ２０において実行される方法選択ステップが説明される。方法選択ステップは、CPU１１によって実行される。CPU１１は、未知の薬物動態パラメータの個数と、測定回数とを比較する。

未知の薬物動態パラメータの個数は、対象となった薬剤によって決まったコンパートメントモデルが有する複数の薬物動態パラメータのうち、当該コンパートメントモデルが有する関係式によって制限されていない範囲で、推定すべき薬物動態パラメータとして選択された薬物動態パラメータの個数である。また、測定回数は、図２のステップＳ１０において入力された測定データの組数である。

CPU１１は、測定回数が未知の薬物動態パラメータの個数よりも多い場合、最尤法を選択し、表示部１６にモデル選択画面３５（図４参照）を表示させる。このような場合に最尤法が選択される理由は、このような場合においては方程式の数を変数の個数以上とすることができるために、連立方程式によって未知の薬物動態パラメータを求めることができるからである。よって、測定回数が未知の薬物動態パラメータの個数よりも多い場合、連立方程式を解くことによって変数の値、つまり薬物動態パラメータを求めることができる。

一方、CPU１１は、測定回数が未知の薬物動態パラメータの個数以下であり、薬物動態パラメータの母集団パラメータがあり、且つ当該母集団パラメータが互いに独立である場合、ベイズ推定を選択し、表示部１６にモデル選択画面３６（図５参照）を表示させる。ここで、薬物動態パラメータの母集団パラメータは、過去における多数の患者の血中濃度の測定から求められた薬物動態パラメータのデータ群であり、大容量記憶部１３に格納されている。

なお、測定回数が未知の薬物動態パラメータの個数以下であるが、薬物動態パラメータの母集団パラメータがない場合、或いは当該母集団パラメータがあるが当該母集団パラメータに含まれる薬物動態パラメータが互いに独立でない場合、CPU１１は、最尤法及びベイズ推定の何れも選択せずにその旨のメッセージを表示部１６に表示させてもよいし、ベイズ推定ではなく最尤法を選択するために更に測定データを入力して測定回数を増やすように要求するメッセージを表示部１６に表示させてもよい。

［式選択ステップ］
以下、図２のステップＳ３０において実行される式選択ステップが説明される。式選択ステップは、ユーザがモデルを選択するステップと、CPU１１がユーザによって選択されたモデルに応じて、後述する推定ステップ（図２のステップＳ４０）において用いる基本式に決定するステップとを含む。ここで、基本式は、薬物動態パラメータx_iのベクトルx=(x₁,x₂,・・・,x_n)^Tの関数s(x)で表される。

ユーザがモデルを選択するステップは、上述したように、ユーザがモデル選択領域８４、８５（図４及び図５参照）に表示された３つのモデルの何れかを選択して「推定実行」ボタン３７、３８をクリックするステップである。

ＣＰＵ１１は、方法選択ステップ（図２のステップＳ２０）で最尤法が選択された場合において、正規分布モデル１が選択された場合、以下の（式１）を基本式とし、正規分布モデル２が選択された場合、以下の（式２）を基本式とし、対数正規分布モデルが選択された場合、以下の（式３）を基本式とする。

ここで、m：血中濃度測定回数、t_j：j回目の測定時刻(j=1,2,・・・,m)
c_j：時刻t_jにおける血中濃度測定値、c(t_j)=c(t_j,x)：時刻t_jにおける血中濃度理論値、σ-_c(但し、"-"は文字σのアッパーラインの意味である。)：jによらず同じ値をとる場合のlog(c_j)の標準偏差、ω_c：σ_j／c(t_j,x)で与えられる変動係数である。

なお、（式１）は、時刻t_jにおける血中濃度測定値c_jが正規分布N(c(t_j),σ_c ²)に従う場合、つまり誤差の分散が一定の場合の最尤法の基本式であり、（式２）は、血中濃度測定値c_jが正規分布N(c(t_j),{ω_c×c(t_j)}²)に従う場合、つまり誤差の分散が血中濃度に比例する場合の最尤法の基本式であり、（式３）は、上記血中濃度測定値c_jの対数値log(c_j)が正規分布N(log{c(t_j)},(σ-_c ²))(但し、"-"は文字σのアッパーラインの意味である。)に従う場合、つまり誤差の分散が一定の場合の最尤法の基本式である。

また、CPU１１は、方法選択ステップ（図２のステップＳ２０）でベイズ推定が選択された場合において、正規分布モデル１が選択された場合、以下の（式４）を基本式とし、正規分布モデル２が選択された場合、以下の（式５）を基本式とし、対数正規分布モデルが選択された場合、以下の（式６）を基本式とする。

ここで、m：血中濃度測定回数、n：薬物動態パラメータ数、t_j：j回目の測定時刻(j=1,2,・・・,m)、c_j：時刻t_jにおける血中濃度測定値、c(t_j)=c(t_j,x)：時刻t_jにおける血中濃度理論値、σ_c：jによらず同じ値をとる場合のc_jの標準偏差、σ-_c(但し、"-"は文字σのアッパーラインの意味である。)：jによらず同じ値をとる場合のlog(c_j)の標準偏差、μ_i：x_i(i=1,2,・・・,n)の平均、μ-_i(但し、"-"は文字μのアッパーラインの意味である。)：log(x_i)の平均、σ_i：x_iの標準偏差、σ-_i(但し、"-"は文字σのアッパーラインの意味である。)：log(x_i)の標準偏差、ω_c：σ_c／c(t_j,x)で与えられる変動係数、ω_i：σ_i／μ_iで与えられる変動係数、ｘ-_i(但し、"-"は文字xのアッパーラインの意味である。)：log(x_i)である。

なお、（式４）は、時刻t_jにおける血中濃度測定値c_jが正規分布N(c(t_j),σ_c ²)に従い且つ薬物動態パラメータx_iが正規分布N(μ_i,σ_i ²)に従う場合、つまり誤差の分散が一定の場合のベイズ推定の基本式であり、（式５）は、血中濃度測定値c_jが正規分布N(c(t_j),{ω_c×c(t_j)}²)に従い且つ薬物動態パラメータx_iが正規分布N(μ_i,(ω_i×μ_i)²)に従う場合、つまり誤差の分散が血中濃度に比例する場合のベイズ推定の基本式であり、（式６）は、血中濃度測定値c_jの対数値log(c_j)が正規分布N(log{c(t_j)},(σ-_c ²))(但し、"-"は文字σのアッパーラインの意味である。)に従い且つ薬物動態パラメータx_iの対数値log(x_i)が正規分布N(μ-_i,σ-_i ²)(但し、"-"は文字μ,σそれぞれのアッパーラインの意味である。)に従う場合のベイズ推定の基本式である。なお、（式６）は、血中濃度測定値c_jの変動が大きい場合に適用される。

［推定ステップ］
以下、図２のステップＳ４０において実行される推定ステップの詳細が、図７のフローチャートに基づいて説明される。推定ステップは、CPU１１が上記プログラムに従って以下に詳述する演算を実行することによって実現される。

推定ステップは、以下に詳述するように、方法選択ステップ（図２のステップＳ２０）において選択された方法（最尤法またはベイズ推定）における基本式（（式１）〜（式６））である薬物動態パラメータx_iのベクトルx=(x₁,x₂,・・・,x_n)^Tの関数s(x)に対して２階偏微分を実行し、当該２階偏微分によって得られた式から同一の薬物動態パラメータx_iによる２階偏導関数の項を削除した近似式に修正マルカート法を適用して関数s(x)を最小化し、最小化した際のベクトルxを求めることによって、薬物動態パラメータx_iを推定する。

上記のような、修正マルカート法による薬物動態パラメータx_iの推定は、図７のフローチャートに従った処理が実行されることによって実現される。つまり、推定ステップは、初期値準備ステップ（Ｓ２１０）、計算ステップ（Ｓ２２０）、変数変更ステップ（Ｓ２３０）、判定ステップ（Ｓ２４０）、及び決定ステップ（Ｓ２７０）を含む。

なお、図７のフローチャートにおいて、計算ステップ（Ｓ２２０）は、式選択ステップ（図２のステップＳ３０）において基本式とされた（式１）〜（式６）に応じて異なる行列による式が導出されるが、ステップＳ２２０以外においては、（式１）〜（式６）の何れが基本式とされても同処理である。よって、以下の説明では、式選択ステップ（図２のステップＳ３０）において（式１）が基本式とされたと仮定して図７のフローチャートの説明を行い、その後で、他の式（（式２）〜（式６））が基本式とされた場合の計算ステップ（Ｓ２２０）について説明する。

初期値準備ステップ（Ｓ２１０）において、CPU１１は、薬物動態パラメータx_iのベクトルxの初期値x⁽⁰⁾=(x₁ ⁽⁰⁾,x₂ ⁽⁰⁾,・・・,x_n ⁽⁰⁾)^Tを準備する。本実施形態では、以下に詳述するように、初期値x⁽⁰⁾の各成分(x₁ ⁽⁰⁾,x₂ ⁽⁰⁾,・・・,x_n ⁽⁰⁾)^T、つまり各薬物動態パラメータx_iは、ボックス・ミューラー法を用いて作成される。なお、初期値x⁽⁰⁾は、以下に詳述する方法以外の公知の方法によって準備されてもよい。

以下に、初期値x⁽⁰⁾を作成する手順について詳述する。最初に、独立な区間[0,1)の一様乱数を複数組用意する。例えば、統計数理研究所によって提供されている多種多様な物理乱数・算術乱数を用意する。次に、２組ずつの一様乱数からボックス・ミューラー法により、正規分布N(0,1)に従う正規乱数を作成する。なお、ボックス・ミューラー法は、u,vを独立な区間[0,1)の一様乱数として、以下の（式７）とする方法である。

このとき、x,yはN(0,1)に従う独立な正規乱数となる。更に、薬物動態パラメータx_i(i=1,2,・・・,n)に対して、適当な平均値μ_xi及び変動係数ω_xiを考え、一様乱数列ｒ_iを用いて、以下の（式８）として初期値を作成する。

また、初期値準備ステップ（Ｓ２１０）において、CPU１１は、最小値が更新されたか否かを表す指標であるパラメータrの初期値を０とする。準備されたベクトルxの初期値x⁽⁰⁾及びパラメータrの初期値はRAM１２に格納される。

初期値準備ステップ（Ｓ２１０）の実行が完了すると、計算ステップ（Ｓ２２０）が実行される。計算ステップ（Ｓ２２０）において、CPU１１は、最初に、式選択ステップ（図２のステップＳ３０）において選択された基本式である薬物動態パラメータx_iのベクトルx=(x₁,x₂,・・・,x_n)^Tの関数s(x)、つまり（式１）に対して２階偏微分を実行する。以下に、（式１）に対して１階偏微分が実行された（式９）と２階偏微分が実行された（式１０）とが示される。

次に、CPU１１は、当該２階偏微分によって得られた式から血中濃度に対する薬物動態パラメータx_iによる２階偏導関数の項を削除することによって、（式１１）に示される近似式を導出する。

次に、CPU１１は、近似式（式１１）を以下の（式１２）に適用する。ここで、（式１２）はニュートン法の式である。

すると、以下の（式１３）が導出される。そして、（式１３）を変形することにより（式１４）が導出される。

次に、CPU１１は、（式１４）から、ベクトルxのｉ（iは０以上の整数）番目の近似値であるx⁽ⁱ⁾=(x₁ ⁽ⁱ⁾,x₂ ⁽ⁱ⁾,・・・,x_n ⁽ⁱ⁾)^T、及びΔx=(Δx₁,Δx₂,・・・,Δx_n）^Tを含む行列による（式１５）を導出する。ここで、（式１５）においてA,Δx,及びb₁は以下の（式１６）、（式１７）、（式１８）の通りである。また、（式１５）において、Iはm次の単位行列である。また、（式１４）、（式１５）において、λは正のパラメータで、2^r-1である。λを大きくすればするほど、Δxは最急降下方向に近づく。

次に、CPU１１は、行列による（式１５）においてΔxに対するｎ元連立一次方程式を解き、x'=x⁽ⁱ⁾+Δxとおく。

変数変更ステップ（Ｓ２３０）において、CPU１１は、計算ステップ（Ｓ２２０）においてｎ元連立一次方程式を解いた結果、s(x')<s(x⁽ⁱ⁾)であることを条件として、パラメータrをmin(0,r-1)とし且つiに１を加えてx⁽ⁱ⁾=x'とする。一方、s(x')>=s(x⁽ⁱ⁾)であることを条件として、パラメータrに１を加える。つまり、変数変更ステップ（Ｓ２３０）では、パラメータr及び偏微分の階数iの値を調整する。

変数変更ステップ（Ｓ２３０）の実行後に、判定ステップ（Ｓ２４０）が実行される。判定ステップ（Ｓ２４０）において、CPU１１は、第１条件及び第２条件を具備するか否かを判断する。ここで、第１条件は、i>=Mまたはr>=Nであることである。このとき、M,Nは予め設定された正の整数定数である。また、第２条件は、|s(x⁽ⁱ⁺¹⁾)-s(x⁽ⁱ⁾)|<ε且つ|Δx|<δ且つr=0であることである。このとき、ε,δは予め設定された小さな実数正数である。

CPU１１は、第１条件を具備することを条件として、s(x)が収束しないと判定する（Ｓ２５０）。このとき、薬物動態パラメータx_iは推定されることなく、推定ステップは終了する。なお、CPU１１は、その旨のメッセージを表示部１６に表示させてもよい。

一方、CPU１１は、第２条件を具備することを条件として、s(x)が収束すると判定する（Ｓ２６０）。このとき、決定ステップ（Ｓ２７０）が実行される。決定ステップ（Ｓ２７０）において、CPU１１は、s(x)が収束すると判定された時点におけるベクトルx⁽ⁱ⁾の各成分を薬物動態パラメータx_iとする。つまり、決定ステップでは、s(x)の最小値を与える薬物動態パラメータx_iが選定される。

また、一方、CPU１１は、第１条件及び第２条件の何れも具備しないことを条件として、再び計算ステップ（Ｓ２２０）を実行する。この計算ステップ（Ｓ２２０）においては、以前に実行された計算ステップ（Ｓ２２０）及び変数変更ステップ（Ｓ２３０）によってパラメータrや偏微分の階数iが変更されている。そのため、次回の判定ステップ（Ｓ２４０）においては、前回と異なる判定結果となり得る。

以上のように、修正マルカート法では、初期値が設定された上で（Ｓ２１０）、判定ステップ（Ｓ２４０）においてs(x)が収束する（Ｓ２６０）或いは収束しない（Ｓ２５０）と判定されるまで、計算ステップ（Ｓ２２０）及び変数変更ステップ（Ｓ２３０）が繰り返し実行される。つまり、推定ステップでは、計算ステップ（Ｓ２２０）及び変数変更ステップ（Ｓ２３０）が反復されることによって、ｎ元連立非線形方程式が解かれる。

以下、式選択ステップ（図２のステップＳ３０）において（式２）が基本式とされた場合の計算ステップ（Ｓ２２０）について説明する。計算ステップ（Ｓ２２０）において、CPU１１は、最初に、式選択ステップ（図２のステップＳ３０）において選択された基本式である（式２）に対して２階偏微分を実行する。次に、CPU１１は、当該２階偏微分によって得られた式から同一の薬物動態パラメータx_iによる２階偏導関数の項を削除することによって、近似式を導出する。以下に、（式２）に対して１階偏微分が実行された（式１９）と２階偏微分が実行され且つ２階偏導関数の項を削除された（式２０）とが示される。

次に、CPU１１は、近似式（式２０）を（式１２）に適用する。これにより、以下の（式２１）が導出される。そして、（式２１）を変形することにより（式２２）が導出される。

次に、CPU１１は、（式２２）から、ベクトルxのi（iは０以上の整数）番目の近似値であるx⁽ⁱ⁾=(x₁ ⁽ⁱ⁾,x₂ ⁽ⁱ⁾,・・・,x_n ⁽ⁱ⁾)^T、及びΔx=(Δx₁,Δx₂,・・・,Δx_n）^Tを含む行列による（式２３）を導出する。ここで、（式２３）においてA,D,Δx,及びb₂は以下の（式２４）、（式２５）、（式２６）、（式２７）の通りである。また、（式２３）において、Iはm次の単位行列とする。

次に、CPU１１は、行列による（式２３）においてΔxに対するｎ元連立一次方程式を解き、x'=x⁽ⁱ⁾+Δxとおく。

以下、式選択ステップ（図２のステップＳ３０）において（式３）が基本式とされた場合の計算ステップ（Ｓ２２０）について説明する。

（式３）が基本式とされた場合は、x,c_j,c(t_j,x)をそれぞれlog(x),log(c_j),log(c(t_j,x))とすることにより、（式１）が基本式とされた場合に対する計算方法をそのまま使用することが可能である。

以下、式選択ステップ（図２のステップＳ３０）において（式４）が基本式とされた場合の計算ステップ（Ｓ２２０）について説明する。計算ステップ（Ｓ２２０）において、CPU１１は、最初に、式選択ステップ（図２のステップＳ３０）において選択された基本式である（式４）に対して２階偏微分を実行する。次に、CPU１１は、当該２階偏微分によって得られた式から同一の薬物動態パラメータx_iによる２階偏導関数の項を削除することによって、近似式を導出する。以下に、（式４）に対して１階偏微分が実行された（式２８）と２階偏微分が実行され且つ２階偏導関数の項を削除された（式２９）とが示される。

次に、CPU１１は、近似式（式２９）を（式１２）に適用する。これにより、以下の（式３０）が導出される。そして、（式３０）を変形することにより（式３１）が導出される。

次に、CPU１１は、（式３１）から、ベクトルxのi（iは０以上の整数）番目の近似値であるx⁽ⁱ⁾=(x₁ ⁽ⁱ⁾,x₂ ⁽ⁱ⁾,・・・,x_n ⁽ⁱ⁾)^T、及びΔx=(Δx₁,Δx₂,・・・,Δx_n）^Tを含む行列による（式３２）を導出する。ここで、（式３２）においてA,B,Δx,b₁,及びb₃は以下の（式３３）、（式３４）、（式３５）、（式３６）、（式３７）の通りである。また、（式３２）において、Iはm次の単位行列とする。

次に、CPU１１は、行列による（式３２）においてΔxに対するｎ元連立一次方程式を解き、x'=x⁽ⁱ⁾+Δxとおく。

以下、式選択ステップ（図２のステップＳ３０）において（式５）が基本式とされた場合の計算ステップ（Ｓ２２０）について説明する。計算ステップ（Ｓ２２０）において、CPU１１は、最初に、式選択ステップ（図２のステップＳ３０）において選択された基本式である（式５）に対して２階偏微分を実行する。次に、CPU１１は、当該２階偏微分によって得られた式から同一の薬物動態パラメータx_iによる２階偏導関数の項を削除することによって、近似式を導出する。以下に、（式４）に対して１階偏微分が実行された（式３８）と２階偏微分が実行され且つ２階偏導関数の項を削除された（式３９）とが示される。

次に、CPU１１は、近似式（式３９）を（式１２）に適用する。これにより、以下の（式４０）が導出される。そして、（式４０）を変形することにより（式４１）が導出される。

次に、CPU１１は、（式４１）から、ベクトルxのi（iは０以上の整数）番目の近似値であるx⁽ⁱ⁾=(x₁ ⁽ⁱ⁾,x₂ ⁽ⁱ⁾,・・・,x_n ⁽ⁱ⁾)^T、及びΔx=(Δx₁,Δx₂,・・・,Δx_n）^Tを含む行列による（式４２）を導出する。ここで、（式３２）においてA,D,B,Δx,b₂,及びb₄は以下の（式４３）、（式４４）、（式４５）、（式４６）、（式４７）、（式４８）の通りである。また、（式４２）において、Iはm次の単位行列とする。

次に、CPU１１は、行列による（式４２）においてΔxに対するｎ元連立一次方程式を解き、x'=x⁽ⁱ⁾+Δxとおく。

以下、式選択ステップ（図２のステップＳ３０）において（式６）が基本式とされた場合の計算ステップ（Ｓ２２０）について説明する。

（式６）が基本式とされた場合は、x,c_j,c(t_j,x),σ_c,x_i,μ_i,σ_iをそれぞれlog(x),log(c_j),log{c(t_j,x)},log(σ_c),log(x_i),log(μ_i),log(σ_i)とすることにより、（式４）が基本式とされた場合に対する計算方法をそのまま使用することが可能である。

［実施形態の作用効果］
本実施形態は、方法選択ステップ（Ｓ２０）を含む。これにより、血中濃度の測定回数が薬物動態パラメータx_iの個数以上の場合、過去のデータから得られた母集団パラメータを用いることなく、最尤法によって薬物動態パラメータx_iを推定することができる。一方、血中濃度の測定回数が薬物動態パラメータx_iの個数未満の場合、母集団パラメータを用いて、ベイズ推定によって薬物動態パラメータx_iを推定することができる。

また、本実施形態では、修正マルカート法が用いられている。このため、ニュートン法よりも収束性を改善することができる。また、修正マルカート法は、ガウス・ニュートン法とは異なり、２乗和以外の項がある関数に適用した場合でも、収束性を改善することができる。

また、本願の出願人は、最尤法及びベイズ推定の基本式を２階偏微分して得られた式から血中濃度に対して２階偏導関数の項を削除した近似式を用いることによって、関数s(x)が発散する可能性を低くすることができることを発見した。

以上より、本実施形態によると、推定する薬物動態パラメータx_iの数にかかわらず、当該薬物動態パラメータx_iを、発散の可能性を低く抑えつつ推定することができる。

また、本実施形態によれば、最尤法が選択された場合において、血中濃度測定値が正規分布に従う場合（（式１）及び（式２））と対数正規分布に従う場合（式３）とに応じて式を使い分けることができる。また、血中濃度測定値が正規分布に従う場合には更に、血中濃度測定値の誤差の分散が一定の場合（式１）と、血中濃度測定値の誤差の分散が血中濃度に比例する場合（式２）とに応じて式を使い分けることができる。

また、本実施形態によれば、ベイズ推定が選択された場合において、血中濃度の測定において混入した誤差が正規分布に従う場合（（式４）及び（式５））と対数正規分布に従う場合（式６）とに応じて式を使い分けることができる。また、上記誤差が正規分布に従う場合には更に、血中濃度測定値の誤差の分散が一定の場合（式４）と、血中濃度測定値の誤差の分散が血中濃度に比例する場合（式５）とに応じて式を使い分けることができる。

また、本実施形態のように、行列による式を用いることによって、計算が容易となる。

また、本実施形態のようにボックス・ミューラー法を用いることにより、作成された初期値x⁽⁰⁾の各成分を正規分布に従う独立な正規乱数とすることができる。

［変形例１］
上述の実施形態では、CPU１１が、未知の薬物動態パラメータの数と測定回数とを比較し、その比較の結果によって、最尤法とベイズ推定とを選択していた。しかし、最尤法とベイズ推定との選択は、ユーザによって行われてもよい。例えば、図８に示されるように、データ入力画面３１に、最尤法における正規分布モデル１、正規分布モデル２、及び対数正規分布モデル、並びに、ベイズ推定における正規分布モデル１、正規分布モデル２、及び対数正規分布モデルを選択するためのモデル選択領域４５を設けてもよい。これにより、ユーザは、データ入力画面３１において、推定方法及び当該推定方法におけるモデルを選択することができる。この場合、測定回数が未知の薬物動態パラメータの個数よりも多い場合であっても、ベイズ推定を選択することができる。

［変形例２］
上述の実施形態では、各推定方法（最尤法及びベイズ推定）におけるモデル選択は、ユーザがモデル選択画面３５において選択することによって行われていた。しかし、上記モデル選択は、CPU１１によって自動的に行われてもよい。例えば、CPU１１は、データ入力画面３１においてユーザに入力された測定データや、大容量記憶部１３に格納された母集団パラメータなどに基づいて、最適と判断されるモデルを選択してもよい。具体的には、CPU１１は、最尤法が選択されている場合において測定データが概ね正規分布に従い且つ誤差が略一定の場合に、正規分布モデル１を選択する。

［変形例３］
上述の実施形態では、初期値準備ステップ（図７のステップＳ２１０）において、CPU１１は、１つの初期値x⁽⁰⁾を準備したが、複数の初期値x⁽⁰⁾を準備してもよい。この際、各初期値x⁽⁰⁾は、上述の実施形態と同様に、ボックス・ミューラー法を用いて作成されてもよい。

変形例３では、推定ステップは、図９のフローチャートに従って実行される。つまり、複数の初期値x⁽⁰⁾が準備された場合、全ての初期値x⁽⁰⁾の各々について図７のステップＳ２２０〜Ｓ２７０が実行される。そして、CPU１１は、各初期値x⁽⁰⁾に基づく演算結果について、s(x)が収束するか否かを判断し、s(x)が収束する場合にはその時点におけるベクトルx⁽ⁱ⁾を求める。全ての初期値x⁽⁰⁾について、s(x)が収束するか否かが判断されると（Ｓ４１０：Ｙｅｓ）、CPU１１は、収束したベクトルx⁽ⁱ⁾の値を比較し、最も小さいベクトルx⁽ⁱ⁾の各成分を薬物動態パラメータｘ_ｉと決定する（Ｓ４２０）。つまり、CPU１１は、複数の初期値x⁽⁰⁾を準備した上で、初期値x⁽⁰⁾毎にs(x)の値を求め、最も小さい値に収束したs(x)の各成分を薬物動態パラメータx_iと決定する。

薬物動態パラメータx_iの推定過程において、関数s(x)が発散してしまうと、薬物動態パラメータx_iを推定することができない。また、関数s(x)が収束した場合でも、当該収束が本来収束して欲しい解以外の局所解に収束した場合には、適切な薬物動態パラメータx_iを推定することができていない。そこで、変形例３のように、初期値x⁽⁰⁾を複数個準備して、複数の関数s(x)の解を求め、当該解が最も小さく収束したときのx⁽ⁱ⁾の各成分を薬物動態パラメータx_iと決定することで、局所解ではない適切な薬物動態パラメータx_iを推定することができる。

また、初期値準備ステップにおいて、ボックス・ミューラー法を用いた場合、上述の実施形態と同様に、作成された初期値x⁽⁰⁾の各成分を正規分布に従う独立な正規乱数とすることができる。その結果、初期値x⁽⁰⁾を複数個準備して、複数の関数s(x)の解を求めた際に、当該解の少なくとも一部が適切に収束する可能性を高くすることができる。

１０・・・電子計算機
１１・・・CPU
３１・・・データ入力画面
３５・・・モデル選択画面
３６・・・モデル選択画面
４１・・・血中濃度推移画面

Claims (12)

1. 患者から採取された血液における薬物の血中濃度から薬物動態パラメータx_i(i=1,2,・・・,n)を推定する方法であって、
   薬物動態パラメータx_iを推定する方法として最尤法またはベイズ推定の何れかが選択される方法選択ステップと、
   上記方法選択ステップにおいて選択された方法における基本式である上記薬物動態パラメータx_iのベクトルx=(x₁,x₂,・・・,x_n)^Tの関数s(x)に対して２階偏微分を実行し、当該２階偏微分によって得られた式から同一の薬物動態パラメータx_iによる２階偏導関数の項を削除した近似式に修正マルカート法を適用して上記関数s(x)を最小化し、最小化した際の上記ベクトルxを求めることによって、上記薬物動態パラメータx_iを推定する推定ステップと、を含む薬物動態パラメータの推定方法。
2. 最尤法における上記関数s(x)は、
   時刻t_jにおける血中濃度測定値c_jが正規分布N(c(t_j),σ_c ²)に従う場合、以下の（式１）であり、
   上記血中濃度測定値c_jが正規分布N(c(t_j),{ω_c×c(t_j)}²)に従う場合、以下の（式２）であり、
   上記血中濃度測定値c_jの対数値log(c_j)が正規分布N(log{c(t_j)},(σ-_c ²))(但し、"-"は文字σのアッパーラインの意味である。)に従う場合、以下の（式３）であり、
   

   

   
   

   

   
   

   

   
   以上において、
   m：血中濃度測定回数
   t_j：j回目の測定時刻(j=1,2,・・・,m)
   c_j：時刻t_jにおける血中濃度測定値
   c(t_j)=c(t_j,x)：時刻t_jにおける血中濃度理論値
   σ-_c(但し、"-"は文字σのアッパーラインの意味である。)：jによらず同じ値をとる場合のlog(c_j)の標準偏差
   ω_c：σ_j／c(t_j,x)で与えられる変動係数
   であり、
   上記方法選択ステップにおいて最尤法が選択された場合に、当該３つの式の何れかが選択される式選択ステップを更に含む請求項１に記載の薬物動態パラメータの推定方法。
3. ベイズ推定における上記関数s(x)は、
   時刻t_jにおける血中濃度測定値c_jが正規分布N(c(t_j),σ_c ²)に従い且つ薬物動態パラメータx_iが正規分布N(μ_j,σ_i ²)に従う場合、以下の（式４）であり、
   上記血中濃度測定値c_jが正規分布N(c(t_j),{ω_c×c(t_j)}²)に従い且つ薬物動態パラメータx_iが正規分布N(μ_j,(ω_i×μ_i)²)に従う場合、以下の（式５）であり、
   上記血中濃度測定値c_jの対数値log(c_j)が正規分布N(log{c(t_j)},(σ-_c ²))(但し、"-"は文字σのアッパーラインの意味である。)に従い且つ薬物動態パラメータx_iの対数値log(x_i)が正規分布N(μ-_i,σ-_i ²)(但し、"-"は文字μ,σそれぞれのアッパーラインの意味である。)に従う場合、以下の（式６）であり、
   

   

   
   

   

   
   

   

   
   以上において、
   m：血中濃度測定回数
   n：薬物動態パラメータ数
   t_j：j回目の測定時刻(j=1,2,・・・,m)
   c_j：時刻t_jにおける血中濃度測定値
   c(t_j)=c(t_j,x)：時刻t_jにおける血中濃度理論値
   σ_c：jによらず同じ値をとる場合のc_jの標準偏差
   σ-_c(但し、"-"は文字σのアッパーラインの意味である。)：jによらず同じ値をとる場合のlog(c_j)の標準偏差
   μ_i：x_i(i=1,2,・・・,n)の平均
   μ-_i(但し、"-"は文字μのアッパーラインの意味である。)：log(x_i)の平均
   σ_i：x_iの標準偏差
   σ-_i(但し、"-"は文字σのアッパーラインの意味である。)：log(x_i)の標準偏差
   ω_c：σ_c／c(t_j,x)で与えられる変動係数
   ω_i：σ_i／μ_iで与えられる変動係数
   ｘ-_i(但し、"-"は文字xのアッパーラインの意味である。)：log(x_i)
   であり、
   上記方法選択ステップにおいてベイズ推定が選択された場合に、当該３つの式の何れかが選択される式選択ステップを更に含む請求項１に記載の薬物動態パラメータの推定方法。
4. 上記推定ステップは、
   上記ベクトルxの初期値x⁽⁰⁾=(x₁ ⁽⁰⁾,x₂ ⁽⁰⁾,・・・,x_n ⁽⁰⁾)^Tを準備し、最小値が更新されたか否かを表す指標であるrの初期値をr=0と準備する初期値準備ステップと、
   上記式選択ステップにおいて選択された式から、上記ベクトルxのi（iは０以上の整数）番目の近似値であるx⁽ⁱ⁾=(x₁ ⁽ⁱ⁾,x₂ ⁽ⁱ⁾,・・・,x_n ⁽ⁱ⁾)^T、及びΔx=(Δx₁,Δx₂,・・・,Δx_n）^Tを含む行列による式を導出し、当該行列による式において上記Δxに対するn元連立一次方程式を解き、x'=x⁽ⁱ⁾+Δxとおく計算ステップと、
   上記計算ステップにおいて上記n元連立一次方程式を解いた結果、s(x')<s(x⁽ⁱ⁾)であることを条件として、rをmin(0,r-1)とし且つiに１を加えてx⁽ⁱ⁾=x'とし、s(x')>=s(x⁽ⁱ⁾)であることを条件として、rに１を加える変数変更ステップと、
   上記変数変更ステップの実行後に、i>=Mまたはr>=N（但し、M,Nは予め設定された正の整数定数）である第１条件を具備することを条件として、s(x)が収束しないと判定し、|s(x⁽ⁱ⁺¹⁾)-s(x⁽ⁱ⁾)|<ε且つ|Δx|<δ且つr=0（但し、ε,δは予め設定された実数正数）である第２条件を具備することを条件として、s(x)が収束すると判定する判定ステップと、
   上記判定ステップにおいて、s(x)が収束すると判定されたことを条件として、そのときのx⁽ⁱ⁾の各成分を薬物動態パラメータx_iとする決定ステップと、を含み、
   上記計算ステップは、上記初期値準備ステップの実行が完了したこと、または、上記判定ステップにおいて上記第１条件及び上記第２条件を具備しないことを条件として実行される請求項２または３に記載の薬物動態パラメータの推定方法。
5. 上記推定ステップは、
   上記初期値準備ステップにおいて、上記ベクトルxの初期値x⁽⁰⁾を複数個準備し、
   上記決定ステップは、複数個の初期値x⁽⁰⁾の各々について上記計算ステップ、上記変数変更ステップ、及び上記判定ステップが実行された結果、s(x)が最も小さく収束したときのx⁽ⁱ⁾の各成分を薬物動態パラメータx_iと決定する請求項４に記載の薬物動態パラメータの推定方法。
6. 上記初期値準備ステップにおいて、上記ベクトルxの初期値x⁽⁰⁾は、ボックス・ミューラー法によって作成される請求項４または５に記載の薬物動態パラメータの推定方法。
7. 患者から採取された血液における薬物の血中濃度から薬物動態パラメータx_i(i=1,2,・・・,n)を推定するためのプログラムであって、
   薬物動態パラメータx_iを推定する方法として最尤法またはベイズ推定の何れかが選択される方法選択ステップと、
   上記方法選択ステップにおいて選択された方法における基本式である上記薬物動態パラメータx_iのベクトルx=(x₁,x₂,・・・,x_n)^Tの関数s(x)に対して２階偏微分を実行し、当該２階偏微分によって得られた式から同一の薬物動態パラメータx_iによる２階偏導関数の項を削除した近似式に修正マルカート法を適用して上記関数s(x)を最小化し、最小化した際の上記ベクトルxを求めることによって、上記薬物動態パラメータx_iを推定する推定ステップと、を実行させる薬物動態パラメータの推定プログラム。
8. 最尤法における上記関数s(x)は、
   時刻t_jにおける血中濃度測定値c_jが正規分布N(c(t_j),σ_c ²)に従う場合、以下の（式１）であり、
   上記血中濃度測定値c_jが正規分布N(c(t_j),{ω_c×c(t_j)}²)に従う場合、以下の（式２）であり、
   上記血中濃度測定値c_jの対数値log(c_j)が正規分布N(log{c(t_j)},(σ-_c ²))(但し、"-"は文字σのアッパーラインの意味である。)に従う場合、以下の（式３）であり、
   

   

   
   

   

   
   

   

   
   以上において、
   m：血中濃度測定回数
   t_j：j回目の測定時刻(j=1,2,・・・,m)
   c_j：時刻t_jにおける血中濃度測定値
   c(t_j)=c(t_j,x)：時刻t_jにおける血中濃度理論値
   σ-_c(但し、"-"は文字σのアッパーラインの意味である。)：jによらず同じ値をとる場合のlog(c_j)の標準偏差
   ω_c：σ_j／c(t_j,x)で与えられる変動係数
   であり、
   上記方法選択ステップにおいて最尤法が選択された場合に、当該３つの式の何れかが選択される式選択ステップを更に実行させる請求項７に記載の薬物動態パラメータの推定プログラム。
9. ベイズ推定における上記関数s(x)は、
   時刻t_jにおける血中濃度測定値c_jが正規分布N(c(t_j),σ_c ²)に従い且つ薬物動態パラメータx_iが正規分布N(μ_j,σ_i ²)に従う場合、以下の（式４）であり、
   上記血中濃度測定値c_jが正規分布N(c(t_j),{ω_c×c(t_j)}²)に従い且つ薬物動態パラメータx_iが正規分布N(μ_j,(ω_i×μ_i)²)に従う場合、以下の（式５）であり、
   上記血中濃度測定値c_jの対数値log(c_j)が正規分布N(log{c(t_j)},(σ-_c ²))(但し、"-"は文字σのアッパーラインの意味である。)に従い且つ薬物動態パラメータx_iの対数値log(x_i)が正規分布N(μ-_i,σ-_i ²)(但し、"-"は文字μ,σそれぞれのアッパーラインの意味である。)に従う場合、以下の（式６）であり、
   

   

   
   

   

   
   

   

   
   以上において、
   m：血中濃度測定回数
   n：薬物動態パラメータ数
   t_j：j回目の測定時刻(j=1,2,・・・,m)
   c_j：時刻t_jにおける血中濃度測定値
   c(t_j)=c(t_j,x)：時刻t_jにおける血中濃度理論値
   σ_c：jによらず同じ値をとる場合のc_jの標準偏差
   σ-_c(但し、"-"は文字σのアッパーラインの意味である。)：jによらず同じ値をとる場合のlog(c_j)の標準偏差
   μ_i：x_i(i=1,2,・・・,n)の平均
   μ-_i(但し、"-"は文字μのアッパーラインの意味である。)：log(x_i)の平均
   σ_i：x_iの標準偏差
   σ-_i(但し、"-"は文字σのアッパーラインの意味である。)：log(x_i)の標準偏差
   ω_c：σ_c／c(t_j,x)で与えられる変動係数
   ω_i：σ_i／μ_iで与えられる変動係数
   ｘ-_i(但し、"-"は文字xのアッパーラインの意味である。)：log(x_i)
   であり、
   上記方法選択ステップにおいてベイズ推定が選択された場合に、当該３つの式の何れかが選択される式選択ステップを更に実行させる請求項７に記載の薬物動態パラメータの推定プログラム。
10. 上記推定ステップは、
    上記ベクトルxの初期値x⁽⁰⁾=(x₁ ⁽⁰⁾,x₂ ⁽⁰⁾,・・・,x_n ⁽⁰⁾)^Tを準備し、最小値が更新されたか否かを表す指標であるrの初期値をr=0と準備する初期値準備ステップと、
    上記式選択ステップにおいて選択された式から、上記ベクトルxのi（iは０以上の整数）番目の近似値であるx⁽ⁱ⁾=(x₁ ⁽ⁱ⁾,x₂ ⁽ⁱ⁾,・・・,x_n ⁽ⁱ⁾)^T、及びΔx=(Δx₁,Δx₂,・・・,Δx_n）^Tを含む行列による式を導出し、当該行列による式において上記Δxに対するn元連立一次方程式を解き、x'=x⁽ⁱ⁾+Δxとおく計算ステップと、
    上記計算ステップにおいて上記n元連立一次方程式を解いた結果、s(x')<s(x⁽ⁱ⁾)であることを条件として、rをmin(0,r-1)とし且つiに１を加えてx⁽ⁱ⁾=x'とし、s(x')>=s(x⁽ⁱ⁾)であることを条件として、rに１を加える変数変更ステップと、
    上記変数変更ステップの実行後に、i>=Mまたはr>=N（但し、M,Nは予め設定された正の整数定数）である第１条件を具備することを条件として、s(x)が収束しないと判定し、|s(x⁽ⁱ⁺¹⁾)-s(x⁽ⁱ⁾)|<ε且つ|Δx|<δ且つr=0（但し、ε,δは予め設定された実数正数）である第２条件を具備することを条件として、s(x)が収束すると判定する判定ステップと、
    上記判定ステップにおいて、s(x)が収束すると判定されたことを条件として、そのときのx⁽ⁱ⁾の各成分を薬物動態パラメータx_iとする決定ステップと、を実行させ、
    上記計算ステップは、上記初期値準備ステップの実行が完了したこと、または、上記判定ステップにおいて上記第１条件及び上記第２条件を具備しないことを条件として実行される請求項８または９に記載の薬物動態パラメータの推定プログラム。
11. 上記推定ステップは、
    上記初期値準備ステップにおいて、上記ベクトルxの初期値x⁽⁰⁾を複数個準備し、
    上記決定ステップは、複数個の初期値x⁽⁰⁾の各々について上記計算ステップ、上記変数変更ステップ、及び上記判定ステップが実行された結果、s(x)が最も小さく収束したときのx⁽ⁱ⁾の各成分を薬物動態パラメータx_iと決定する請求項１０に記載の薬物動態パラメータの推定プログラム。
12. 上記初期値準備ステップにおいて、上記ベクトルxの初期値x⁽⁰⁾は、ボックス・ミューラー法によって作成される請求項１０または１１に記載の薬物動態パラメータの推定プログラム。
    

Images