基礎

コンパートメントモデル解析

核医学検査におけるコンパートメントモデル解析について
コンパートメントモデル解析を論じた書籍や論文等は、モデルを解析する式が導出される過程が詳細に記述されていないものが多く、理解が困難である。
コンパートメントモデルの式を解く方法は、２種類ある。
1. ラプラス変換の公式を用いてコンパートメントモデルの微分方程式を解く。
2. 積分法を用いてコンパートメントモデルの微分方程式を漸化式に変換して近似解を算出する。
ラプラス変換とは、電気回路の解析などに用いられる演算法である。
電気回路学の書籍をみると、電気回路の解析式は、コンパートメントモデル
解析の式と同じ型をしているものが多くある。
関数 f(t) のラプラス変換 L( f(t) ) を F(s) とする。ここで t は時間、s は複素数
で、複素平面上の点を表す。
L( f(t) ) = F(s) = ∫ f(t) exp(-st) dt（ ラプラス変換の公式 ）ラプラス変換はフーリエ変換に類似している。畳込みの定理も成り立つ。
微分や指数関数を四則演算に変換できる。
畳込みの定理指数関数

L(f) L(g) = L( f＊g )L( exp(at) ) = 1 / ( s – a)

微分
定数関数

L( f’ ) = s L(f) – f(0)L( C ) = C / s

ラプラス変換を用いて、シングルコンパートメントモデルの式を導出してみる。

抵抗とコンデンサの結合回路の式と同じである。

基礎

ラプラス変換の公式

基礎

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

ラプラス変換を用いて、２コンパートメントモデルの式を導出する。

基礎

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

ラプラス変換を用いて、FDGの３コンパートメントモデルの式を導出する。

Xe(0) = 0
Xm(0) = 0
Xe’ = K1 Cp - (k2+k3) Xe + k4 Xm

Xm’= k3 Xe – k4 Xm

s L(Xe) = K1 L(Cp) – (k2+k3) L(Xe) + k4 L(Xm)s L(Xm) = k3 L(Xe) – k4 L(Xm)

--- ①
--- ②

( L(f’) = sL(f) – f(0) )②→
L(Xm) = (k3/(s+k4)) L(Xe)--- ③①③ → (s+k2+k3-k3k4/(s+k4)) L(Xe) = K1 L(Cp)
L(Xe) = ( (s+k4) / ( s2+(k2+k3+k4)s+k3k4) ) K1 L(Cp)= K1 ( (s+k4)/((s+α1)(s+α2)) ) L(Cp)--- ④④ → (s+k4)/((s+α1)(s+α2)) = P/(s+α1) + Q/(s+α2) --- ⑤⑤x(s+α1) → (s+k4)/(s+α2) = P + Q (s+α1)/(s+α2)
s = -α1 → (-α1+k4)/(-α1+α2) = P--- ⑥⑤x(s+α2) → (s+k4)/(s+α1) = P (s+α2)/(s+α1) + Q
s = -α2 → (-α2+k4)/(-α2+α1) = Q--- ⑦④⑤⑥⑦ → L(Xe) = K1 ( P/(s+α1) + Q/(s+α2) ) L(Cp)
L(Xe) = (K1/(α2-α1)) ( (k4-α1)(1/(s+α1)+(α2-k4)(1/(s+α2) ) L(Cp)= (K1/(α2-α1)) ( (k4-α1)L(exp(-α1t)+(α2-k4)L(exp(-α2t) ) L(Cp)
Xe = (K1/(α2-α1)) ( (k4-α1)exp(-α1t)+(α2-k4)exp(-α2t) )＊Cp(＊ convolution )② → L(Xm) = (k3/(s+k4)) L(Xe)= (k3/(s+k4)) K1 ( (s+k4)/ ((s+α1)(s+α2)) ) L(Cp)= K1k3 ( 1/ ((s+α1)(s+α2)) ) L(Cp)--- ⑧
1/ ((s+α1)(s+α2)) = P/ (s+α1) + Q/(s+α2)⑨x(s+α2) → 1/(s+α1) = P (s+α2)/(s+α1) + Q
s = -α2 → 1/(-α2+α1) = Q⑨x(s+α1) → 1/(s+α2) = P + Q (s+α1)/(s+α2)
s = -α1 → 1/(-α1+α2) = P

--- ⑨

--- ⑩

--- ⑪

⑧⑨⑩⑪→ L(Xm) = K1k3 (1/(α2-α1)) ( 1/(s+α1) - 1/(s+α2) ) L(Cp)= (K1k3/(α2-α1)) ( L(exp(-α1t) – L(exp(-α2t) ) L(Cp)
Xm = (K1k3/(α2-α1)) ( exp(-α1t) – exp(-α2t) )＊Cp(＊ convolution )基礎

脳のダイナミックFDG PET解析

X(t) = Xe(t) + Xm(t)= (K1/(α2-α1)) ((K1+k4-α1) exp(-α1t) + (–k3-k4-α2) exp(-α2t) )＊Cp(t)組織の時間放射能曲線 = 組織の応答関数＊動脈血漿の時間放射能曲線（ 脳のFDGに対する応答関数は、Bi-exponential 曲線になる。)
FDGダイナミック脳PET検査で測定されるデータは、採血による動脈
血漿中時間放射能曲線 Cp(t) とPET
画像から得られる脳内の任意の部
位における組織時間放射能曲線 X(t)である。
ラプラス変換から算出された式に測定をされた曲線を入力して脳内の速度定数 K1,k2,k3,k4の分布画像が得られる。
K1の速度でXeに流入するが、FDGを欲しがらない組織が 速度k2 でCpに戻る。そのため Xe からXm へ移動する割合が k3/(k2+k3)になるのでCpからXmへの速度定数は、
K1 x k3/(k2+k3) となる。
左図でも頭蓋骨や頭皮等の組織はブドウ糖を欲しがっていないので（脂肪酸で生きている）k2値が高いことが示されている。
Xmはブドウ糖を必要とする組織を示す。脳組織100gが1分間で何molのブドウ糖を摂取するか（CMRGlc;
Cerebral Metabolic Rate of Glucose）
を以下の式で定量計算ができる。
CMRGlc = 血糖値 x K1x k3/(k2+k3) / (Lumped constant)
Lumped constant は、ブドウ糖よりもFDGのほうが取込みが低下する比率で、脳組織では42%程度とされており、その値でFDGから算出した摂取率を補正している。

基礎

Patlak 解析

Patlak解析は、３コンパートメントモデルを簡易的に解く方法である。
３コンパートメントモデルで算出される速度定数の中で、k4の値は非常に小さい場合が多く、k4を 0 と設定したモデルで簡易的な解法を求める。

Xe’ = K1 Cp –(k2+k3) Xe --- ①sL(Xe) = K1 L(Cp) – (k2+k3) L(Xe)(s+(k2+k3)) L(Xe) = K1 L(Cp)
Xe(t) = K1 Cp(t)＊exp( -(k2+k3)t )Xm’ = k3 Xe--- ②

X = Xe + Xm--- ③（＊ convolution）

Xe ≒ K1 Cp ∫ exp(-(k2+k3)t) dt = (K1/(k2+k3)) Cp (t→∞ の場合 Cp ≒ const)①→ (k2+k3)Xe = K1Cp –Xe’②→ Xm’ = k3/(k2+k3)( K1 Cp –Xe’)= (K1k3/(k2+k3)) Cp – (k3/(k2+k3)) Xe’
Xm(t) = (K1k3/(k2+k3))∫Cp(τ) dτ – (k3/(k2+k3))Xe(t)③→ X(t) = Xe(t) + Xm(t) = (1 – k3/(k2+k3)) Xe(t) * (K1k3/(k2+k3)) ∫Cp(τ) dτ= (K1k2/(k2+k3)2) Cp(t) + (K1k3/(k2+k3)) ∫Cp(τ) dτ
ここで、Vd (分布容積) = (K1k2/(k2+k3)2) 、Ki (総合速度定数) = (K1k3/(k2+k3))とすると、X(t) / Cp = Vd + (Ki / Cp) ∫Cp(τ) dτ
18F-FDG脳PETなどでは、精度の良いダイナ

ミック断層画像が得られるので、連立微分方程式をそのまま解いて定量画像を算出す
るが、99mTc-ECDなどのSPECT薬剤では画像の精度が悪いので計算誤差を生じにくい簡便法（Patlak plot 法）を用いる。
Ki （薬剤動態の速度定数 K1, k2, k3 から求められる値で、薬剤を能動的に取り込む組織への集積速度）および Vd （K1, k2, k3 から求められる値で、薬剤の体内の分布容積
Distribution Volume を意味する）を用いて、３コンパートメントモデルの式を変形すると、左図に示すPatlak plot の式 が得られる。

基礎

123I-IMPによる脳血流定量123I-IMP は、脳SPECTに使用される薬剤で、脳組織への集積量が脳血流量に良く比例することから、脳血流検査に適する。
99mTc-ECDや 99mTc-HMPAOなどは、緊急検査には適し投与量も多くできるので脳SPECT像の画質は良いが、左図のように軽度の脳血流低下部位のカウントはあまり低下しないため、認知症などの検査には123I-IMPで脳血流を定量すべきである。
123I-IMP は、脳血流定量のためにシリンジポンプを使い正確に１分間で静脈投与する。 123I-IMP は、 肺組織に摂取されてから動脈系に入る。そのため他の薬剤と異なり、投与20分後程度まで動脈血中濃度が高い。その特徴を利用して、投与10分後に橈骨動脈などから１点動脈採血を行う。複数例の連続動脈採血検査から得られた標準動脈血漿放射能曲線があるので、 １点動脈採血をオクタノール抽出補正した血漿放射能で、各症例の動脈血漿放射能曲線 Cp(t) を導出する。
123I-IMP は脂溶性だが、検査中、肝などで代謝されて３割ほどが水溶性の放射性薬剤に変化する。水溶性の放射性薬剤は
脳組織には入らないので、その放射能を入力関数から除外したい。そのため脂溶性成分だけが溶け込むオクタノールに血漿を入れて正確に123I-IMPだけの放射能を測定する。
SPECT撮像は投与20分後から20分間実施する。
得られた脳SPECT像は、30分後の像と解釈する。
123I-IMPは、2コンパートメントモデルに従うが、k2は小さいとみなし、下記の式による概算値が脳血流定量に使用される。この123I-IMPによる脳血流定量法を、
オートラジオグラフィ法（ARG法）という。

SPECT画像の画素値はカウントなので、放射能濃度が判っている 123I 水溶液を入れたファントムをSPECT装置で20分間撮像し、さらにSPECT画像を吸収補正し（Chang補正）、SPECTカウント値（count）と放射能（Bq/ml）の関係を測定しておく必要がある（クロスキャリブレーション）。

基礎

大脳基底核の機能測定

認知症の検査には123I-IMP脳SPECTが有効だが、脳血流分布だけでは鑑別が困難な場合がある。例えば、アルツハイマー型認知症とレビー小体型認知症は、どちらも左右頭頂葉あたりに血流低下を示し、他の薬剤による鑑別が望まれていた。
123I-イオフルパン（Dat scan®）は、ドーパミン受容体に集積する薬剤で、投与３時間後に３０分間ＳＰＥCT撮像する。正常例では、基底核のドーパミン受容体に集積する。正常集積は、基底核/バックグラウンド（小脳など）比 が ６ 以上。
パーキンソン病やレビー小体型認知症では、集積低下を示す。

123I-イオフルパン

SPECT画像

大脳基底核

Basal ganglia

大脳基底核という部位は、解剖学的構造や、機能的な内容も理解が困難である。
概要は、大脳皮質の脳神経細胞の下に白質があり、その奥にある細胞核集団。
線条体 （被殻 と 尾状核）、淡蒼球、黒質、視床下核 の総称。
視床や大脳皮質、小脳へ情報を送り、円滑な運動を調整する。
周囲の大脳辺縁系（海馬、扁桃体）は、記憶や自律神経を調整する。
基底核障害 があると、パーキンソン病、チック（突発的な瞬動や発声）を生じる。

基礎

積分法、台形公式

コンパートメントモデルを解析する場合、式の中に補正したいパラメータを追加することがある。そのような場合にはラプラス変換で解けない式が出現する。
その場合には、積分法を用いて数学的に解けない式の解を導出できる。
コンパートメントモデルの式は微分方程式が連立することが多いが、これらをラプラス変換しないで、下記に示す台形公式を連立微分方程式に加えて積分し、漸化式の型にして解析する方法である。

以下に示す３コンパートメントモデルの式(1),(2),(3)は、連立微分方程式で表されている。式(4),(5),(6)で、それぞれの微分方程式を積分形に変換し、さらに上図に示した台形公式を使い、式(7),(8)を作成している。
--------------------------------------------------------------------------------------------------------------3 compartment model analysis--------------------------------------------------------------------------------------------------------------issue x(t) = x1(t) + x2(t)_________________________________________| Plasma| K1 || k3 ||||------->|x1(t)|------>|x2(t)|| Cp(t)|||||||<-------|(Bq/ml) |<------|(Bq/ml) || ( Bq/ml)| k2 || k4 |||______________||___________ _||____________|

dCp/dt = k2 x1(t) - K1 Cp(t)dx1/dt = K1 Cp(t) + k4 x2(t) - (k2+k3) x1(t)dx2/dt = k3 x1(t) - k4 x2(t)---(1)---(2)---(3)(1)---> Cp(i) = k2 Sx1(i) - K1 SCp(i)---(4)(2)---> x1(i) = K1 SCp(i) + k4 Sx2(i) -(k2+k3) Sx1(i) ---(5)(3)---> x2(i) = k3 Sx1(i) - k4 Sx2(i)---(6)( i : データ測定時間 )( SCp(i), Sx1(i), Sx2(i) : Cp(i), x1(i), x2(i) の時間0からiまでの積分値 )Sx1(i) = Sx1(i-1) + ( x1(i-1)+x1(i) )*d
Sx2(i) = Sx2(i-1) + ( x2(i-1)+x2(i) )*d( d: データ測定時間間隔 / 2 )---(7)---(8)

基礎

積分法、漸化式

以下に示す式の変形で、式(13)から x1(i) が実測された採血データCp(i)と漸化式で算出できることがわかる。x1(i) の初期値は 0 として、帰納法の考え方から
x1(0)からx1(i-1) は x1(i) を算出する際には既知なのでx1(i) は式(13)から計算できる。さらに式(10) から、x2(i) も既知データCp(i)と x2(i-1) と既知になった x(i) を用いて計算できる。 PETから得た脳などの時間放射能曲線 x(i) は、x1(i) + x2(i) と同じ曲線になるはずなので、その曲線間の差が最小になるように、４個の速度
定数 K1, k2, k3, k4 を非線形最小二乗法のプログラムを用いて導出する。
(5)(6)---> x1(i) + x2(i) = K1 SCp(i) - k2 Sx1(i)(7)(9)---> x1(i) + x2(i) = K1 SCp(i) - k2 Sx1(i-1) - k2 d x1(i-1) - k2 d x1(i)---(9)

x2(i) = K1 SCp(i) - k2 Sx1(i-1) - k2 d x1(i-1) - (1+k2d) x1(i)---(10)..................................................................................................................................
(5)(7)(8)---> x1(i) = K1 SCp(i) + k4 Sx2(i-1) + k4d x2(i-1) + k4d x2(i)- (k2+k3)Sx1(i-1) - (k2+k3)d x1(i-1) -(k2+k3)d x1(i)(1+(k2+k3)d) x1(i) = K1SCp(i) + k4 Sx2(i-1) + k4d x2(i-1) + k4d x2(i)- (k2+k3)Sx1(i-1) - (k2+k3)d x1(i-1)---(11)(6)(7)(8)--->
x2(i) = k3 Sx1(i-1) + k3d x1(i-1) + k3d x1(i) - k4 Sx2(i-1) - k4d x2(i-1) - k4d x2(i) ---(12)(11)(12)--->(1+(k2+k3)d) x1(i) = K1SCp(i) + k4Sx2(i-1) + k4dx2(i-1) - (k2+k3)Sx1(i-1) - (k2+k3)x1(i-1)+ k4d ( k3Sx1(i-1) + k3dx1(i-1) + k3dx1(i) - k4Sx2(i-1) - k4dx2(i-1) ) /(1+k4d)(1+(k2+k3)d)(1+k4d) x1(i) = K1(1+k4d)SCp(i) + k4(1+k4d)Sx2(i-1) + k4d(1+k4d)x2(i-1)- (1+k4d)(k2+k3)Sx1(i-1) - (1+k4d)(k2+k3)x1(i-1)+ k3k4d Sx1(i-1) + k3k4dd x1(i-1) + k3k4dd x1(i)- k4k4d Sx2(i-1) - k4k4dd x2(i-1)(1+(k2+k3+k4)d+k2k4dd)x1(i) = K1(1+k4d)SCp(i) - (k2+k3+k2k4d)(Sx1(i-1)+dx1(i-1))+ k4((Sx2(i-1)+dx2(i-1))
x1(i) = (K1(1+k4d)SCp(i)-(k2+k3+k2k4d)(Sx1(i-1)+dx1(i-1))+k4((Sx2(i-1)+dx2(i-1))/ (1+(k2+k3+k4)d+k2k4dd)---(13)..................................................................................................................................
(13) ---> 曲線 x1(i) は、SCp(i), Sx1(i-1), x1(i-1), Sx2(i-1), x2(i-1) から算出される。
(10) ---> 曲線 x2(i) は、SCp(i),Sx1(i-1),x1(i-1),x1(i) から算出される。
PETで得た組織時間放射能曲線 x(i) と 曲線 ( x1(i) + x2(i) ) の誤差二乗和を最小にする K1,k2,k3,k4 の値を、非線形最小二乗法( Powell法 ) で求める。

基礎

積分法、非線形最小二乗法

積分法による3コンパートメントモデル解析の具体的なC言語プログラムを示す。
関数 func_3comp は、与えられた速度定数 k1,k2,k3,k4 から得る曲線 x1(i)+x2(i) と、測定された組織放射能曲線 x(i) との誤差二乗和 error を出力する。
main関数では、error を最小にする k1,k2,k3,k4 の数値を算出する関数 Powell を用意する。Powell法は、複数の変数を持つ関数が最小値になるときの変数の値を推定する非線形最小二乗法で、多くの書籍にプログラムが掲載されている。
#include <stdio.h>#include <Powell.h>double cp[3600], x[3600], x1[3600], x2[3600] ;double K[5] ;double interval = 1.0 ;/* 関数 Powell の設定/* 曲線を入れる配列の設定/* 速度定数を入れる配列の設定/* 放射能曲線の測定間隔の設定

*/
*/
*/
*/

double func_3comp( double *k ) /* 3-compartment model analysis to calculate k1,k2,k3,k4 */{int i ;
double k1, k2, k3, k4, scp, sx1, sx2, d, error ;/* 変数の設定*/
d = interval / 2.0 ;/* 積分法の計算間隔設定 */k1 = k[1] ; k2 = k[2] ; k3 = k[3] ; k4 = k[4] ;/* 速度定数の取得*/
x1[0]=0.0 ; x2[0] = 0.0 ; scp = cp[0] ; sx1 = 0.0 ; sx2=0.0 ;/* 初期値の設定*/
for( i=1 ; i<3600 ; i++ ){/* この例では1時間(1秒から3600秒まで)の曲線データを解析 */scp = scp + ( cp[i-1] + cp[i] )*d;/* 台形公式 */x1[i] = ( k1*(1+k4*d)*scp - (k2+k3+k2*k4*d)*(sx1+x1[i-1]*d) + k4*(sx2+x2[i-1]*d) )/ (1.0 + (k2+k3+k4)*d + k2*k4*d*d ) ;/* 式(13) */x2[i] = k1*scp - k2*sx1 - k2*d*x1[i-1] - (1+k2*d)*x1[i] ;/* 式(10) */sx1 = sx1 + ( x1[i-1] + x1[i] )*d ; sx2 = sx2 + ( x2[i-1] + x2[i] )*d ; /* 台形公式 */}error =0. 0 ;/* 曲線 (x1[i]+x2[i]) と曲線 x[i] の誤差二乗和の算出 */for( i=0 ; i<3600 ; i++ ){ error = error + ( x1[i]+x2[i] – x[i] ) * ( x1[i]+x2[i] – x[i] ) ; }return error ;/* 誤差二乗和 error の出力 */}void main(void){
int i ; char f[100] ; double data ; FILE *fp ;/* 変数の設定*/

f = “血漿放射能曲線 Cp” ; fp = fopen( f,“rb”);/* 血漿放射能曲線 Cp の入力 */for( i=0 ; i<3600 ; i++ ){ fread( &data, sizeof(double), 1, fp ) ; cp[i]=data; } fclose(fp);
f = “組織放射能曲線 X” ; fp = fopen( f,“rb”);/* 組織放射能曲線 X の入力
for( i=0 ; i<3600 ; i++ ){ fread( &data, sizeof(double), 1, fp ) ; x[i]=data; } fclose(fp);K[1] = 0.1; K[2] = 0.1; K[3] = 0.1; K[4]=0.01;/* 速度定数の初期値設定*/*/

Powell( K, 4, func_3comp ) ; /* func_3comp が最小値になる４個の速度定数を算出 */printf(“ K1, k2, k3, k4 = %lf, %lf, %lf, %lf “, K[1], K[2], K[3], K[4] ) ; /* 速度定数の表示 */}基礎

Fick の原理と15O-H2O PET

15O-H O を用いてダイナミックPETを実施すると、心筋などの組織血流量を Fick2の原理から理論上正確に測定できる。Fickの原理とは、単位時間に微小組織が摂取した量を dCt/dt、動脈血中濃度をCa、静脈血中濃度をCvとし、微小組織内では流入する動脈血流量と流出する静脈血流量は等しく、それをｆとみなすと、
d Ct/dt = f (Ca - Cv) となる。

15O-H O ダイナミックPETから、心筋組織への水の流入速度定数を算出すれば、
2
その値は心筋への動脈血流量に等しい。13N-NH3や 82Rb などを用いたPETによる

心筋血流測定は、心筋組織への流入速度定数K1は心筋血流量とは比例せず、Renkin-Croneの式による曲線関係を示すので、その補正式を導いて心筋血流量を計算する必要がある。
15O-H Oを用いればFickの式から心筋血流量を算出できるが、心筋静脈血放射能
2
Cv を測定することは不可能である。そこで、組織中の血管は静脈の方が多く
内径も広いので、組織中放射能 Ct は、静脈血放射能 Cv と比例すると考える。
静脈血放射能に対する組織中放射能の比率を分配係数 p（partition coefficient）と定義し（心筋では p = 0.91 (mL/g) ）、p = Ct / Cv の式をFickの式に代入すると、
d Ct(t) / dt = f Ca(t) – (f / p) Ct(t) の式が導かれる。
左室心筋は当然、拍動している組織なので、PET画像で左室内腔や心筋壁に関心領域を設定しても正確な動脈血や心筋組織の時間
放射能曲線にはならない。
そこで、左室内腔、心筋に設定した関心領
域の時間放射能曲線を、LV(t)、R(t) とする。
求めたい動脈血と心筋の時間放射能曲線を、
Ca(t)、Ct(t) とする。
これらの曲線の関係は以下の式で示される。
R(t) = PTF Ct(t) + va Ca(t)LV(t) = b Ca(t) + (1-b)m Ct(t)心筋壁放射能 R(t) は、真の心筋組織放射能 Ct(t) の割合が PTF で、左室内腔の動脈血放射能混入の割合が va となる。PTF（ perfusable tissue fraction）は、
心筋関心領域内1mL 中での灌流のある組織が占める割合 (g/mL)である。PTFと va は未知の変数で、コンパートメントモデル解析によって算出される。
左室内腔放射能 LV(t) は、動脈血放射能 Ca(t) の割合が b で、心筋組織放射能Ct(t) の割合が (1-b)とする。b の値は 0.85 になるように左室内腔の関心領域設定プログラムを作成する（領域内平均値が領域内最大値の85%になるように関心領域の位置を設定する）。Ct(t)に m (myocardial tissue density = 1.04 g/mL)が掛かるのは、Ct(t) の単位 (Bq/g) を mCt(t) (Bq/mL) に変換するためである。

基礎

15O-H O PETによる心筋血流測定
2

15O-H O ダイナミック心筋PETで得る測定値は、心筋ROI放射能曲線 R(t)、左室
2

内腔ROI放射能曲線 LV(t)である。さらに、心筋組織分配係数 p( = 0.91 (mL/g) ) 、
心筋組織密度m( = 1.04 (g/mL) ) 、左室内腔ROIの放射能リカバリー係数 b( = 0.85)も既知である。
未知データは、心筋組織放射能曲線 Ct(t)、動脈血放射能曲線 Ca(t)、PTF、va 、そして心筋血流量ｆの値である。
これらの既知データと未知データの区別に注意しながら、積分法と台形公式を用いて、未知データを導出できる漸化式を、下記の式 (1), (2), (3) から作成する。
( i ; 測定時間 (sec) )( PTF ; 心筋ROI内の灌流組織量 (g/mL) )( va ; 心筋ROI内への動脈血放射能混入率 )( m ; 心筋組織密度 = 1.04 (g/mL))( b ; LV ROI の放射能リカバリー係数=0.85 )( Fick の原理 )( p ; 心筋組織分配係数 = 0.91 (mL/g) )R[i] = PTF*Ct[i] + va*Ca[i]----(1)

LV[i] = b*Ca[i] + (1-b)*m*Ct[i]----(2)d Ct[i]/dt = f*Ca[i] – (f/p)*Ct[i]----(3)(3)--> Ct[i] = f*SCa[i] - (f/p)*SCt[i]---(3)’( SCa[i], SCt[i] ; Ca[i], Ct[i] の 0 から i 秒までの積分 )SCa[i] = SCa[i-1] + (Ca[i-1]+Ca[i])*d ---(4)SCt[i] = SCt[i-1] + (Ct[i-1]+Ct[i])*d( 台形公式 ) ( d ; 測定時間間隔 / 2 )---(5)(3)'(4)(5) --->
Ct[i] = f*SCa[i-1]+f*d*Ca[i-1]+f*d*Ca[i]-(f/p)*SCt[i-1]-(f/p)*d*Ct[i-1]-(f/p)*d*Ct[i] ---(6)(2)(6) --> (1+(f/p)*d+f*d*m*(1-b)/b)*Ct[i]= f*SCa[i-1]+f*d*Ca[i-1]-(f/p)*SCt[i-1]-(f/p)*d*Ct[i-1]+(f/b)*d*LV[i]---(6)’(6)’ --> Ct[i] = (f*SCa[i-1]+f*d*Ca[i-1]-(f/p)*SCt[i-1]-(f/p)*d*Ct[i-1]+(f/b)*d*LV[i])/(1+(f/p)*d+f*d*m*(1-b)/b)---(7)(2) ---> Ca[i] = (LV[i]-(1-b)*m*Ct[i])/b---(8)

漸化式 (4)(5)(7)(8) から 曲線 Ct[i], Ca[i] が算出される。
(1) --> 曲線 ( PTF*Ct[i] + va*Ca[i]) と 曲線 R[i] の誤差二乗和が最小になる
f, PTF, va の値を、非線形最小二乗法（Powell法）で算出する。

基礎

積分法による心筋血流量算出

Fickの原理に基づくシングルコンパートメントモデル解析の積分法によるC言語
プログラムを示す。関数 func_MBF は、与えられた変数 f, va, PTF から得る曲線
PTF*Ct(i)+va*Ca(i) と、心筋ROI放射能曲線 R(i) との誤差二乗和 error を出力する。
main関数では、error を最小にする3個の変数 f, va, PTF の値を算出する関数
Powell を用意する。
#include <stdio.h>#include <Powell.h>double R[360], LV[360], Ct[360], Ca[360] ;double K[4] ;double interval = 1.0 ;/* 関数 Powell の設定/* 曲線を入れる配列の設定/* 変数を入れる配列の設定/* 放射能曲線の測定間隔の設定

*/
*/
*/
*/

double func_MBF( double *k ) /* 1-compartment model analysis to calculate MBF*/{
int i ; double f, va, PTF, SCa, SCt, d, z, error, p=0.91, b=0.85, m=1.04 ; /* 変数の設定 */f = k[1] ; va = k[2] ; PTF = k[3] ;/* 変数 f, va, PTF の取得 */Ct[0] = 0.0 ; Ca[0] = 0.0 ; SCt = Ct[0] ; SCa = Ca[0] ;/* 初期値の設定*/
d = interval / 2.0 ;/* 積分法の計算間隔設定 */for( i=1 ; i<360 ; i++ ){/* この例では6分間 (1秒から360秒まで)の曲線データを解析 */z = (1.0+f*d/p+f*d*m*(1.0-b)/b) ; if(z<1.0e-20)z=1.0e-20;/* 0での割算防止 */Ct[i] = ( f*SCa+f*d*Ca[i-1]-SCt*f/p-Ct[i-1]*f*d/p+LV[i]*f*d/b ) / z ;/* 式 (7) */Ca[i] = ( LV[i]-m*(1.0-b)*Ct[i] ) / b;/* 式 (8) */SCa = Sca + ( Ca[i-1] + Ca[i] )*d;SCt = SCt + ( Ct[i-1] + Ct[i] )*d;/* 式 (4) *//* 式 (5) */}error =0.0 ;/* 曲線 ( PTF*Ct[i]+va*Ca[i] ) と曲線 R[i] の誤差二乗和の算出 */for( i=0 ; i<360 ; i++ ){ error = error + ( PTF*Ct[i]+va*Ca[i] – R[i] ) * (PTF*Ct[i]+va*Ca[i] – R[i] ) ; }return error ;/* 誤差二乗和 error の出力 */}void main(void){
int i ; char f[100] ; double data ; FILE *fp ;/* 変数の設定*/

f = “心筋ROI放射能曲線 R” ; fp = fopen( f,“rb”);/* 心筋ROI放射能曲線 R の入力 */for( i=0 ; i<360 ; i++ ){ fread( &data, sizeof(double), 1, fp ) ; R[i]=data; } fclose(fp);
f = “左室内腔放射能曲線 LV” ; fp = fopen( f,“rb”); /* 左室内腔放射能曲線 LV の入力 */for( i=0 ; i<360 ; i++ ){ fread( &data, sizeof(double), 1, fp ) ; LV[i]=data; } fclose(fp);K[1] = 1.0 ;K[2] = 0.3 ;Powell( K, 3, func_MBF ) ;K[3] = 0.8 ;/*変数 f, va, PTF の初期値設定 *//* func_MBF が最小値になる 3個の変数を算出 */printf(“ MBF , va , PTF = %lf, %lf, %lf “, K[1], K[2], K[3] ) ;}/*変数 f, va, PTF の表示 */13NH 11C-acetate, 82Rbによる
3,
基礎

心筋血流定量

放射性アンモニア 13NH3 は、初回循環で心筋への取込率が80%と高く、心筋細胞内でグルタミンのアミノ基に代謝固定される。正確には2コンパートメン
トモデルに従うが、取込率が高いので、２分間のダイナミックPET画像を用いてマイクロスフェアモデルで簡便に以下の式で心筋血流量を算出できる。

11C-acetateや82Rb は、心筋への取込率が高くないので、左図のように2コンパートメントモデルで、動脈から心筋組織への流入速度定数K1を算出し、心筋血流量MBFをRenkin-Crone の式から導く。
我々の施設で15O-H2O PETから求めたMBFと20分
間ダイナミック収集した11C-acetate のPETから算出したK1の関係式は、
K1 = MBF (1- 0.764 exp( -1.001 / MBF))であった(1)。
また、15O-H2O PETから求めたMBFと8分間ダイナ
ミック収集した82Rb PETから算出したK1の関係は以下の式であった(2)。
K1 = MBF (1- 0.860 exp( -0.543 / MBF))
これらのRenkin-Crone の式にK1の値を代入してもMBFは単純には算出できないので、Powell法やsimplex法などを用いて、非線形最小二乗法を利用したプログラムを作成してMBFを求める必要がある。
(1) Improved spillover correction model to quantify
myocardial blood flow by 11C-acetate PET: comparison
with 15O-H2O PET.
Mori Y, Manabe O, Naya M, Tomiyama Y, Yoshinaga K,
Magota K, Oyama-Manabe N, Hirata K, Tsutsui H, Tamaki
N, Katoh C. Ann Nucl Med. 2015 Jan;29(1):15-20.
(2) Quantification of regional myocardial blood flow
estimation with three-dimensional dynamic rubidium-82
PET and modified spillover correction model.
Katoh C, Yoshinaga K, Klein R, Kasai K, Tomiyama Y,
Manabe O, Naya M, Sakakibara M, Tsutsui H, deKemp RA,
Tamaki N. J Nucl Cardiol. 2012 Aug;19(4):763-74.

基礎

２コンパートメントモデル
による心筋血流測定

11C-acetateや82Rbダイナミック心筋PETで得る測定値は、心筋ROI放射能曲線 R(t)、左室内腔ROI放射能曲線 LV(t)である。心筋組織密度m( = 1.04 (g/mL) ) 、左室内腔
ROIの放射能リカバリー係数 b( = 0.85) は既知である。
未知データは、心筋組織放射能曲線 Ct(t)、動脈血放射能曲線 Ca(t)、PTF、va 、そして心筋への流入速度定数 K1 と流出速度定数 k2 の値である。
これらの既知データと未知データの区別に注意しながら、積分法と台形公式を用いて、未知データを導出できる漸化式を、下記の式 (1), (2), (3) から作成する。
( i ; 測定時間 (sec) )( PTF ; 心筋ROI内の灌流組織量 (g/mL) )( va ; 心筋ROI内への動脈血放射能混入率 )( m ; 心筋組織密度 = 1.04 (g/mL))( b ; LV ROI の放射能リカバリー係数=0.85 )( 2 compartment model )R[i] = PTF*Ct[i] + va*Ca[i]----(1)

LV[i] = b*Ca[i] + (1-b)*m*Ct[i]----(2)d Ct[i]/dt = K1*Ca[i] – k2*Ct[i]----(3)(3)--> Ct[i] = K1*SCa[i] – k2*SCt[i]---(3)’( SCa[i], SCt[i] ; Ca[i], Ct[i] の 0 から i 秒までの積分 )SCa[i] = SCa[i-1] + (Ca[i-1]+Ca[i])*d ---(4)SCt[i] = SCt[i-1] + (Ct[i-1]+Ct[i])*d( 台形公式 ) ( d ; 測定時間間隔 / 2 )---(5)(3)'(4)(5) --->
Ct[i] = f*SCa[i-1]+f*d*Ca[i-1]+f*d*Ca[i]-(f/p)*SCt[i-1]-(f/p)*d*Ct[i-1]-(f/p)*d*Ct[i] ---(6)(2)(6) --> (1+k2*d+K1*d*m*(1-b)/b)*Ct[i]= K1*SCa[i-1]+K1*d*Ca[i-1]-k2*SCt[i-1]-k2*d*Ct[i-1]+(K1/b)*d*LV[i] ---(6)’(6)’ --> Ct[i] = (K1*SCa[i-1]+K1*d*Ca[i-1]-k2*SCt[i-1]-k2*d*Ct[i-1]+(K1/b)*d*LV[i])/(1+k2*d+K1*d*m*(1-b)/b)---(7)(2) ---> Ca[i] = (LV[i]-(1-b)*m*Ct[i])/b---(8)

漸化式 (4)(5)(7)(8) から 曲線 Ct[i], Ca[i] が算出される。
(1) --> 曲線 ( PTF*Ct[i] + va*Ca[i]) と 曲線 R[i] の誤差二乗和が最小になる
K1, k2, PTF, va の値を、非線形最小二乗法（Powell法）で算出する。

積分法による２コンパートメント
モデル心筋血流量算出

基礎２コンパートメントモデル解析の積分法によるC言語プログラムを示す。関数
func_2comp は、与えられた変数 K1, k2, va, PTF から得る曲線 PTF*Ct(i)+va*Ca(i) と、心筋ROI放射能曲線 R(i) との誤差二乗和 error を出力する。
main関数では、error を最小にする3個の変数 K1, k2, va, PTF の値を算出する関数
Powell を用意する。
#include <stdio.h>#include <Powell.h>double R[1200], LV[1200], Ct[1200], Ca[1200] ;double K[4] ;double interval = 1.0 ;/* 関数 Powell の設定/* 曲線を入れる配列の設定/* 変数を入れる配列の設定/* 放射能曲線の測定間隔の設定

*/
*/
*/
*/

double func_comp2( double *k ) /* 2-compartment model analysis to calculate K1*/{
int i ; double K1, k2, va, PTF, SCa, SCt, d, z, error, b=0.85, m=1.04 ; /* 変数の設定 */K1 = K[1] ; k2 = K[2] ; va = K[3] ; PTF = K[4] ;/* 変数 K1,k2, va, PTF の取得 */Ct[0] = 0.0 ; Ca[0] = 0.0 ; SCt = Ct[0] ; SCa = Ca[0] ;/* 初期値の設定*/
d = interval / 2.0 ;/* 積分法の計算間隔設定 */for( i=1 ; i<1200 ; i++ ){ /* この例では20分間 (1秒から1200秒まで)の曲線データを解析 */z = (1.0+k2*d+K1*d*m*(1.0-b)/b) ; if(z<1.0e-20)z=1.0e-20;/* 0での割算防止 */Ct[i] = ( K1*SCa+K1*d*Ca[i-1]-k2*SCt-k2*d*Ct[i-1]+K1*(d/b)*LV[i] ) / z ; /* 式 (7) */Ca[i] = ( LV[i]-m*(1.0-b)*Ct[i] ) / b;/* 式 (8) */SCa = Sca + ( Ca[i-1] + Ca[i] )*d;SCt = SCt + ( Ct[i-1] + Ct[i] )*d;/* 式 (4) *//* 式 (5) */}error =0.0 ;/* 曲線 ( PTF*Ct[i]+va*Ca[i] ) と曲線 R[i] の誤差二乗和の算出 */for( i=0 ; i<1200 ; i++ ){ error = error + ( PTF*Ct[i]+va*Ca[i] – R[i] ) * (PTF*Ct[i]+va*Ca[i] – R[i] ) ; }return error ;/* 誤差二乗和 error の出力 */}void main(void){
int i ; char f[100] ; double data ; FILE *fp ;/* 変数の設定*/

f = “心筋ROI放射能曲線 R” ; fp = fopen( f,“rb”);/* 心筋ROI放射能曲線 R の入力 */for( i=0 ; i<1200 ; i++ ){ fread( &data, sizeof(double), 1, fp ) ; R[i]=data; } fclose(fp);
f = “左室内腔放射能曲線 LV” ; fp = fopen( f,“rb”); /* 左室内腔放射能曲線 LV の入力 */for( i=0 ; i<1200 ; i++ ){ fread( &data, sizeof(double), 1, fp ) ; LV[i]=data; } fclose(fp);K[1] = 1.0 ;K[2] = 0.1 ;K[3] = 0.2 ; K = 0.8 ;Powell( K, 4, func_comp2 ) ;/* 変数 K1, k2, va, PTF の初期値設定 *//* func_comp2 が最小値になる 4個の変数を算出 */printf(“ K1 , k2, va , PTF = %lf, %lf, %lf, %lf “, K[1], K[2], K[3], K[4] ) ; /*変数 K1, k2, va, PTF の表示 */}基礎

Renkin-Crone 式の解法

２コンパートメントモデル解析から心筋への薬剤流入速度K1が得られたら、
Renkin-Crone の式で心筋血流量MBFを算出するが、K1 = MBF (1- α exp( β / MBF))という形式なので直接MBFを算出できない。そこで、解けない式の近似解を導く方法として、Powell法などの最小二乗法を、ここでも利用する。

#include <stdio.h>#include <Powell.h>/* 関数 Powell の設定*/

double K[2] , K1;/* 変数を入れる配列, 変数の設定*/

double Renkin-Crone( double *K ){double a, b, MBF, error ;/* Renkin-Crone to calculate MBF from K1 *//* 変数の設定

MBF = K[1];*//* MBF値の受取り */a = 0. 764 ; b = -1.001 ; /* 変数 α, β の設定（ここでは11C-acetate PET の変数を設定) */error = K1 – MBF * ( 1.0 – a * exp( b / MBF ) ) ;return error * error ;}void main(void){
double MBF ;printf(“ Input K1 value = “); scanf(“%lf”, &K1 );

}

変数の設定

*/

/* K1 値 の入力*//*

K[1] = 1.0 ;/* MBF値が渡される変数 K[1] の初期値設定 */Powell( K, 1, Renkin-Crone ) ;/* Renkin-Crone が最小値になる変数 MBF を算出 */MBF = K[1] ;printf(“ MBF = %lf “, MBF ) ;/*Renkin-Crone が最小値になる変数 MBF を受取る *//*変数 MBF の表示 */基礎

Reference tissue model

11C-フルマゼニルは、大脳皮質脳神経細胞のベンゾジアゼピン受容体に特異的に結合する薬剤で、定量したい値は大脳皮質ベンゾジアゼピン受容体結合能である。このように特定の受容体に結合する薬剤のPET検査では、受容体結合能
BP（Binding Potential）が定量できればよい場合がある。そのような検査では、
動脈採血が不要な便利な解析法があり、Reference tissue model という。
ただし、分布容積 DV（distribution volume）の定量には動脈採血が必要である。
受容体結合能 BPは
BP = k3 / k4の式で算出される。

11C-フルマゼニルでは、受容体がない組織（Reference tissue）として脳幹部の橋 ponsに関心領域を設定し、Reference tissue 曲線を入力関数として使用する。
ただし単純に動脈血放射能曲線として代用できないので、Reference tissue 曲線
を入力関数にすることができる式を下図のように作成する。

受容体のある組織での K1 / k2 比と受容体のない組織
での K1 / k2 比（K1’ / k2’）は同じという式K1’ / k2’ = K1 / k2
を使用する点がこの解析法のポイントである。

基礎

積分法によるReference tissue model

Reference tissue model を、積分法を用いて解析する方法を示す。
既知曲線データ（ダイナミックPET画像から得る放射能曲線）
Ct(t) ; (= Cf(t) + Cb(t) ) 組織放射能曲線
Cr(t) ; reference tissue の放射能曲線
未知曲線データ（漸化式から算出する曲線）
Cf(t) ; 組織内の free tissue （受容体のない組織）の放射能曲線
Cb(t) ; 組織内の bound tissue （受容体のある組織）の放射能曲線
未知変数
測定対象組織の 3コンパートメントモデル速度定数 K1, k2, k3, k4
Reference組織の2コンパートメントモデル速度定数 K1’ , k2’
K1 / K1’ 比 R1 （ R1 = K1 / K1’）
dCr(t)/dt = K1’Cp(t) - k2’Cr(t)dCf(t)/dt = K1 Cp(t) – k2 Cf(t) – k3 Cf(t) + k4 Cb(t)dCb(t)/dt = k3 Cf(t) – k4 Cb(t)R1 = K1/K1‘ ----(4)K1'/k2'= K1/k2 ----(5)----(1)----(2)----(3)(4)(5) ---> k2'= k2/R1 ----(6)式 (1)から(6) から、曲線 Cf(t), Cb(t) を算出する以下の漸化式が導かれる。
式 (7), (8), (9), (10), (11) の順序で解いていけば、左辺が未知数で、右辺が既知数で構成された式になっている。
また、漸化式の中に動脈血漿曲線 Cp(t) が表れていないことに注目する。
つまり、動脈採血が不要な漸化式になっている。
SCr = SCr + ( Cr(t-1)+Cr(t) )*d ;----(7)( 台形公式 ) ( d ; 測定時間間隔 / 2 )( SCr ; Cr(t) の 0 から t 秒までの積分値 )Cf(t) = (1.0+k4*d)*(R1*Cr(t)+k2*CSr)-(k2+k3+k2*k4*d)*SCf+Cf(t-1)*d)+k4*(SCb+Cb(t-1)*d)/ ( 1.0+(k2+k3+k4+k2*k4*d)*d) ;----(8)Cb(t) = R1*Cr(t)+k2*SCr-k2*SCf-k2*d*Cf(t-1)-(1.0+k2*d)*Cf(t) ;SCf = SCf + ( Cf(t-1) +Cf(t) )*d ; ----(10)SCb = SCb + ( Cb(t-1)+Cb(t) )*d; ----(11)----(9)( SCf ; Cf(t) の 0 から t 秒までの積分値 )( SCb ; Cb(t) の 0 から t 秒までの積分値 )漸化式 (7), (8), (9), (10), (11) から、曲線 Cf(t), Cb(t) が算出される。

曲線 ( Cf(t) + Cb(t) ) と 曲線 Ct(t) の誤差二乗和が最小になる R1, k2, k3, k4 の値を、非線形最小二乗法（Powell法）で算出する。
受容体結合能 BP は、BP = k3/k4 から算出される。
Reference tissue model は、K1 は算出されず、動脈血流量や分布容積は測定できない。

基礎

Reference tissue model のプログラム

積分法による Reference tissue model 解析の具体的なC言語プログラムを示す。
Reference tissue model for D1,D2 receptor, dopamine transporter, FMZ#include <stdio.h>#include <Powell.h>double Ct[5400], Cf[5400], Cb[5400], Cr[5400] ;double K[5] ;double interval = 1.0 ;/* 関数 Powell の設定*//* 曲線を入れる配列の設定*//* 速度定数を入れる配列の設定 *//* 放射能曲線の測定間隔の設定 */double func_Reference( double *k ) /* Reference tissue model to calculate R1,k2,k3,k4 */{int i ;
double R1, k2, k3, k4, SCr, SCf, SCb, d, error ;/* 変数の設定*/
R1 = k[1] ; k2 = k[2] ; k3 = k[3] ; k4 = k[4] ;/* 速度定数の取得 */Cf[0]=0.0 ; Cb[0] = 0.0 ; SCf = Cf[0] ; SCb = Cb[0] ; SCr = Cr[0] ; /* 初期値の設定 */d = interval / 2.0 ;/* 積分法の計算間隔設定 */for( i=1; i<5400; i++ ){ /* この例では90分間(1秒から5400秒まで)の曲線データを解析 */SCr = SCr + ( Cr[i-1]+Cr[i] )*d ;Cf[i] = (1.0+k4*d)*(R1*Cr[i]+k2*CSr)-(k2+k3+k2*k4*d)*SCf+Cf[i-1]*d)+k4*(SCb+Cb[i-1]*d)/ ( 1.0+(k2+k3+k4+k2*k4*d)*d) ;Cb[i] = R1*Cr[i]+k2*SCr-k2*SCf-k2*d*Cf[i-1]-(1.0+k2*d)*Cf[i] ;SCf = SCf + ( Cf[i-1]+Cf[i] )*d ;SCb = SCb + ( Cb[i-1]+Cb[i] )*d ;}error =0. 0 ;/* 曲線 (Cf[i]+Cb[i]) と曲線 Ct[i] の誤差二乗和の算出 */for( i=0 ; i<5400 ; i++ ){ error = error + ( Cf[i]+Cb[i] –Ct[i] ) * ( Cf[i]+Cb[i] – Ct[i] ) ; }return error ;/* 誤差二乗和 error の出力 */}void main(void){
int i ; char f[100] ; double data ; FILE *fp ;/* 変数の設定*/

f = “組織放射能曲線 Ct” ; fp = fopen( f,“rb”);/* 組織放射能曲線 Ct の入力 */for( i=0 ; i<5400 ; i++ ){ fread( &data, sizeof(double), 1, fp ) ; Ct[i]=data; } fclose(fp);

f = “Reference放射能曲線 Cr” ; fp = fopen( f,“rb”); /* Reference組織放射能曲線 Cr の入力
for( i=0 ; i<5400 ; i++ ){ fread( &data, sizeof(double), 1, fp ) ; Cr[i]=data; } fclose(fp);K[1] = 1.0; K[2] = 0.1; K[3] = 0.1; K[4]=0.01;/* 速度定数の初期値設定*/*/

Powell( K, 4, func_Reference ) ; /* func_Reference が最小値になる４個の速度定数を算出 */printf(“ R1, k2, k3, k4 = %lf, %lf, %lf, %lf “, K[1], K[2], K[3], K[4] ) ; /* 速度定数の表示
printf(“ Binding Potential BP = %lf “, K[3] / K[4] ) ;}*//* BP (= k3/k4) の表示 */