玄関雑学の部屋雑学コーナー統計学入門

第14章 薬物動態学

この章では薬物動態学の基礎的な概念と、薬物動態学において中心的な役割を担うコンパートメントモデルについて解説します。

14.1 コンパートメントモデル

(1) 薬物動態学と薬力学

薬物が生体に吸収され、生体内で分布し、酵素等の働きで代謝され、腎臓等から体外に排泄される過程をADME(Absorption・Distribution・Metabolism・Elimination/Excretion、アドメ)といい、薬物の投与(Liberation)まで含めた過程をLADME(ラドメ)といいます。 このADMEについて研究する時、主としてその速度に注目して速度論的なアプローチをする研究法を薬物動態学(pharmacokinetics)といいます。 つまり薬物動態学とは薬物が生体内でどのような影響を受け、どのような挙動をするのかということを研究する学問です。

それに対して薬物が生体に与える影響を、その物質的な活性、受容体や作用部位との相互干渉、用量−反応関係などに注目してアプローチする研究法を薬力学(pharmacodynamics)といいます。 これらは薬理学に現代科学の方法論を取り入れたものであり、数学的なアプローチを多用するので数学アレルギーの人にとっては頭が痛くなるような研究方法です。

(2) コンパートメントモデルの原理

薬物を生体に投与した時、生体全体を薬物処理系として捉え、生体内をいくつかの処理区画に分けて考えると便利です。 例えば薬物を内服投与した時は胃や腸を吸収区画、血液全体を分布区画、肝臓を代謝区画、腎臓を排泄区画と考え、薬物が尿中に排泄されることによって系の外に出ていくと考えるわけです。 これをコンパートメントモデル(compartment model、区画モデル)といいます。

図14.1.1 コンパートメントモデル

コンパートメントモデルでは、薬物の移行速度について次のような物理的仮定をします。

  1. 区画から薬物が排出される速度は、その区画の薬物濃度に比例する。 これを一次速度過程と呼ぶ。
  2. 薬物の移行速度が時間によって変化せず、一定の時は0次過程と呼ぶ。

1番目の仮定を簡単に言えば「濃いものほど速く薄くなる」ということです。 単位時間あたりに区画から排出される血液量が常に一定とすれば、薬物の排出速度はその血中濃度に比例することになります。 つまり薬物はその血中濃度が濃いほど速く排出され、薄くなるにつれてゆっくり排出されるわけです。 2番目の仮定は点滴静注のように一定速度で薬物を投与する場合に相当します。 どちらの仮定も、ある程度妥当な仮定でしょう。

(3) 静注1コンパートメントモデル

最初に一番単純なモデルとして血液区画しかない場合を考えてみましょう。 薬物を血液区画に一気に投与し、瞬間的に一様に分布して、徐々に区画外に排出されるとします。 これは薬物を急速静注(bolus)した場合に相当するので静注1コンパートメントモデルまたは開いた1コンパートメントモデルと呼ばれます。 「開いた」とは薬物が区画に自由に出入りできる状態のことです。

図14.1.2 静注1コンパートメントモデルの模式図

投与直後における薬物の濃度を初期濃度といい、普通はC0と書きます。 薬物の投与量をq0、分布区画の容積をVd(distribution volume)とすると、これらとC0との間には次のような関係があります。

C 0 = q 0 V d または V d = V c = q 0 C 0

この場合は血液区画が分布区画なのでVdはVcとも書かれ、血液区画の容積になります。 ほとんどの場合、現実に測定できるのは血中濃度であり、分布容積を直接測定することはできません。 そこで普通は上の右の式のように初期濃度と投与量から分布容積を逆算します。 ところが薬物は血液中のタンパク質と結合したり血管壁に吸着されたりするので実際の血中濃度は多少薄くなり、初期濃度もそれに応じて薄くなります。 そのため上式から求めた分布容積は初期濃度が薄くなった分だけ多めに計算されることになります。 これを見かけの分布容積といいます。

見かけの分布容積は薬物の性質によってかなり差があり、本当の区画容積の1.5〜2倍程度になります。 人間の場合、血液量は体重の1/12〜1/13なので体重60kgの人なら5リットル程度になり、そのうちの2リットル程度が血漿量になります。 したがって見かけの分布容積は薬物が血液全体に分布するなら8〜10リットル程度、血漿だけに分布するなら3〜4リットル程度になります。

血中濃度が時間とともに変化する様子を表すには、血中濃度Cを時間tの関数C(t)と考え、この関数をグラフ化します。 そのためには、まず関数C(t)の姿を決定する必要があります。 血液区画から薬物が排出される速度vは、単位時間Δtあたりに血中濃度が低下した量ΔCとして表すことができます。 そしてコンパートメントモデルの1番目の仮定から、この排出速度vは血中濃度Cに比例するので次のような関係があることになります。

v = ⊿C ⊿t = -kC (k:比例定数)

この関係式と初期濃度C0から時間−血中濃度関数C(t)の姿を数学的に導くことができます。 それは次のような指数関数であり、図14.1.3のようなグラフになります。 (注1)

静注1コンパートメントモデル関数:C(t) = C0 e-kt = C0 exp(-kt)
図14.1.3 静注1コンパートメントモデル関数(実数単位) 図14.1.4 静注1コンパートメントモデル関数(対数単位)

比例定数kは速度定数(rate constant)と呼ばれていて、[時間単位]-1(通常はhr-1)という単位を持っています。 この値は単位時間あたりに排出される薬物量の割合を表し、ある時間の血中濃度にこの値をかけると、その時の排出速度になります。 このモデルでは排出時の速度定数だけですが、一般には吸収時や分布時の速度定数もあります。 そこでそれらの速度定数と区別するために、普通は「ke」または「kel」と書いて排出速度定数と呼びます。 排出速度定数が大きいほど薬物が速やかに体外に排出されます。

排出速度定数は実際に薬物を投与し、時間を追って血中濃度を測定したデータに基いて計算します。 図14.1.3の黒丸がそのデータをプロットしたものです。 これらのデータから排出速度定数を求めるには非線形最小2乗法という非常に面倒な手法を利用します。 この手法は原理的には回帰直線を求める時の最小2乗法と同じであり、関数がデータに最もフィットするようにパラメータの値を決定します。 しかし関数が直線ではなく曲線のため「非線形」と呼ばれ、そのせいで計算が非常に複雑になるのです。 (注2)

非線形最小2乗法が開発されたのはかなり昔ですが、何しろ計算が恐ろしく複雑なので一部の専門家にしか知られていませんでした。 そこでひと昔前までは近似的な簡便法として対数変換を利用した直線回帰分析が用いられていました。 また現在でも場合によってはこの簡便法が用いられることがあります。 血中濃度関数を自然対数変換すると次のようになります。 そのため血中濃度関数を片対数で描くと図14.1.4のように直線になり、血中濃度データをプロットするとほぼ直線状に並びます。

ln{C(t)} = ln(C0 e-kt) = ln(C0) - kt

そこで測定された血中濃度データを自然対数変換し、そのデータに回帰直線を当てはめます。 すると次のように回帰直線の傾きが速度定数kの推測値になり、定数を指数変換した値が初期濃度C0の推測値になります。

ln{C(t)} = a - bt = ln(C0) - kt   C0 = exp(a)  k = b

こうして求めた関数が図14.1.4の青色の破線で描いた直線です。 この図では、この青色の破線の方が非線形最小2乗法で求めた黒色の実線よりもデータにフィットしているように見えます。 しかしこのグラフは片対数表示ですから、血中濃度の大きい部分が縮小され、小さい部分が拡大されています。 そのため血中濃度の大きい部分における青色の破線とデータとのズレが実測値ではかなり大きく、実際には黒色の実線の方がデータによりフィットしています。

このように、この簡便法で求めた関数はデータに最もフィットするものではないものの、近似関数としてなら利用することができます。 この方法は計算が簡単なので以前はよく用いられていました。 しかしコンピュータを手軽に利用できる現在では、この方法で求めたパラメーターの値を初期値とした非線形最小2乗法によってデータに最もフィットする関数を求めることができます。 (注3)

(4) 薬物動態学的指標

血中濃度関数C(t)の具体的な姿が決まると、薬物の作用を反映すると思われる次のような薬物動態学的指標を求めることができます。 (注4)

○曲線下面積(AUC:Area Under the concentration Curve)
 無限時間後のAUC:
 t時間後のAUC:
○半減期(half time):
○最大血中濃度:Cmax = C0
○最大血中濃度時間:tmax = 0
○有効血中濃度持続時間:
 CED:閾血中濃度(有効血中濃度の下限)

最も直接的に薬効を反映すると考えられるのは有効血中濃度持続時間tEDです。 上式から、この指標は初期血中濃度C0の対数変換値に正比例し、排出速度定数kに反比例することがわかります。 そして初期血中濃度は薬物の投与量q0と比例するので、薬効は投与量の対数変換値つまり対数用量と比例することになります。 第13章で説明したように薬物の対数用量と反応がほぼ正比例するのは、これが理由のひとつと考えられています。 (→第13章 用量反応解析)

有効血中濃度持続時間を求めるには有効血中濃度の下限つまり閾血中濃度CEDを決定する必要があります。 しかしこの値が判明している薬物は残念ながらほとんどありません。 そこで代わりに曲線下面積AUCを薬効の指標にすることが多くなります。 この値は薬物の直接的な作用を反映するわけではないものの、あえて解釈すれば薬物が生体に与える全体的な作用を反映すると考えられます。

半減期t1/2は薬物の血中濃度または体内の薬物量が半分に減るまでの時間を表します。 しかしこの値が文字通りの意味を持つのは、実は静注1コンパートメントモデルだけです。 そのため薬効の目安というよりも、むしろ排出速度定数kを感覚的に理解しやすい時間単位の値に変換した指標と解釈する方が良いと思います。 つまり半減期が長いほど排出速度定数が小さくなり、その結果、区画からの排出速度が遅くなると解釈することができるわけです。

最大血中濃度Cmax最大血中濃度時間tmaxは薬物の急性毒性的な作用を反映すると考えられています。 静注モデルの場合、最大血中濃度は初期濃度C0になり、最大血中濃度時間は0つまり投与直後になります。

上式からわかるように、これらの指標は全て初期濃度と排出速度定数に依存します。 しかし初期濃度は薬物の投与量によって決まるので一番重要な値は排出速度定数ということになります。

以上のように、コンパートメントモデルでは血中濃度関数と速度定数が主役として活躍し、薬理作用を反映する指標は全てこれらに基いて数学的に導くことができます。 コンパートメントモデルを用いず、例えば実際の血中濃度データを折れ線で結び、その折れ線グラフに基づいてAUCを近似的に計算することがあります。 しかしそれはあくまでも近似的な代用品にすぎません。 例えば検量線を引く時、回帰直線を用いずに、実際のデータとデータを折れ線で結んで検量線として用いるようなものであり、正確さを欠き、しかも普遍性がありません

コンパートメントモデルを用いると、全てのデータを効率的に使用して誤差の少ない指標を求めることができます。 また血中濃度関数を利用すれば実際にデータが測定されていない時間についても血中濃度を推測することができる上、反復投与した時の血中濃度の変化をシミュレートすることもできます。

そしてコンパートメントモデルを用いる最も重要な点は、

「観測データを総合して普遍的な原理を帰納的に推理・洞察し、それに基いて仮定を立て、その仮説が諸々の事実を説明できるかどうかを演繹的に検証する」

という仮説演繹法と呼ばれる現代科学の基本的な手法を、物理学などと同じように数学という科学言語を使って厳密に行う点にあります。 単なるカンや経験——実はこれらも科学には大切な要素なのですが——だけで物を言ったり、実験結果の現象論的な羅列だけで終始しやすい薬学を現代的な学問に変革する意味で、これは大きな進歩と言えるでしょう。 (→1.8 科学的研究の種類)


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

 (k:比例定数)

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

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

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

C(0) = C0 = exp(K)exp(-k×0) = exp(K)
∴C(t) = C0 exp(-kt)

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

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

[mg/dl]/hr = k[mg/dl]
∴k = 1/hr = hr-1

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

[mg/dl] = [mg/dl] × exp(hr-1×hr) = [mg/dl]

一般に微分方程式を解くのは非常に難しいことです。 しかしコンパートメントモデルで扱うような比較的単純な微分方程式は、ラプラス変換(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) f1(t)とf2(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) ステップ関数Ha(t)に関して
 t < aの時:Ha(t) = 0  t ≧ aの時:Ha(t) = 1
 
13) パルス関数Uab(t)に関して
 t < a または t > bの時:Uab(t) = 0   a ≦ t ≦ bの時:Uab(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) = C0 exp(-kt) (C(0) = C0)

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

(注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章 重回帰分析)

y = f(x) = α + βx + γx2 (x:変数 α、β、γ:パラメータ)

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

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

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

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

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

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

図14.1.5 滑降シンプレックス法の概念図

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

○パラメーターが1個の時の初期値:Q(k0) = Q0
○パラメーターがp個の時の初期値:Q(0) = Q0   

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

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

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

○パラメーターが1個の時
Q2 < Q1ならQ(k)上で谷底に向かっているので、k0をk1で置き換えて2)に戻る。
これはQ1 = Q1 = Q(k1)   Q2 = Q2 = Q(k2=k1 + s)と置き換えて、同じ方向に探索を続けることに相当する。
○パラメーターがp個の時
3)-1 谷底方向を探索するためにQh反射または折り返す(refrection)
 r = (1+α)m - αh = m + α(m - h) (α > 0)
 例えばα = 1とするとr = 2m - h
 Qr = Q(r):シンプレックスの最大点Qhを重心Qmを中心にして折り返した点
3)-2 Qr < Qlなら谷底方向に折り返したので、その方向にシンプレックスを拡大(expansion)する。
 e = γr + (1-γ)m = m + γ(r - m) (γ > 1)
 例えばγ = 2とするとe = 2r - m
 Qe = Q(e):Qrをその方向に拡大した点。
3)-2-1 Qe < Qrなら拡大が成功したので、heで置き換えてシンプレックスを更新し、2)に戻る。
3)-2-2 Qe ≧ Qrなら行き過ぎなので、hrで置き換えてシンプレックスを更新し、2)に戻る。

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

○パラメーターが1個の時
Q2 > Q1ならQ(k)上で斜面を登ったので、sの符号を反対にして2)に戻る。
これはとしてQ2 = Q2 = Q(k2 = k1 - 2s = k0 - s)と置き換えて、反対方向に探索を続けることに相当する。
○パラメーターがp個の時
Qs ≧ Qr ≧ Qlなら谷底に対して斜め方向に折り返したので、hrで置き換えてシンプレックスを更新し、2)に戻る。

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

○パラメーターが1個の時
Qi+1 > Qi < Qi-1つまり前回はQ2<Q1で今回はQ2>Q1なら谷底をまたいだので、歩幅をs=0.5×sと半分にして2)に戻る。
これはシンプレックス{Q1,Q2}を縮小して探索を続けることに相当する。
○パラメーターがp個の時
5)-1 Qr ≧ Qhなら折り返す前のシンプレックスの中に谷底が入っているので、シンプレックスを縮小(contraction)する。
 c = βh + (1-β)m = m + β(h - m) (0 < β < 1)
 例えばβ = 0.5とするとc = 0.5(m + h)
 Qc=Q(c):QhをQm方向に縮小した点
5)-1-1 Qc < Qhなら縮小が成功したので、hcで置き換えてシンプレックスを更新し、2)に戻る。
5)-1-2 Qc ≧ Qhなら縮小が失敗したので、l以外の点をl方向に近づけた小さなシンプレックスに更新し、2)に戻る。
 j = 0.5(j + l) (j=1,…,p+1)
5)-2 Qh > Qr > Qsなら折り返した後のシンプレックスの中に谷底が入っているので、hrで置き換えて、あらためてhにしてからシンプレックスを縮小する。
 5)-1と同じ手順でQc = Q(c)を求める。
5)-2-1 Qc < Qhなら縮小が成功したので、hcで置き換えてシンプレックスを更新し、2)に戻る。
5)-2-2 Qc ≧ Qhなら縮小が失敗したので、l以外の点をl方向に近づけた小さなシンプレックスに更新し、2)に戻る。
 j = 0.5(j + l) (j=1,…,p+1)

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

○・パラメーターが1個の時
収束条件:|Qi+1-Qi| < |Qi+1|×ε=|Qi+1|×10-3〜-6 → 収束後のQi+1をQminとし、ki+1をkminとする。
○パラメーターがp個の時
収束条件:|Qh-Ql| < |Ql|×ε=|Ql|×10-3〜-6 → 収束後のQlをQminとし、klをkminとする。
どちらの場合もQminの相対誤差は10-3〜-6未満になる。

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

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

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

図14.1.6 ガウス・ニュートン法の概念図

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


l + l[ - l] + ε = l + l - ll + ε
        
ヤコビニアン行列(Jacobian matrix):

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

Q() = ε'ε = [ - l + ll - l]'[ - l + ll - l] → 最小

= [l'l]-1l'[ - l + ll] = [l'l]-1l'[ - l] + l

ここでl = - lと置くと、

l = [l'l]-1l'[-l]
l+1 = l + l

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

・レーベンバーグ・マーカート法(Levenberg-Marquardt method)
 μ>0の適当な値をl'lの対角要素に上乗せして発散しにくくする。
 l = [l'l + μ]-1l'[-l]
・レーベンバーグ・マーカート・杉本変法
 μの初期値を大きな値にし、収束するにしたがって小さくしていき、安定したところでμ=0として本来のガウス・ニュートン法につなぐ。
 l = [l'l + μl]-1l'[-l]
・重み付けガウス・ニュートン法
 適当な重み行列を用いて発散しにくくする
 l = [l'l]-1l'[-l]

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

0のi番目の要素:ci0 = C(0|ti)
0の[i,j]要素:
C(0|t)を線形近似した時のC(0|t)の近似値:0 = 0 + 0]

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

0 = [0'0]-10'[0 - 0]
1 = 0 + 0

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

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

収束条件:|Qi+1-Qi| < |Qi+1|×ε = 10-3〜-6
Qminの相対誤差は10-3〜-6未満になる。

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

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

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

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

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

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

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

誤差平方和関数: → 最小

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

誤差平方和関数: → 最小

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

また最小2乗法の一種にデータに適当な重みを付けて回帰を行う重み付け最小2乗法という手法があります。 通常、この方法では各データの信頼度つまり分散σi2の逆数を重みwiにし、次のような最小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は次のように排泄速度定数だけで決まります。 つまり結局のところ、モーメント解析でも排泄速度定数が最も重要な値なのです。