第14章　薬物動態学

(注1)　時間tにおける血中濃度の変化速度vは、厳密にはΔtを無限小にしたdtあたりの濃度変化dcで表され、これは血中濃度関数C(t)を時間tで微分したものになります。 したがって血中濃度の変化速度は、厳密には次のように表すことができます。

　(k：比例定数)

このような式を微分方程式といいます。 そして微分方程式を満足するような関数を見つけることを「微分方程式を解く」といいます。 微分方程式を解くには微分の逆演算である積分を利用します。 そこで積分を利用して上の微分方程式を解いてみましょう。 まず関数C(t)を左辺に集め、それから両辺をtで積分します。

　　
← 置換積分： より
ln{C(t)} = -kt + K　(K：両辺の積分定数を合計したもの)
C(t) = e^K-kt = exp(K - kt) = exp(K)exp(-kt)

ここで初期条件「t=0の時、C(t)=C(0)=C₀」から積分定数Kを決定することができます。

検算のためにこの関数をtで微分すると、次のように確かに最初の微分方程式を満足しています。

実は自然界にはこの形式の関数で表される現象が大変多く、自然対数の底eはこの形式の微分方程式に関連して発見されました。 参考までに速度定数の単位を求めてみましょう。 濃度の単位をmg/dl、時間の単位をhrとすると、微分方程式からkの単位は次のようになります。

この単位を血中濃度関数に当てはめると次のようになり、濃度の単位がどんなものでも速度定数の単位は常に[時間単位]^-1になります。

一般に微分方程式を解くのは非常に難しいことです。 しかしコンパートメントモデルで扱うような比較的単純な微分方程式は、ラプラス変換(Laplace transform)という技法を利用すれば簡単に解くことができます。 ラプラス変換はフーリエ変換の親戚であり、微分方程式を積分しやすい形に変換してから積分し、また逆変換して微分方程式の解を得る技法です。 ラプラス変換の解釈は非常に複雑ですが、内容を知らなくても利用できるので安心してください。 (→12.5 スペクトル解析)

ラプラス変換は次のようなものです。

ラプラス変換：
ラプラス逆変換：L^-1{f^*(s)} = f(t)
f(t)：原関数、表　　f^*(s)：像関数、裏　　t ≧ 0 　　s：実数または虚数
○ラプラス変換の性質と例
1) L{Af(t) + Bg(t)} = AL{f(t)} + BL{g(t)} = Af^*(s) + Bg^*(s)：線形性
　
2) 　
3) s = 0の時：
4)
5) 、F(0) = 0の時：
6)
7) a ＞ 0の時： 　　
8) d ≧ 0の時：L{f(t-d)} = exp(-ds)f^*(s)
9) f₁(t)とf₂(t)の合成であるたたみ込み(convolution)に関して
　
　
10) 　　 　　 　　
11) デルタ関数δ(t)(Dirac δ-function、単位インパルス)に関して
　t = 0の時：δ(0) = 1　　t ≠ 0の時：δ(t) = 0
　 　　δ(t) * f(t) = f(t) * δ(t) = f(t)　　L{δ(t)} = 1
12) ステップ関数H_a(t)に関して
　t ＜ aの時：H_a(t) = 0　　t ≧ aの時：H_a(t) = 1
　
13) パルス関数U_ab(t)に関して
　t ＜ a または t ＞ bの時：U_ab(t) = 0 　　a ≦ t ≦ bの時：U_ab(t) = 1
　
14)
15) 　　
16) 　　
17) 　　
18)
19) 　　
20)
21)

ラプラス変換を利用して前述の微分方程式を解いてみましょう。 両辺をラプラス変換すると性質1)と4)から次のようになります。

　　sC^*(s) - C(0) = -k C^*(s) 　　

ここで両辺をラプラス逆変換すると性質14)から次のようになります。

C(t) = C(0)exp(-kt) = C₀ exp(-kt)　(C(0) = C₀)

このように、ラプラス変換を利用すると半ば機械的に微分方程式を解くことができます。

(注2)　一般に関数f(x)が次のような性質を持っている時、その関数のグラフが直線的になることから、その関数は線形(linear)または1次であるといわれます。

1. f(cx) = cf(x)
2. f(x+y) = f(x) + f(y)　(x、y：変数　c：定数)

直線回帰分析で用いる最小2乗法では、関数f(x)はパラメータa(定数)とb(回帰係数)に関して線形です。 そのため誤差の平方和関数Q(a,b)はパラメータに関して2次関数つまり放物線的になり、この関数が最小になる時のパラメータ推定値を解析的に——つまり微分・積分を利用して——求めることができます。 このような場合を線形最小2乗法または単に最小2乗法といいます。 (→第5章　相関分析と回帰分析)

関数f(x)がxに関して線形ではなくてもパラメータに関して線形なら、やはり線形最小2乗法と呼ぶので注意してください。 例えば次のような2次関数はxに関して線形ではなく、グラフは放物線になります。 しかしパラメーターに関しては線形であり、線形最小2乗法によってパラメータの推定値を求めることができます。 (→第7章　重回帰分析)

これに対してコンパートメントモデルの関数C(t)は初期濃度C₀に関しては線形ですが、変数tに関しても速度定数kに関しても線形ではありません。 そのためC(t)に最小2乗法を適用した時の誤差平方和関数Q(C₀,k)は放物線的にならず、この関数が最小になる時のパラメータ推定値を解析的に求めることはできません。 このような場合を非線形最小2乗法といいます。

非線形最小2乗法において、関数が最小になる時のパラメーター推定値を求める方法として次のような手法があります。

1. 滑降シンプレックス法(滑降単体法、Downhill simplex method、Nelder-Mead method)
2. ガウス・ニュートン法(Gauss-Newton method)
3. ニュートン・ラプソン法(Newton-Raphson method)

これらの手法を厳密に説明すると複雑すぎて頭が痛くなるので、簡略かつ概念的に説明してみましょう。

(1) 滑降シンプレックス法

滑降シンプレックス法は、いわば”誤差平方和関数Qの上に乗って”関数の最小値Q_minを直接探索する最も単純な手法です。 図14.1.5を参考にしながら、その手順を説明しましょう。

1) Q(k)が最小値Q_minになる時のパラメーター推定値に近いと思われる適当な初期値k₀を何らかの方法——単なる山カンでも良い——で推定し、その時の関数値Q(k₀) = Q₀を求める。

○パラメーターが1個の時の初期値：Q(k₀) = Q₀
○パラメーターがp個の時の初期値：Q(ｋ₀) = Q₀ 　　

2) k₀を元にして初期シンプレックスを作る。

シンプレックス(Simplex、単体)：{Q¹,Q²,Q³,…,Q^p+1} … p次元空間上に(p + 1)個の点を一定の規則で並べたもの
○パラメーターが1個の時
Q¹ = Q₀=Q(k₀) 　　Q² = Q₁ = Q(k₁ = k₀ + s)
※微小距離sについては、例えばs = k₀×0.1などを用いる。
○パラメーターがp個の時
　　として
Q¹ = Q₀ = Q(ｋ₀) = Q(ｋ¹) 　　Q² = Q(ｋ²) 　　Q³ = Q(ｋ³) 　　…　　Q^p+1 = Q(ｋ^p+1)
　　 　　 　　…　　
※微小距離sについては、例えば絶対値が最小のパラメーターk_j0を用いてs = k_j0×0.1などを用いる。
{Q¹,Q²,Q³,…,Q^p+1}の中から次のような点を探索する。
Q^h = Q(ｋ^h)：Qが最大値になる点 　　Q^s = Q(ｋ^s)：Qが2番目に大きな値の点 　　Q^l = Q(ｋ^l)：Qが最小値になる点
：Q^hを除いたp個の点の重心
どちらの場合もここで6)の収束条件をチェックし、収束条件を満足していれば探索を終了する。

3) シンプレックスを利用してQ(k)の谷底方向を探索する。

4) Q(k)上で斜面を登ってしまった時は別の方向を探索する。

5) こうして谷底に向かって探索を続ければ、いつかはシンプレックスの中に谷底が入る時が来る。

6) こうして谷底を狭い範囲に追い詰めていき、次のようになった時に探索を終了する。

滑降シンプレックス法には次のような特徴があります。

• 安定に収束する。
• 収束速度が遅い。
• Q(k)に局所的な極小点があると、初期値によってはその極小点に収束してしまう。

(2) ガウス・ニュートン法

ガウス・ニュートン法は関数C(ｋ|t)をパラメーターｋに関して線形近似し、線形最小2乗法を用いて近似解を求める手法です。 図14.1.6を参考にしながら、その手順を説明しましょう。

0) 関数C(ｋ|t)をパラメーターｋに関して線形近似するにはC(ｋ|t)をｋ_lでテーラー展開(Tayler expansion)し、高次項を省いて1次式で近似します。

ｙ ≒ ｃ_l + Ｊ_l[ｋ - ｋ_l] + ε = ｃ_l + Ｊ_lｋ - Ｊ_lｋ_l + ε
　　 　　 　　
ヤコビニアン行列(Jacobian matrix)：

この式はｋに関して1次式ですから、この式の誤差平方和関数Q(ｋ)は2次式つまり放物線になります。 そこでQ(ｋ)を最小にする時のｋを線形最小2乗法によって求めると次のようになります。

Q(ｋ) = ε'ε = [ｙ - ｃ_l + Ｊ_lｋ_l - Ｊ_lｋ]'[ｙ - ｃ_l + Ｊ_lｋ_l - Ｊ_lｋ] → 最小

∴ｋ = [Ｊ_l'Ｊ_l]^-1Ｊ_l'[ｙ - ｃ_l + Ｊ_lｋ_l] = [Ｊ_l'Ｊ_l]^-1Ｊ_l'[ｙ - ｃ_l] + ｋ_l

ここでｄ_l = ｋ - ｋ_lと置くと、

と更新すれば最小値ｋ_minに近づくことになります。 ただしこの計算法はデータによって発散しやすい時があるので、次のような方法が提唱されています。

1) 適当な初期値ｋ₀——滑降シンプレックス法の解を用いることが多い——を推定し、その時のC(ｋ₀|t)の値とその傾きつまりＪ₀を求める。 Ｊ₀は適当な微小単位Δｋを用いた差分法によって求める。

ｃ₀のi番目の要素：c_i0 = C(ｋ₀|t_i)
Ｊ₀の[i,j]要素：
C(ｋ₀|t)を線形近似した時のC(ｋ₀ + Δｋ|t)の近似値：ｙ₀ = ｃ₀ + Ｊ₀[Δｋ]

2) Ｊ₀を用いてC(ｋ₀|t)を直線によって近似し、Q(ｋ)を最小にする時の近似解k₁を求める。

3) こうして求めたｋ₁はｋ_minに近づくはずだから、次はｋ₁を用いて1)と2)の手順を繰り返す。

4) 以後も同様の計算を続けていき、次のような収束条件を満足した時にQ_i+1をQminとし、ｋ_i+1をｋ_minとする。

ガウス・ニュートン法には次のような特徴があります。

• 精度良く高速に収束する。
• C(ｋ|t)が直線と大きく異なっていたり、データに誤差が多かったりすると収束せずに発散してしまう。
• 初期値が解と離れすぎていると収束せずに発散してしまう。

(3) ニュートン・ラプソン法

ニュートン・ラプソン法はQ(k)を放物線によって近似し、解析的な近似解を求める手法です。 この手法については第10章で説明したので、計算原理についてはそちらを見てください。 (→10.3 ロジスティック回帰分析の計算方法 (注2))

ニュートン・ラプソン法には次のような特徴があります。

• 精度良く高速に収束する。
• Q(k)が放物線と大きく異なっていたり、データに誤差が多かったりすると収束せずに発散してしまう。
• 初期値が解と離れすぎていると収束せずに発散してしまう。
• Q(k)に局所的な極小点があると初期値によってはその極小点に収束してしまう。

どの手法についてもその欠点を補う変法が考案されていますが、やはり一長一短があります。 そこで普通はこれらの手法を併用して安定かつ精度の良い推定値を得るようにします。 例えば滑降シンプレックス法で一度収束させた値を初期値とし、ガウス・ニュートン法またはニュートン・ラプソン法で再び収束させるといった方法を用いるわけです。

非線形最小2乗法は良い初期値を与えることが重要なポイントです。 コンパートメントモデルの場合、対数変換による回帰直線を利用した近似的な簡便法によって得られたパラメーターを初期値として用いるのが普通です。

(注3)　血中濃度を対数変換してから回帰直線を当てはめる簡便法は次のような最小2乗基準を用います。 この方法で求めた関数は図14.1.3(b)のように血中濃度を片対数でプロットしたデータに最も良く適合します。

誤差平方和関数： → 最小

これに対して非線形最小2乗法による回帰は次のような最小2乗基準を用います。 この方法で求めた関数は図14.1.3(a)のように血中濃度を実測値のままプロットしたデータに最も良く適合します。

誤差平方和関数： → 最小

血中濃度を片対数でプロットするのは、あくまでも計算をしやすくするための技法にすぎず、実際の血中濃度の変化を正しく反映しているわけではありません。 例えば内服モデルの開始時のように、実際の血中濃度が0になるポイントは片対数ではプロットできません。 血中濃度は実測値でプロットするのが本来ですから、対数変換を利用した簡便法よりも非線形最小2乗法によって求めた関数の方が正確です。

また最小2乗法の一種にデータに適当な重みを付けて回帰を行う重み付け最小2乗法という手法があります。 通常、この方法では各データの信頼度つまり分散σ_i²の逆数を重みw_iにし、次のような最小2乗基準を用います。 (→第7章　重回帰分析)

誤差平方和関数： → 最小

この方法で求めた関数のグラフは信頼度の高いデータに近付き、信頼度の低いデータからは離れます。 そのため見かけ上はデータにあまり適合しなくなります。

同一条件で何度も実験を繰り返して、その平均値をデータにした時は平均値の分散つまり標準誤差の平方を求めることができます。 しかし1人の被験者を用いて実験を1回だけ行った時はデータの分散を求めることができません。 そこで分散の平方根つまり標準偏差がデータに比例すると仮定して、データの平方の逆数を重みにした次のような最小2乗基準を用いることがあります。

誤差平方和関数： → 最小

この方法で求めた関数のグラフは値が小さい部分ほどデータに近付き、値が大きい部分はデータから離れます。 そのため対数変換を利用した簡便法で求めた関数と似たものになります。 しかし血中濃度のデータと標準偏差が比例するとは限らないので、何らかの方法で標準偏差が求められない限り重み付けはしない方が無難でしょう。

(注4)　AUCは医学・薬学分野独特の用語であり、数学的には単なる積分値のことです。 また半減期は、本来は物理学において放射性元素が崩壊して放射能が半分になる時間を表す用語です。 静注1コンパートメントモデルの血中濃度関数が放射性元素の崩壊関数と同じ形式になることから、この用語が転用されたのです。 ちなみに物理学では速度定数のことを崩壊定数(decay constant)とか時定数(time constant)などと呼んでいます。

各種の指標は血中濃度関数から次のようにして求めることができます。

　　

∴　(一般に1/m減期：)
　　
∴

これ以外にもモーメント解析(moment analysis)で用いられる指標もあります。 モーメント解析は薬物の体内動態を時間的な広がりを持った確率過程と考え、血中濃度関数を確率分布のように捉えて平均値や分散に相当する指標を求めます。

○AUC：0次モーメント、確率分布における積分値つまり1に相当する

：確率分布における確率密度関数に相当する
○AUMC(Area Under the Moment Curve)：モーメント曲線下面積

○MRT(Mean Residence Time、平均滞留時間)：1次モーメント、確率分布における期待値つまり平均値に相当する

○VRT(Variance of Residence Time、平均滞留時間の分散)：2次モーメント、確率分布における分散に相当する

ただし

静注1コンパートメントモデルの場合、MRTとVRTは次のように排泄速度定数だけで決まります。 つまり結局のところ、モーメント解析でも排泄速度定数が最も重要な値なのです。