1995年[ 技術開発研究助成 ] 成果報告 ： 年報第09号
相関スペクトル解析法による局部微小網膜電位の計測
概要
1.まえがき
網膜電位(ERG)は光刺激に対して網膜細胞層から発生する電位で,臨床眼科において診断及び治療評価に応用されるとともに,視覚機能の初期過程に関連するものとして,生体工学の立場からも関心がもたれている。現在,臨床眼科におけるERG計測にはフラッシュ刺激や市松模様のパターン刺激が使われているが,網膜全体を刺激しているため,局部的特性を知ることができない。このような状況に対し最近,ランダムドットパターン刺激により局部網膜電位を計測する方法が開発され(1),またERG計測結果から網膜機能の解析が試みられている(2)が,網膜部位における光刺激ERGの変換過程の全貌が解明されるところまでは至っていない。
このような問題に一石を投じるために,発光ダイオード(LED)を用いて網膜をランダム刺激し,発生するERGの周波数特性を相関解析によって求める微小網膜電位(BERG)計測法を開発した(3)。また,微小網膜電位はフラッシュERGについて知られているa波,b波,OP波等の原波形であることを示すとともに,最小2乗法を用いて網膜部位による各成分の特徴パラメータの分布状態を解析した(4)一(7)。しかし,a波,b波,OP波等は臨床眼科において便宜上に付けられた名称で,網膜各部位からどれほどの波形成分が発生しているのか,明かになっていない。しかし一方において,解析において多数の成分を考慮することは,それに応じて各成分を特徴付ける特徴パラメータの増加をきたし,これを計測された微小網膜電位の周波数特性から推定することが困難になる。
そこで,本研究では多解像度変換を用いて前処理を行った微小網膜電位に対して最小2乗法を適用することにより特徴パラメータを推定する方法を開発した。すなわち多解像度変換によって微小網膜電位の各成分波形は,不完全ながら対応するチャンネルに分離されるので,チャンネルごとに最小2乗法を適用すれば,特徴パラメータの推定が容易になると期待される。なお,この手法を適用するためには,従来周波数領域で行っていた相関スペクトル解析を等価的に時間領域で行う必要があり,解析手法を変更した(8)。
2.微小網膜電位の解析法
2.1刺激方法
まず図1に示すような実験システムを作成した。コンピュータからランダム刺激信号を発生し,V/F変換器(1V/1KHz)によってFM変調してLEDを発光させる。原刺激信号はDC成分6Vに振幅4Vのランダム信号を重畳しているので,2V以下にはならない。したがって,FM変調後の周波数は2KHz以上であり,LED光はこの周波数で点滅しているが,網膜はこのような高周波成分に応答できないので,網膜はランダム刺激信号に比例した輝度で刺激されていると考えられる。
また,商用周波数とともに,コンピュータやLED駆動回路等から発生するノイズを防止するために,今回新たに総銅板製のシールド暗室を作成し,刺激光を光ファイバを通してシールド内の被験者の眼前300mmの位置に導入する方法をとった。ファイバは8個所(網膜中心および周辺7個所)に配置し,その単独,あるいは任意の2個所を選択刺激できるようにした。なお,ファイバの発光端面と試験眼からの距離より,眼球内での光散乱を無視すれば,網膜照射面の直径は1deg程度になる。なお,周辺刺激時には刺激光とは異なる色の微弱なLEDで正面に注視点を作った。今回は刺激光は赤色,注視点は緑色である。
2.2計測方法
電極としてERG用コンタクトレンズ電極を試験眼(右眼)の角膜上に装着し,リットマン電極を不関電極として額中央,接地電極として試験眼側の耳たぶに装着する。試験眼には散瞳剤を点眼し,シールド暗室に入り約30分暗順応した後,計測を行った。
ERGをコンピュータに取込む際のサンプル周波数は1KHzとし,サンプル定理に基づいて,カットオッフ500Hzのバターワース型2次フィルタを2段前置して高周波成分を除去した。また眼球運動に伴う筋電位など低周波のドリフトを除去するために,0.2KHzを折点周波数とする不完全微分回路を挿入し,コンピュータに取込んだ後に逆周波数特性のディジタルフィルタで波形再生する方法を採った。計測装置の総合増幅率は0.2KHzにおいて最大100,000倍で調整できる。
なお,1回の計測時間は6秒間で,再現性のある疑似白色ノイズ刺激を用いて30回程度繰返し刺激し,得られた応答を加算平均した。この6000点の時系列データを基に以下に述べる相関分析を行った。
2.3相関スペクトル解析
相関スペクトル解析はノイズに埋没した信号を検出する周知の手法であるが,その原理を簡単に説明する。
刺激をu(t),電極電位をy(t),ノイズをn(t)とし,ERG発生機構の伝達特性をg(t)とすれば,電極電位は
と表される。両辺にu(σ一λ)を掛けて期待値をとれば
となる。ここで,疑似白色ノイズ刺激とノイズとは相関がなく,両者間の相互相関関数Φun(λ)は0に収束すると期待されるので,刺激,電極間の相互相関関数は刺激の自己相関関数とERG発生機構の特性との畳込み積分で近似できる。
従来は,(2)式をフーリエ変換することにより,ERG発生機構の伝達特性G(jω)を
として求めていた。すなわち,疑似白色ノイズ刺激とERGとの相互スペクトル密度関数Φxy(jω)の疑似白色ノイズ刺激の自己スペクトル密度関数Φxy(jω)に対する比より,ERGの周波数特性を求めた。また,5種類の成分を仮定し,各成分がそれぞれむだ時間をもつ2次振動波形
で近似できるものと仮定して,その合成波形により微小網膜電位のモデルを作成し,計測データとの周波数領域(Bode線図)における誤差を最小化することにより,各成分の特徴パラメータPn=(Κn,τn,ζn,Τn)を推定していた。
今回は多解像度変換を利用するために,上記の相関スペクトル解析を以下のように時間領域において行った。また,成分の数を増してより精密な解析を行った。
3. 多解像度変換による解析
時間領域における相関スペクトル解析の概念図を図2に示す。ここでCは相関演算,LSMは最小2乗法を表している。すなわち,長さ2Lのランダム信号uを用いて計測した微小網膜電位yとこのランダム信号の前半の長さLの部分U*との相互相関関数φ。yを求めれば,その結果は(3)式のΦ。yGω)と対応することが知られている。ここでランダム信号の前半部分の信号u*を用いるのは有限長データにおいて発生する折り返し現象を避けるための技法である。同様にランダム信号u,u*の自己相関関数より(3)式のΦuuGω)に対応するφuuが得られる。そこでこのφ。、を入力として,(4)式で与えられる微小網膜電位モデルGm(S,p,)の応答φ、mを求める。ここで計測した微小網膜電位G(s)とモデルGm(s,p。)が一致していれば,φuyとφ。mとは等しくなる。したがってφ。yとφ。mとの誤差を最小化することにより特徴パラメータp.が推定できる。
以上は2.3に述べた相関スペクトル法を時間領域で行ったことに相当するが,このφ。y,φ。mをそれぞれ多解像度変換した上で最小2乗法を適用することにより成分の特徴パラメータを推定した。多解像度変換は線形変換であるから,数学的には多解像度変換によって最小2乗推定の精度が向上することはない。しかし,微小網膜電位の成分はむだ時間をもつ振動波形で,いま簡単のために1つの成分波形を考えた場合,計測波形の第n波頭とモデルの第n±1波頭が一致するような時間シフトが淀こっても2乗誤差値が減少する。すなわち,最小2乗手定の評価関数は多峰性をもち,初期値の設定によって言った局所解へ収束する可能性がある。適正な初期値をπログラムによって自動的に設定することは極めて困鮪ディスプレイに表示された波形を観察して設定する必aがあるが,複数の成分が混合した微小網膜電位の場創はこれも容易ではない。多解像度変換を行えば,各成分が周波数の異なるチャンネルに不完全ながら分離され尾ので,観察による初期値設定が容易になることが期待冤きる。φ(t)の多解像度変換ψ(σ,d)は
で与えられるが,これは一種の並列フィルタと考えられ,dは各フィルタの解像度(分析周波数)を決めるパラメータである。代表的な重み関数wとしては矩形パルス状のHarr関数があり,パラメータdを1/2ずつ減少して解像度を上げながら分析する。重み関数にはdの異なる2つの重み関数間に直交性が要求され,Harr関数はこの直交関数の中で計算が容易な利点をもつものであるが,位相のシフトを伴うので,ここでは図3に示すようなHarr関数を対称化した関数を考え,関数幅dを1/2つつ縮小した複数のチャンネルに分析した。この重み関数は連続関数では直交性をもっているが,数値計算のためにサンプルした場合,負値をとるサンプル点数2nに対し,正値をとるサンプル点数が原点の1個だけ多くなるために直交性が失われる。そこで,原点における関数値を0とすることにより直交性を保存した。
4.成分分析
この多解像度変換後の各データΨiuyについて,それぞれ分析区間をとり,その区間でのΨiuy とΨiumとの2乗誤差の和
を評価関数とし,最小2乗法を用いて成分分析を行った。
最小2乗法に従えば,(4)式のモデルにおける特徴パラメータの修正量dp,は,評価関数Jのp,による偏微分
で与えられる。これに(6)式を代入すれば
となるが, は
と表される。すなわち,モデルの伝達関数Gm(s,pn)をpnで偏微分したフィルタをプログラムし,これに多解像度変換したランダム信号を入力すれば∂Ψium / ∂Pnが得られる。また, ∂Ψium / ∂Pnと誤差(Ψiuy - Ψium)との相関演算より特徴パラメータの修正量が求められる。
なお,評価関数計算における誤差2乗積分の上下限の設定については,計測データ(あるいはランダム信号との相互相関関数φuy)のどの区間を有意とみなすかということに帰着する。フラッシュERGにおける成分波形は刺激後100ms以内に含まれているので,これを目安に100msまでの区間とした。また,波形成分の数の設定についても現時点では必然的な根拠はなく,従来はフラッシュERGの成分数より5個としていた。今回の多解像度変換を用いた成分分析では従来より詳細な解析ができるため,5個の成分数では計測データとモデルとの充分な一致が見られない場合が多く,微小網膜電位にはさらに多数の成分が含まれていることが示唆された。本来,網灰刺激部位に含まれる多数の細胞がそれぞれが独自の活動をしており,その細胞数の成分が必要ということにもなるが,現実にはそれが何個に種別できるかということであり,成分の増加によって最小2乗推定に要する演算時間が急激に増加するため,今回の解析では9個に設定した。これらの成分の内には後述のように網膜部位に局限して発生しているものと,網膜の比較的広い領域から共通して発生しているものとがあり,後者の普遍的成分に限れば,やはり5成分程度であるように思われる。
5.特徴パラメータの位置連続性の判定
網膜水平軸に沿った何点かの部位で計測した微小網膜電位から推定した特徴パラメータには後述のように局在的なものと,比較的広い領域に一定の傾向で変化しているものがあるようにみえる。しかし,隣接した部位間でどの成分同士が連続しているか,データを観察するだけで判断するのは容易ではない。そこで,隣接部位間について以下のような評価関数を考え,連続性の判定の目安にした。
特徴pn(たとえばむだ時間)について,部位kから発生したERGのi番目の成分p。kiと隣接部位k+1から発生したERGのj番目の成分p。(k+1)jとの2乗誤差の和
を評価関数とし,部位kとk+1との間の隣接行列を作る。ここでw,は特徴パラメータについての荷重である。iを固定したときの最小値を与えるj(部位kから見て最も近い成分番号)を探し,隣接行列にマークする。同様にjを固定したときの最小値を与えるi(部位k+1から見て最も近い成分番号)をマークする。これらのマークが一致したときには,それは最も連続性の高い組合わせと考える。また,マークが全く付かないものの間には連続性がないと判断する。隣接行列から連続性の高い組合わせを除き,行,列のどちらか一方にマークが付いたものが残っていれば,その成分について同様の評価を行い,つぎに連続性のある組合わせを判定する。この判定結果は荷重の値によって左右されるが,特徴パラメータの間に連続性をもつものができるだけ多くなるように調整することで,連続性の判定を行った。また,部位iのERGに近接した特徴をもつ2成分A,A'が含まれているとき,その一方Aが部位jと接続関係をもつと判断されるとA'が孤立する可能性がある。この場合,特徴が部位iにおいてA,A'に分岐していると考えるべきかも知れないが,現時点では特徴パラメータはそれぞれ1本の鎖で繋がれると仮定して解析を行った。
6.解析結果
解析例として網膜中心部刺激に対する微小網膜電位(実線)とモデル応答(破線)の多解像度変換パターンの例(網膜中心部刺激時)を図4に示す。
今回の解析では多解像変換のチャンネル数は8とし,最低次のチャンネルの荷重関数幅4dを1024点として,次数の増加に伴い荷重関数幅を1/2ずつ縮小し,解像度を上げた。したがって,多解像度変換パターンは次数が上がる程,原点付近の状況を拡大して表示している。最高次のチャンネルには最も周波数の高い成分波形のみが抽出されていると考えられ,これをモデルで近似することにより,この成分の特徴パラメータが推定される。より低い次数のチャンネルには複数の成分が含まれているが,周波数の高い成分がすでに推定されているので,他の成分の特徴パラメータも比較的容易に推定される。計算上は各チャンネルについて同時に評価関数の最小化を行っているが,定性的には以上のように高次のチャンネルから推定が進行していくように思われる。したがって,低次数のチャンネルに見られる細かいリップル様の変動まで近似することができる。このような精度の高い近似は従来の手法では得られなかったもので,多解像度変換の効果を示すものである。
一方,表1は最小2乗法で推定された特徴パラメータの接続関係を5.に述べた手法により解析した結果である。ここで,fは中心周波数,ζは減衰係数,τはむだ時間である。また反応強度は,減衰係数が小さい成分では直流ゲイン係数Kのみでは表現できないので,1V*secの仮想インパルス刺激に対する応答の最大振幅を,
により計算し,これを反応強度gmとした。刺激部位は中心部(CENTER),こめかみ側10,20,45deg(T10,T20,T45)および鼻側10,20,30,45deg(N10,N20,N30,N45)の8部位である。各部位において網膜電位は一旦9個の成分に分離されるが,5.に述べた隣接行列を適用すると12成分に分けられた。隣接行列法のもつ前記の問題点も考慮して連続性の高い(全網膜に共通する)6成分を示すと表1のようになった。
これらの他に中心周波数20～25HzでT45を除く網膜部位で発生している成分がみられたが,その他は以下のように局在性か,特徴の分岐による断片であった。すなわち,中心周波数30～40Hz発生部位T20,T45,中心周波数60～75Hz発生部位CENTERからT45までのこめかみ側,中心周波数220Hz発生部位CENTER,N10が見られた。
表1に示した6成分は網膜全野に共通する成分で,その周波数帯域より1はc波,III,IV波はb波,Vはa波,VI波はOP波と対応すると考えられる。これらの結果は先の報告(4)の結果にも矛盾しないが,II波はb波としてはやや周波数帯域が低いように思われる。フラッシュERGの対応関係は明かでないが,今回の解析による新しい知見ともいえる。
7.むすび
微小網膜電位は数Hzから数100Hzまで102の周波数帯域にわたり多数の成分を含んでおり,その解析は極めて困難である。本研究では微小網膜電位の解析に多解像度変換を適用することにより,従来より詳細な成分分折ができることを示した。従来はフラッシュ網膜電位に関する知見から5つの波形成分を仮定したが,多解像度変換を用いた解析により,さらに多くの成分が含まれていることが示唆された。また連続性の解析により,微小網膜電位には網膜全野位から普遍的に発生する成分と,局在して発生する成分があり,また,特徴が分岐している可能性もあることが示唆された。異なった部位から発生する網膜電位をそれぞれ独立に解析した結果から,部位問で接続関係をもつ特徴パラメータが得られたことは,本研究のパラメータ推定法の妥当性を示すものと考えられ,従来の解析では得られなかった詳細な知見が示されたといえる。
網膜電位(ERG)は光刺激に対して網膜細胞層から発生する電位で,臨床眼科において診断及び治療評価に応用されるとともに,視覚機能の初期過程に関連するものとして,生体工学の立場からも関心がもたれている。現在,臨床眼科におけるERG計測にはフラッシュ刺激や市松模様のパターン刺激が使われているが,網膜全体を刺激しているため,局部的特性を知ることができない。このような状況に対し最近,ランダムドットパターン刺激により局部網膜電位を計測する方法が開発され(1),またERG計測結果から網膜機能の解析が試みられている(2)が,網膜部位における光刺激ERGの変換過程の全貌が解明されるところまでは至っていない。
このような問題に一石を投じるために,発光ダイオード(LED)を用いて網膜をランダム刺激し,発生するERGの周波数特性を相関解析によって求める微小網膜電位(BERG)計測法を開発した(3)。また,微小網膜電位はフラッシュERGについて知られているa波,b波,OP波等の原波形であることを示すとともに,最小2乗法を用いて網膜部位による各成分の特徴パラメータの分布状態を解析した(4)一(7)。しかし,a波,b波,OP波等は臨床眼科において便宜上に付けられた名称で,網膜各部位からどれほどの波形成分が発生しているのか,明かになっていない。しかし一方において,解析において多数の成分を考慮することは,それに応じて各成分を特徴付ける特徴パラメータの増加をきたし,これを計測された微小網膜電位の周波数特性から推定することが困難になる。
そこで,本研究では多解像度変換を用いて前処理を行った微小網膜電位に対して最小2乗法を適用することにより特徴パラメータを推定する方法を開発した。すなわち多解像度変換によって微小網膜電位の各成分波形は,不完全ながら対応するチャンネルに分離されるので,チャンネルごとに最小2乗法を適用すれば,特徴パラメータの推定が容易になると期待される。なお,この手法を適用するためには,従来周波数領域で行っていた相関スペクトル解析を等価的に時間領域で行う必要があり,解析手法を変更した(8)。
2.微小網膜電位の解析法
2.1刺激方法
まず図1に示すような実験システムを作成した。コンピュータからランダム刺激信号を発生し,V/F変換器(1V/1KHz)によってFM変調してLEDを発光させる。原刺激信号はDC成分6Vに振幅4Vのランダム信号を重畳しているので,2V以下にはならない。したがって,FM変調後の周波数は2KHz以上であり,LED光はこの周波数で点滅しているが,網膜はこのような高周波成分に応答できないので,網膜はランダム刺激信号に比例した輝度で刺激されていると考えられる。
また,商用周波数とともに,コンピュータやLED駆動回路等から発生するノイズを防止するために,今回新たに総銅板製のシールド暗室を作成し,刺激光を光ファイバを通してシールド内の被験者の眼前300mmの位置に導入する方法をとった。ファイバは8個所(網膜中心および周辺7個所)に配置し,その単独,あるいは任意の2個所を選択刺激できるようにした。なお,ファイバの発光端面と試験眼からの距離より,眼球内での光散乱を無視すれば,網膜照射面の直径は1deg程度になる。なお,周辺刺激時には刺激光とは異なる色の微弱なLEDで正面に注視点を作った。今回は刺激光は赤色,注視点は緑色である。
2.2計測方法
電極としてERG用コンタクトレンズ電極を試験眼(右眼)の角膜上に装着し,リットマン電極を不関電極として額中央,接地電極として試験眼側の耳たぶに装着する。試験眼には散瞳剤を点眼し,シールド暗室に入り約30分暗順応した後,計測を行った。
ERGをコンピュータに取込む際のサンプル周波数は1KHzとし,サンプル定理に基づいて,カットオッフ500Hzのバターワース型2次フィルタを2段前置して高周波成分を除去した。また眼球運動に伴う筋電位など低周波のドリフトを除去するために,0.2KHzを折点周波数とする不完全微分回路を挿入し,コンピュータに取込んだ後に逆周波数特性のディジタルフィルタで波形再生する方法を採った。計測装置の総合増幅率は0.2KHzにおいて最大100,000倍で調整できる。
なお,1回の計測時間は6秒間で,再現性のある疑似白色ノイズ刺激を用いて30回程度繰返し刺激し,得られた応答を加算平均した。この6000点の時系列データを基に以下に述べる相関分析を行った。
2.3相関スペクトル解析
相関スペクトル解析はノイズに埋没した信号を検出する周知の手法であるが,その原理を簡単に説明する。
刺激をu(t),電極電位をy(t),ノイズをn(t)とし,ERG発生機構の伝達特性をg(t)とすれば,電極電位は
と表される。両辺にu(σ一λ)を掛けて期待値をとれば
となる。ここで,疑似白色ノイズ刺激とノイズとは相関がなく,両者間の相互相関関数Φun(λ)は0に収束すると期待されるので,刺激,電極間の相互相関関数は刺激の自己相関関数とERG発生機構の特性との畳込み積分で近似できる。
従来は,(2)式をフーリエ変換することにより,ERG発生機構の伝達特性G(jω)を
として求めていた。すなわち,疑似白色ノイズ刺激とERGとの相互スペクトル密度関数Φxy(jω)の疑似白色ノイズ刺激の自己スペクトル密度関数Φxy(jω)に対する比より,ERGの周波数特性を求めた。また,5種類の成分を仮定し,各成分がそれぞれむだ時間をもつ2次振動波形
で近似できるものと仮定して,その合成波形により微小網膜電位のモデルを作成し,計測データとの周波数領域(Bode線図)における誤差を最小化することにより,各成分の特徴パラメータPn=(Κn,τn,ζn,Τn)を推定していた。
今回は多解像度変換を利用するために,上記の相関スペクトル解析を以下のように時間領域において行った。また,成分の数を増してより精密な解析を行った。
3. 多解像度変換による解析
時間領域における相関スペクトル解析の概念図を図2に示す。ここでCは相関演算,LSMは最小2乗法を表している。すなわち,長さ2Lのランダム信号uを用いて計測した微小網膜電位yとこのランダム信号の前半の長さLの部分U*との相互相関関数φ。yを求めれば,その結果は(3)式のΦ。yGω)と対応することが知られている。ここでランダム信号の前半部分の信号u*を用いるのは有限長データにおいて発生する折り返し現象を避けるための技法である。同様にランダム信号u,u*の自己相関関数より(3)式のΦuuGω)に対応するφuuが得られる。そこでこのφ。、を入力として,(4)式で与えられる微小網膜電位モデルGm(S,p,)の応答φ、mを求める。ここで計測した微小網膜電位G(s)とモデルGm(s,p。)が一致していれば,φuyとφ。mとは等しくなる。したがってφ。yとφ。mとの誤差を最小化することにより特徴パラメータp.が推定できる。
以上は2.3に述べた相関スペクトル法を時間領域で行ったことに相当するが,このφ。y,φ。mをそれぞれ多解像度変換した上で最小2乗法を適用することにより成分の特徴パラメータを推定した。多解像度変換は線形変換であるから,数学的には多解像度変換によって最小2乗推定の精度が向上することはない。しかし,微小網膜電位の成分はむだ時間をもつ振動波形で,いま簡単のために1つの成分波形を考えた場合,計測波形の第n波頭とモデルの第n±1波頭が一致するような時間シフトが淀こっても2乗誤差値が減少する。すなわち,最小2乗手定の評価関数は多峰性をもち,初期値の設定によって言った局所解へ収束する可能性がある。適正な初期値をπログラムによって自動的に設定することは極めて困鮪ディスプレイに表示された波形を観察して設定する必aがあるが,複数の成分が混合した微小網膜電位の場創はこれも容易ではない。多解像度変換を行えば,各成分が周波数の異なるチャンネルに不完全ながら分離され尾ので,観察による初期値設定が容易になることが期待冤きる。φ(t)の多解像度変換ψ(σ,d)は
で与えられるが,これは一種の並列フィルタと考えられ,dは各フィルタの解像度(分析周波数)を決めるパラメータである。代表的な重み関数wとしては矩形パルス状のHarr関数があり,パラメータdを1/2ずつ減少して解像度を上げながら分析する。重み関数にはdの異なる2つの重み関数間に直交性が要求され,Harr関数はこの直交関数の中で計算が容易な利点をもつものであるが,位相のシフトを伴うので,ここでは図3に示すようなHarr関数を対称化した関数を考え,関数幅dを1/2つつ縮小した複数のチャンネルに分析した。この重み関数は連続関数では直交性をもっているが,数値計算のためにサンプルした場合,負値をとるサンプル点数2nに対し,正値をとるサンプル点数が原点の1個だけ多くなるために直交性が失われる。そこで,原点における関数値を0とすることにより直交性を保存した。
4.成分分析
この多解像度変換後の各データΨiuyについて,それぞれ分析区間をとり,その区間でのΨiuy とΨiumとの2乗誤差の和
を評価関数とし,最小2乗法を用いて成分分析を行った。
最小2乗法に従えば,(4)式のモデルにおける特徴パラメータの修正量dp,は,評価関数Jのp,による偏微分
で与えられる。これに(6)式を代入すれば
となるが, は
と表される。すなわち,モデルの伝達関数Gm(s,pn)をpnで偏微分したフィルタをプログラムし,これに多解像度変換したランダム信号を入力すれば∂Ψium / ∂Pnが得られる。また, ∂Ψium / ∂Pnと誤差(Ψiuy - Ψium)との相関演算より特徴パラメータの修正量が求められる。
なお,評価関数計算における誤差2乗積分の上下限の設定については,計測データ(あるいはランダム信号との相互相関関数φuy)のどの区間を有意とみなすかということに帰着する。フラッシュERGにおける成分波形は刺激後100ms以内に含まれているので,これを目安に100msまでの区間とした。また,波形成分の数の設定についても現時点では必然的な根拠はなく,従来はフラッシュERGの成分数より5個としていた。今回の多解像度変換を用いた成分分析では従来より詳細な解析ができるため,5個の成分数では計測データとモデルとの充分な一致が見られない場合が多く,微小網膜電位にはさらに多数の成分が含まれていることが示唆された。本来,網灰刺激部位に含まれる多数の細胞がそれぞれが独自の活動をしており,その細胞数の成分が必要ということにもなるが,現実にはそれが何個に種別できるかということであり,成分の増加によって最小2乗推定に要する演算時間が急激に増加するため,今回の解析では9個に設定した。これらの成分の内には後述のように網膜部位に局限して発生しているものと,網膜の比較的広い領域から共通して発生しているものとがあり,後者の普遍的成分に限れば,やはり5成分程度であるように思われる。
5.特徴パラメータの位置連続性の判定
網膜水平軸に沿った何点かの部位で計測した微小網膜電位から推定した特徴パラメータには後述のように局在的なものと,比較的広い領域に一定の傾向で変化しているものがあるようにみえる。しかし,隣接した部位間でどの成分同士が連続しているか,データを観察するだけで判断するのは容易ではない。そこで,隣接部位間について以下のような評価関数を考え,連続性の判定の目安にした。
特徴pn(たとえばむだ時間)について,部位kから発生したERGのi番目の成分p。kiと隣接部位k+1から発生したERGのj番目の成分p。(k+1)jとの2乗誤差の和
を評価関数とし,部位kとk+1との間の隣接行列を作る。ここでw,は特徴パラメータについての荷重である。iを固定したときの最小値を与えるj(部位kから見て最も近い成分番号)を探し,隣接行列にマークする。同様にjを固定したときの最小値を与えるi(部位k+1から見て最も近い成分番号)をマークする。これらのマークが一致したときには,それは最も連続性の高い組合わせと考える。また,マークが全く付かないものの間には連続性がないと判断する。隣接行列から連続性の高い組合わせを除き,行,列のどちらか一方にマークが付いたものが残っていれば,その成分について同様の評価を行い,つぎに連続性のある組合わせを判定する。この判定結果は荷重の値によって左右されるが,特徴パラメータの間に連続性をもつものができるだけ多くなるように調整することで,連続性の判定を行った。また,部位iのERGに近接した特徴をもつ2成分A,A'が含まれているとき,その一方Aが部位jと接続関係をもつと判断されるとA'が孤立する可能性がある。この場合,特徴が部位iにおいてA,A'に分岐していると考えるべきかも知れないが,現時点では特徴パラメータはそれぞれ1本の鎖で繋がれると仮定して解析を行った。
6.解析結果
解析例として網膜中心部刺激に対する微小網膜電位(実線)とモデル応答(破線)の多解像度変換パターンの例(網膜中心部刺激時)を図4に示す。
今回の解析では多解像変換のチャンネル数は8とし,最低次のチャンネルの荷重関数幅4dを1024点として,次数の増加に伴い荷重関数幅を1/2ずつ縮小し,解像度を上げた。したがって,多解像度変換パターンは次数が上がる程,原点付近の状況を拡大して表示している。最高次のチャンネルには最も周波数の高い成分波形のみが抽出されていると考えられ,これをモデルで近似することにより,この成分の特徴パラメータが推定される。より低い次数のチャンネルには複数の成分が含まれているが,周波数の高い成分がすでに推定されているので,他の成分の特徴パラメータも比較的容易に推定される。計算上は各チャンネルについて同時に評価関数の最小化を行っているが,定性的には以上のように高次のチャンネルから推定が進行していくように思われる。したがって,低次数のチャンネルに見られる細かいリップル様の変動まで近似することができる。このような精度の高い近似は従来の手法では得られなかったもので,多解像度変換の効果を示すものである。
一方,表1は最小2乗法で推定された特徴パラメータの接続関係を5.に述べた手法により解析した結果である。ここで,fは中心周波数,ζは減衰係数,τはむだ時間である。また反応強度は,減衰係数が小さい成分では直流ゲイン係数Kのみでは表現できないので,1V*secの仮想インパルス刺激に対する応答の最大振幅を,
により計算し,これを反応強度gmとした。刺激部位は中心部(CENTER),こめかみ側10,20,45deg(T10,T20,T45)および鼻側10,20,30,45deg(N10,N20,N30,N45)の8部位である。各部位において網膜電位は一旦9個の成分に分離されるが,5.に述べた隣接行列を適用すると12成分に分けられた。隣接行列法のもつ前記の問題点も考慮して連続性の高い(全網膜に共通する)6成分を示すと表1のようになった。
これらの他に中心周波数20～25HzでT45を除く網膜部位で発生している成分がみられたが,その他は以下のように局在性か,特徴の分岐による断片であった。すなわち,中心周波数30～40Hz発生部位T20,T45,中心周波数60～75Hz発生部位CENTERからT45までのこめかみ側,中心周波数220Hz発生部位CENTER,N10が見られた。
表1に示した6成分は網膜全野に共通する成分で,その周波数帯域より1はc波,III,IV波はb波,Vはa波,VI波はOP波と対応すると考えられる。これらの結果は先の報告(4)の結果にも矛盾しないが,II波はb波としてはやや周波数帯域が低いように思われる。フラッシュERGの対応関係は明かでないが,今回の解析による新しい知見ともいえる。
7.むすび
微小網膜電位は数Hzから数100Hzまで102の周波数帯域にわたり多数の成分を含んでおり,その解析は極めて困難である。本研究では微小網膜電位の解析に多解像度変換を適用することにより,従来より詳細な成分分折ができることを示した。従来はフラッシュ網膜電位に関する知見から5つの波形成分を仮定したが,多解像度変換を用いた解析により,さらに多くの成分が含まれていることが示唆された。また連続性の解析により,微小網膜電位には網膜全野位から普遍的に発生する成分と,局在して発生する成分があり,また,特徴が分岐している可能性もあることが示唆された。異なった部位から発生する網膜電位をそれぞれ独立に解析した結果から,部位問で接続関係をもつ特徴パラメータが得られたことは,本研究のパラメータ推定法の妥当性を示すものと考えられ,従来の解析では得られなかった詳細な知見が示されたといえる。