(57)【要約】  本発明は、高周波の深部脳刺激療法（ＤＢＳ）中の神経細胞活性の信号をフィルタリングして、観測された信号中の刺激アーチファクトを除去するための方法に関し、該方法は、位相空間中で、観測された信号軌跡を近似する工程であって、該観測信号は、刺激信号によって誘起された刺激アーチファクトの和として見なされ、該刺激信号は、安定リミット・サイクルによる自励振動システムを含む常微分方程式の解であると仮定される、工程と、該観測信号およびその微分係数をセグメントにスライスする工程であって、各セグメントは、刺激期間に対応する、工程と、学習セットに、Ｎ個の選択された刺激期間を集める工程と、該自励振動システムの該リミット・サイクルを見積もる工程と、該観測信号の各アーチファクトを該見積もられたリミット・サイクルに同期させる工程と、該同期によって、位相空間中で各アーチファクトから該見積もられたリミット・サイクルを減じる工程と、フィルタリングされた信号全体を得るために、すべてのセグメントを集める工程と、最後に、時間領域中でその結果を表現する工程とを含む。【選択図】図１
2010524612.tif 000002
【特許請求の範囲】
【請求項１】
高周波の深部脳刺激療法（ＤＢＳ）中の神経細胞活性の信号をフィルタリングして、観測される信号中の刺激アーチファクトを除去するための方法であって、
ｉ）位相空間中で、観測された信号軌跡を近似する工程であって、観測される信号は、刺激信号によって誘起された刺激アーチファクトの和として見なされ、前記刺激信号は、安定リミット・サイクルによる自励振動システムを含む常微分方程式の解であると仮定される、工程と、
ｉｉ）前記観測信号およびその微分係数をセグメントにスライスする工程であって、各セグメントは、刺激期間に対応する、工程と、
ｉｉｉ）学習セットにＮ個の選択された刺激期間を集める工程と、
ｉｖ）前記自励振動システムの前記リミット・サイクルを見積もる工程と、
ｖ）前記観測信号の各アーチファクトを前記見積もられたリミット・サイクルに同期させる工程と、
ｖｉ）前記同期によって、位相空間中で各アーチファクトから前記見積もられたリミット・サイクルを減じる工程と、
ｖｉｉ）フィルタリングされた信号全体を得るために、すべてのセグメントを集める工程と、
ｖｉｉｉ）時間領域中でその結果を表現する工程と、
を含む方法。
【請求項２】
位相空間中の前記信号軌跡は、ガウス畳み込みを使用して、平滑された信号と、信号の微分係数の近似を計算することによって、近似される、請求項１に記載の方法。
【請求項３】
各タイム・セグメントの中心が、前記観測信号中の最高の信号ピークに近い、請求項１または２のいずれかに記載の方法。
【請求項４】
前記セグメントは、ランダムに選択される、請求項１乃至３のいずれかに記載の方法。
【請求項５】
前記自励振動システムの前記リミット・サイクルは、近傍で確率密度の最大値を与える、学習セットからの要素を使用して見積もられる、請求項１乃至４のいずれかに記載の方法。
【請求項６】
安定リミット・サイクルによる前記自励振動システムは、４次元の自励振動システムである、請求項１乃至５のいずれかに記載の方法。
【請求項７】
前記ガウス畳み込みは、サンプリング・レートが４８ｋＨｚの場合、３〜４ｋＨｚの範囲に含まれる特性周波数について計算される、請求項２に記載の方法。
【請求項８】
前記特性周波数は、３．５２ｋＨｚに等しい、請求項７に記載の方法。
【請求項９】
前記学習セットは、Ｎ個の選択された期間から集められ、少なくとも１０００個のシミュレーションのアーチファクトを含む、請求項１乃至８のいずれかに記載の方法。
【請求項１０】
実体のある媒体中に組み込まれ、その命令がマシーンにオペレーションを実施させるコンピュータ・プログラムであって、
ｉ）位相空間中で、観測された信号軌跡を近似する工程であって、観測される信号は、刺激信号によって誘起された刺激アーチファクトの和として見なされ、前記刺激信号は、安定リミット・サイクルによる自励振動システムを含む常微分方程式の解であると仮定される、工程と、
ｉｉ）前記観測信号およびその微分係数をセグメントにスライスする工程であって、各セグメントは、刺激期間に対応する、工程と、
ｉｉｉ）学習セットにＮ個の選択された刺激期間を集める工程と、
ｉｖ）前記自励振動システムの前記リミット・サイクルを見積もる工程と、
ｖ）前記観測信号の各アーチファクトを前記見積もられたリミット・サイクルに同期させる工程と、
ｖｉ）前記同期によって、位相空間中で各アーチファクトから前記見積もられたリミット・サイクルを減じる工程と、
ｖｉｉ）フィルタリングされた信号全体を得るために、すべてのセグメントを集める工程と、
ｖｉｉｉ）時間領域中でその結果を表現する工程と、
を含むコンピュータ・プログラム。
【請求項１１】
位相空間中の前記信号軌跡は、ガウス畳み込みを使用して、平滑された信号と、信号の微分係数の近似を計算することによって、近似される、請求項１０に記載のコンピュータ・プログラム。
【請求項１２】
各タイム・セグメントの中心が、前記観測信号中の最高の信号ピークに近い、請求項１０または１１のいずれかに記載のコンピュータ・プログラム。
【請求項１３】
前記セグメントは、ランダムに選択される、請求項１０乃至１２のいずれかに記載のコンピュータ・プログラム。
【請求項１４】
前記自励振動システムの前記リミット・サイクルは、近傍で確率密度の最大値を与える、学習セットからの要素を使用して見積もられる、請求項１０乃至１３のいずれかに記載のコンピュータ・プログラム。
【請求項１５】
安定リミット・サイクルによる前記自励振動システムは、４次元の自励振動システムである、請求項１０乃至１４のいずれかに記載のコンピュータ・プログラム。
【請求項１６】
前記ガウス畳み込みは、サンプリング・レートが４８ｋＨｚの場合、３〜４ｋＨｚの範囲に含まれる特性周波数について計算される、請求項１１に記載のコンピュータ・プログラム。
【請求項１７】
前記特性周波数は、３．５２ｋＨｚに等しい、請求項１６に記載のコンピュータ・プログラム。
【請求項１８】
前記学習セットは、Ｎ個の選択された期間から集められ、少なくとも１０００個のシミュレーションのアーチファクトを含む、請求項１０乃至１７のいずれかに記載のコンピュータ・プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、脳信号の処理に関し、より詳しくは複数の信号を表す脳信号を処理するための方法およびシステムに関する。
【背景技術】
【０００２】
高周波（１００〜３００Ｈｚ）の深部脳刺激療法（ＤＢＳ：Deep Brain Stimulation）は、たとえばパーキンソン病の衰弱性症状など、様々な障害を引き起こす神経症状を治療するための外科的処置である。
【０００３】
ＤＢＳの間、電気的な刺激インパルスが、患者の脳中のターゲット領域に送られる。
【０００４】
ＤＢＳから生じる問題は、患者の神経細胞活性（神経細胞活動電位ＡＫＡスパイクによって伝えられる）を刺激期間中に測定することができないということである。大きな刺激アーチファクトによって、スパイク検出が妨げられ、刺激中の単細胞活性を解析する上での問題を引き起こしている。
【０００５】
神経細胞活性の特異的信号および刺激アーチファクトが、記録上に互いに積み重なっている。
【０００６】
刺激アーチファクトは、波形が通常であるが、サンプリング・エラーおよび刺激生成の不規則性によって、同一でない。アーチファクト信号比（ＡＳＲ：Artefact to Signal Ratio）は、刺激前に観測された神経細胞活性のスパイクの平均振幅に対する刺激アーチファクトの振幅の平均値の比であり、５〜２０の範囲内に含まれる。
【０００７】
電気生理学の分野では繰り返される、この問題を回避するために、いくつかの方法が、このアーチファクトを抑制するために開発されてきた。それらの中で、記録中の刺激アーチファクトを減じるための最近の方法をHashimoto他が述べている（Hashimoto T、Elder C. M.、Okun M. S.、Patrick S. K.、Vitek J. L.、Stimulation of the subthalamic nucleus changes the firing pattern of pallidal neurons、J. Neurosci.、 2003、 23:1916〜1923）。
【０００８】
この方法は、以下の一般的なプロトコルによって、刺激周囲セグメント（peri-stimulus segments）のセットの平均化を介してアーチファクト・テンプレートを作成することに基づく。アーチファクトのばらつきを考慮に入れるために、いくつかのテンプレートが生成された。例として示した信号は、ＡＳＲがほぼ６に等しく、それは、アーチファクト振幅が、スパイクより６倍大きいことを意味する。
【０００９】
テンプレート減算後、信号は、スパイクより１．５〜２倍大きいアーチファクト残渣をなお含んでいる。時間領域中の残渣のサイズが、約０．８ｓであった。これらのウィンドウは、さらに解析する前に、信号から除去された。したがって、この部分（記録時間の約１１％）は、スパイクの分類またはさらなるどのような解析にも、もう利用できない。そのような切り捨てによって、解析される記録から得られるスパイク列が、著しく損なわれることになり得る。
【００１０】
本明細書に組み込まれて、その一部分を構成する添付図面は、本発明の様々な実施形態を示し、そして記述とともに、本発明の原理を説明するのに役立つ。
【図面の簡単な説明】
【００１１】
【図１】本発明による、高周波の深部脳刺激療法（ＤＢＳ）中の神経細胞活性の信号をフィルタリングして、観測される信号中の刺激アーチファクトを除去するためのシステムの概略図である。
【図２】本発明による、高周波の深部脳刺激療法（ＤＢＳ）中の神経細胞活性の信号をフィルタリングするための方法のブロック図である。
【図３】ＤＢＳの前およびその間の神経細胞活性の第１の信号を示す図である。
【図４】フィルタリング後の図３のＤＢＳの前およびその間の神経細胞活性の第１の信号を示す図である。
【図５】ＤＢＳの前およびその間の神経細胞活性の第２の信号を示す図である。
【図６】フィルタリング後の図５のＤＢＳの前およびその間の神経細胞活性の第２の信号を示す図である。
【図７】フィルタリング前の刺激期間のウィンドウにスライスされた第１の信号を示す図である。
【図８】フィルタリング後の刺激期間のウィンドウにスライスされた、図７の第１の信号を示す図である。
【図９】フィルタリング前の刺激期間のウィンドウにスライスされた第２の信号を示す図である。
【図１０】フィルタリング後の刺激期間のウィンドウにスライスされた、図７の第２の信号を示す図である。
【図１１】安定リミット・サイクルの近傍中の局所座標を示す図である。
【図１２】第１の記録について、座標が、一次、二次および三次の微分係数に等しい、位相空間射影中のアーチファクト軌跡を示す図である。
【図１３】第２の記録について、座標が、一次、二次および三次の微分係数に等しい、位相空間射影中のアーチファクト軌跡を示す図である。
【図１４】信号フィルタリングの前およびその後のアーチファクト・スパイク比（ＡＳＲ：Artefact Spike Ratio）を示す表である。
【図１５】ピーク領域中のアーチファクト残渣の標準偏差（ＳＴＤ：standard deviation）を示す図である。
【図１６】加法ノイズを伴うモデル（破線）および非線形の振動モデル（実線）による信号の表現に基づく、位相に依存したアーチファクトのＳＴＤを示す図である。
【図１７】アーチファクトのピーク領域中のアーチファクト残渣のＳＴＤを示す表である。
【図１８】フィルタリング後のスパイク検出のエラーのパーセンテージを示す図である。
【図１９】検出されたスパイクについて、第１の２つの主な成分の分散を示す図である。
【図２０】テストのために使用される合成信号中のＳＮＲ（Spike-Noise Ratio：スパイク・ノイズ比）を示す図である。
【図２１】合成データおよび実データに対する時間領域アルゴリズムと位相空間アルゴリズムを比較した、フィルタリング前のＡＳＲに依存した残渣−スパイク比を示し、８つの実際の記録（マーカー）および合成信号（破線および実線）についての結果が良好にフィットしていることを示す図である。
【発明の概要】
【課題を解決するための手段】
【００１２】
ここで、添付の図面に示す、本発明の例示の実施形態を詳しく述べる。可能な場合、同じ参照番号が、同じ、または同様の部分に言及するために、図面すべてにわたって使用される。
【００１３】
図１を参照すると、実験の状況の実施例が示されており、記録および刺激が、同時に実施され、装置は、ガラスが被覆されたプラチナ・イリジウムの微小電極または同様のもの１００を含み、それは、単一（または複数）単位の記録のために、脳１０５の大脳半球の構造中に挿入される。たとえば、微小電極は、内部淡蒼球、他の任意の細胞核、皮質または白質中に挿入される。微小電極１００を介して、刺激電極１１０（これは、他の微小電極または常在的なマクロ電極とすることができる）が、脳構造（たとえば、視床下核）内に埋め込まれる。電極１１０のコンタクトは、一般に刺激装置と言われる、プログラム可能なパルス発生器１１５に接続される。
【００１４】
増幅器１２５、たとえばフィルタ１３５に付属したＡ／Ｄ変換ユニット１３０を含む電気生理学的収集システム１２２を使用して、信号１２０が、電極１００によって記録された。これは、前記フィルタリングされた信号１２０を生成する。
【００１５】
さらに、装置は、前記フィルタリングされた信号１２０を解析するために、解析ユニット１３８を含む。解析ユニット１３８は、コンピュータ・プログラム、すなわちアルゴリズム（通常、コンピュータ中にインストールされて実行される、あるいは収集システム１２２中にモジュールとして含まれる）であり、刺激アーチファクトが信号１２０から除去された信号１４０を生成する。
【００１６】
ここで、図２を参照すると、高周波の深部脳刺激療法（ＨＦ−ＤＢＳ：high frequency deep brain stimulation）中の神経細胞活性の信号をフィルタリングして、観測された信号１２０から刺激アーチファクトを除去するための方法が、観測信号１２０を記録する工程２００を含む。
【００１７】
工程２１０で、位相空間中で観測信号の軌跡を近似するための畳み込み、特にガウス畳み込みを含む方法の１つを使用して、観測信号の一次、二次および三次の微分係数が計算される。微分係数は、位相空間中で信号軌跡を表現するのに役立つ。ガウス畳み込みが、３〜４ｋｈｚの範囲に含まれる特性周波数について、優先的に３．５２ｋｈｚに等しい特性周波数について計算される。観測信号は、電気生理学的な信号、刺激信号（刺激の信号によって誘起されるアーチファクト）およびバックグランド・ノイズの和であると考えられる。
【００１８】
刺激信号は、安定リミット・サイクルによる自励振動システムを記述する常微分方程式の解と仮定される。
【００１９】
工程２１５で、観測信号およびその微分係数は、セグメントにスライスされ、各セグメントは、刺激期間に対応する。さらに、各タイム・セグメントの中心は、観測信号中の最高の信号ピークに近い。
【００２０】
さらに、工程２２０で、Ｎ個の選択された刺激期間が、学習セットに集められ、前記学習セットは、少なくとも１０００個のシミュレーションのアーチファクトを含み、工程２２５で、自励振動システムのリミット・サイクルが見積もられる。自励振動システムのリミット・サイクルは、近傍で確率密度の最大値を与える、学習セットからの要素を使用して、見積もられる。学習セットの要素を使用する、他のどのような一貫性のある統計学も適用することができる。好ましくは、安定リミット・サイクルによる自励振動システムは、４次の自励振動システムとする。
【００２１】
次いで、工程２３０で、観測信号の各アーチファクトは、見積もられたリミット・サイクルに同期され、工程２３５で、見積もられたリミット・サイクルは、同期によって位相空間中で各アーチファクトから減じられる。
【００２２】
最終的に、工程２４０で、フィルタリングされた信号全体を得るために、すべてのセグメントが集められ、工程２４５で、その結果が、時間空間中で表現される（信号１４０）。
【００２３】
方法の原理
観測信号x(t)、t=1, 2, …は、刺激アーチファクトxSt(t)および神経細胞活性の信号xNr(t)の和、すなわちx(t)=xSt(t)+xNr(t)と考える。摂動による常微分方程式の解と仮定される刺激信号は、
【数１】
2010524612.tif 000003
ただし、ｎは、方程式の次数であり、F()は、摂動関数であり、次の式は、座標を
【数２】
2010524612.tif 000004
とする位相空間中で、安定リミット・サイクルによる自励振動システムx0(t)=(x01(t),…, x0n(t))'、0<t≦Tを記述する。
【数３】
2010524612.tif 000005
ここで、Ｔは、刺激期間である安定振動の期間である。摂動関数F()は、微小値によって束縛され、平均値がゼロであり、Δt>Δtcorr、Δtcorr<<Tの場合、安定振動相関時間：Ｂ（Ｆ( ,ｔ)は，F( ,t+Δt)≒0に比較して小さい、ランダム・プロセスである。
【００２４】
安定振動（式３）の場合、信号の軌跡は、近傍で見られるときはいつも、リミット・サイクルになる傾向がある。それは、初期点がリミット・サイクル上に位置する場合、周期解を与える。式（２）の摂動関数Ｆ()は、リミット軌跡から外に信号の軌跡を移動させる傾向がある。しかし、摂動が十分小さい場合、その軌跡は、リミット・サイクルX0(t)、0<t≦Tの近傍に留まる、すなわち式（２）の解は、互いに相似であるが、それらは、決して一致しない。
【００２５】
リミット・サイクルの近傍で局所座標を、（Gudzenko, L. I.、 Statistical method for self-oscillating system characteristics detection、Izvestiia Vuzov Radiophysics、vol. 5、 pp. 573〜587、1962）に従い、導入しよう。リミット・サイクル上の任意点Ｐ0を開始点として決定しよう。
【００２６】
図１１を参照すると、リミット軌跡上のいかなる任意点Ｐの位置も、その位相θによって記述することができ、それは、P0（θ=0）として定義される開始点Ｐ0からのリミット・サイクルに沿った時間運動である。式（２、３）の関数ｆ()が、関数の独立変数すべてに関し２回連続的に微分可能であり、したがって必要な滑らかさが与えられると仮定しよう。位相がθの点Ｐで、リミット・サイクルに垂直である超平面（そして、そのような超平面１つだけ）を構築することが可能である。リミット・サイクルの近傍に式（２）の任意軌跡を考えよう。
【００２７】
位相θの超平面に対するその交差点をＭ（θ）で表し、位相がゼロの点M（0）、θ=0を解析する軌跡の初期点として設定する。いかなる軌跡も、変数ｎ（θ）およびｔ（θ）（図１１）によって記述することができる、ただし、ｎ（θ）は、Ｍ（θ）およびリミット・サイクル上のその正射影P（θ）によって定義される垂直偏差のベクトルである。第２の変数ｔ（θ）は、初期点Ｍ（０）から解析点Ｍ（θ）への軌跡に沿った時間運動である。したがって、リミット軌跡は、n（θ）≡O、t（θ）≡θによって定義される、ただし、Oは、すべての成分が0に等しいベクトルである。
【００２８】
γ（θ）を位相偏差γ（θ）=t（θ）-θとしよう。式（２）は、偏差ｎ（θ）およびγ（θ）によって、線形近似で以下のように書き換えることができる（Gudzenko 1962）。
【数４】
2010524612.tif 000006
ここで、Ｎ（θ）およびΘ（θ）は、パラメータの関数である。その結果、位相空間中の信号軌跡は、リミット・サイクルの周期的成分と偏差の関数の和x(t(θ)=x0(θ)+n(θ)として、線形近似で表される、ただし、ｎ（θ）およびt（θ）=γ（θ）-θは、式（４）に従う。なお、本発明の方法によれば、摂動による非線形振動のモデル（式２）は、振幅および位相の両方のゆがみを説明する（式４）。
【００２９】
以下のアーチファクト・フィルタリングのために、信号の周期的成分を表すリミット・サイクルを見積もるべきである。リミット・サイクルは、x0(θ)、0<θ≦Tとして記述され、位相は、０からＴまでに及ぶ。位相θが０からＴまでに及ぶサイクルとして言及される、任意の信号軌跡xi(t(θ))、0<θ≦Tのセグメントを考えよう。一般集団X={xi(t(θ))、0<θ≦T}が、サイクルによって形成され、刺激アーチファクトの認識に対応する。リミット軌跡x0(θ)、0<θ≦Tが、理想的な「乱されていない刺激アーチファクト」に対応する。任意サイクルは、以下の式のように表される。
xi(t(θ))＝x0(θ)+ni(θ)、0<θ≦T                            式（５）
ただし、ni(θ)は、式（４）によって決定される。平均軌跡は、平均軌跡数が無限に増加する場合、線形近似でリミット・サイクルx0(θ)に収束する（Gudzenko, L. I.、Statistical method for self-oscillating system characteristics detection、Izvestiia Vuzov Radiophysics、vol. 5、pp. 573〜587、1962）。
【００３０】
それによって、位相空間中で、観測されたサイクルxi(t(θ))、0<θ≦Tの平均値を計算してリミット・サイクルx0(θ)を見積もることが可能になる。
【数５】
2010524612.tif 000007
【００３１】
ni(θ)およびγi（θ）が、無相関ノイズF()の場合、任意のθについて漸近的なガウス分布によって特徴づけられること（Gudzenko, L. I.、Statistical method for self-oscillating system characteristics detection、Izvestiia Vuzov Radiophysics、vol. 5、 pp. 573〜587、 1962）を考慮に入れる、他の統計学も使用することができる。
【００３２】
位相空間中での刺激アーチファクトの平均値
【数６】
2010524612.tif 000008
の見積もりの後、それは、信号から減じられる。これに役立つように、局所座標
【数７】
2010524612.tif 000009
が、各サイクルについて計算される。次いで、垂直偏差のベクトルの関数
【数８】
2010524612.tif 000010
が、位相空間中でのフィルタリングの後の信号を表す。垂直偏差の第１の座標が、時間領域中のアーチファクトの残渣xfilt(t)=n1(t)+xNr(t)、t=1, 2,…を表す。スプライン内挿法が、時間の一定間隔の分割のために使用される。
【００３３】
フィルタのアルゴリズムおよびその実現
実現されたアルゴリズムは、以下の工程を含む。
１．位相空間中で信号軌跡を近似する。この目的のために、ガウス畳み込みを使用して、平滑された信号を計算し、信号の微分係数の近似を計算する。したがって、微分係数x(k)(t)、k=0, 1,…, n-1は、以下の式によって近似される。
【数９】
2010524612.tif 000011
２．信号およびその微分係数を刺激期間にスライスする。セグメントのそれぞれは、別のアーチファクトを含む。各時間セグメントの中心は、最高の信号ピークに近い。
３．Ｎ個のランダムに選択された刺激期間を学習セットXNに集める。
４．リミット・サイクルを見積もる。対称的な単一モードの分布密度の場合、数学的期待値x0は、それらの近傍で確率の最大値を与える。
P(d(x、x0)<R)→最大                                      式（８）
式（８）によれば、近傍で確率密度の最大値を与える、学習セットからの要素は、リミット・サイクルを見積もるために使用された。
【数１０】
2010524612.tif 000012
ただし、｜.｜は、セットの要素数である。
サイクルx(θ)、0<θ≦Tは、
【数１１】
2010524612.tif 000013
に属し、
【数１２】
2010524612.tif 000014
中のユークリッド距離が、d(x、y)として用いられた。ベクトルxi(θ)、0<θ≦Tは、それ自体得られない。したがって、2つのサイクルxi(θ)とxj(θ)、0<θ≦Tの間の距離は、信号のスプライン内挿法を使用して、次の式に近似した。
【数１３】
2010524612.tif 000015
ここで、τmaxは、位相偏差の最大値である。局所座標の定義によって、見積もられたリミット・サイクルの時間は、位相と見なした。
５．次の式を計算して、各アーチファクトをリミット・サイクルに同期させる。
【数１４】
2010524612.tif 000016
６．同期によって位相空間中で各アーチファクトから平均サイクルを減じる。
７．フィルタリングされた信号全体を得るために、すべてのセグメントを集める。
８．時間領域中でその結果を表現する。
【００３４】
アルゴリズムのパラメータが次のように選択された
モデルの次数
位相空間中でリミット・サイクル上の自己との交差点を避けるために、図４および５を参照して、4次n=4のモデルを使用する。
【００３５】
ガウス畳み込みの滑らかさのパラメータ
式（７）のすべての微分係数近似について、同じパラメータαを適用することを決定した。この戦略によって、ＯＤＥの近似が、2次の項までオリジナルの式に同一であることを保証することができる。適切なαの判定基準として伝達関数を使用する。演算子Ｄkについて、伝達関数は、次の式の形をとる。
【数１５】
2010524612.tif 000017
その最大値の位置は、
【数１６】
2010524612.tif 000018
である。最大値の平均から計算される特性周波数
【数１７】
2010524612.tif 000019
を使用する。畳み込みの計算を簡単化するために、αのためのサンプリングきざみ（tick）の数を整数に限定する。サンプリング・レートが４８ｋＨｚであることを考慮して、これは、最適のαがα=6.25×10-5 sec=3きざみであり、fc=3.52 kHzである。fcのこの値は、スパイクおよびアーチファクトの観測されたスペクトル帯域と最良に一致する。解析される信号では、ほとんどのエネルギーが、０．５〜８ｋＨｚの帯域中に集中する。
【００３６】
学習データ・セットのサイズ
学習セットに対して、N=1000の刺激アーチファクト（7秒間の観測）を取った。
【００３７】
最大位相偏差
τmaxを見積もるために、位相偏差が、アーチファクトのピークと交差して最大値に達すると仮定し、アーチファクトの最大ピークの時間のサンプリング分布を使用する。２個の続いて起きるピークの時期の間の差は、次の式に等しい。
ΔtP(k)=T+γk(θP)-γk-1(θP)                              式（１３）
ただし、θPは、ピークの位相であり、γk(θP)は、k番目のアーチファクトについての位相偏差である。摂動関数の相関時間が短いことを考慮して、分散がσ2P=σ2(θP)である独立したランダム変数として、γk(θP)を考えた。この場合、σ2(ΔtP) =2σ2Pが成立する。ピーク時間のセットから直接σ(ΔtP)を見積もった。次いで、最大位相偏差は、τmax=3σPとして固定された。漸近的な正規分布のγk(θP)の場合、それは、信頼度レベル０．９９７に相当する。たとえば、信号＃Ｓ１５２９について、σP=2.5927×10-5 secが得られた。最大位相偏差は、τmax=3σP=7.7782×10-5 secとして固定された。次いで、τmaxは、サンプリングきざみをもっとも近い整数（サンプリング・レート４８ｋHzの場合、４きざみ）に丸めて、調節された。
【００３８】
基本半径R
半径のパラメータは、アーチファクトの観測されたサイクルの分散によって計算される。すなわち、サンプリング・フラクタイル（fractile）が使用される。上記に述べたように、サイクルx(θ)、0<θ≦Tは、ユークリッド空間
【数１８】
2010524612.tif 000020
中に正規に分布するとして考えることができる。２つの独立したｘとｙの差は、またガウス分布である。したがって、ランダム変数
【数１９】
2010524612.tif 000021
およびz=x-yは、完全に同じように分布し、さらに、ユークリッド距離
【数２０】
2010524612.tif 000022
と
【数２１】
2010524612.tif 000023
は同様である。サイクルの期待値が分からないので、d2(x、y)/2の分布密度関数を見積もる。これに役立つように、学習セットからサイクルのランダムなペアを生成し、ヒストグラムを構築した。最大値の位置をフラクタイルq0.34〜q0.35にほぼ相当する基本半径Rとして取った。フラクタイルまたは変位値（quantile）は、基本的に、ランダム変数の累積分布関数から一定間隔で取られた点である。順序付けたデータを基本的に等しいｑ個のサイズのデータ・サブセットに分割するのは、ｑ変位値のために動機づけるものである。フラクタイルまたは変位値は、連続したサブセット間の境界を作るデータ値である。より正式には、母集団パラメータXのk番目の「ｑ」変位値またはフラクタイルは、次の式を満たすような「x」として定義することができる。
P(X≦x)≧p、およびP(X≧x)≧1-p、ただし
【数２２】
2010524612.tif 000024
フラクタイルを計算するために、ヒストグラムをχ2d分布によって近似する。有効自由度dが、上記に述べたフラクタルに繋がる７〜８として見積もられた。したがって、サイクルの３４．１〜３５．４％が、リミット・サイクルのRのボール内に含まれる必要がある。実験では、半径Rのボール内に、サイクルの３５〜４０％が、直接計数して得られた。
【００３９】
アルゴリズムのテスト
手術室からのいくつかの記録が、アルゴリズムをテストするために使用された。図３〜６に、オリジナル信号およびフィルタリングされた信号を示す。図７〜１０に、刺激期間のウィンドウのセットにスライスされた、フィルタリングの前およびその後の信号を示す。これらの図は、オリジナルのアーチファクト、およびすべての処理後に残された残渣の形状を明らかにする。また、神経細胞のスパイクを見ることができる。
【００４０】
図１４に、位相空間および時間領域中の信号フィルタリングの前およびその後のアーチファクト・スパイク比（ＡＳＲ：Artefact Spike Ratio）を示す。記録に関するノイズ・スパイク比（ＮＳＲ：Noise Spike Ratio）は、フィルタリング後のＡＳＲと比較するために、さらに計算される。ＮＳＲは、理想的なフィルタリングの場合、ＡＳＲによって得ることができる値を表す。
【００４１】
図１２および１３に、しかるべきアーチファクトのセットの相図を３Ｄで表す。図１２および１３に、図３からの第１の記録のスパイク軌跡（より小さい軌道）を見ることができる。
【００４２】
従来技術による加法ノイズを伴う標準モデルと位相モデル（式２）による信号記述の性質を比較するために、残渣のしかるべき標準偏差を計算した。モデル（式１）の加法ノイズζk(t)の標準偏差σζ(t)、0<t≦T、および垂直偏差のベクトルの第１の座標（式５）n1(t)の標準偏差σn1(θ)、0<θ≦Tを見積もった。値n1(t)は、モデル（式２）について、時間領域中のアーチファクトの残渣を表す。位相モデルが、特にピークの近くで、さらに一層小さい分散を生成することを見ることができる。図１５に、σカーブの正確な特性を示す。
【００４３】
本発明によれば、非線形振動モデルによって、加法ノイズを伴う標準モデルに比べて、刺激アーチファクトの残渣の標準偏差を２〜５倍減少させることが可能になり、次いで、前の研究で示したフィルタリング結果より、さらに良好な結果が得られる。
【００４４】
実施例
ＤＢＳ手術中に５人のヒトの患者から収集された８つの記録が、アルゴリズムをテストするために使用された。神経細胞活性が、視床下核、淡蒼球または黒質内の高周波刺激の前、その間およびその後に記録されている（いくつかの記録については、刺激後のセグメントが、利用できない）。メドトロニック（Medtronic）社製の外部刺激装置を使用して生成された刺激および神経細胞信号（ボルトを単位にして増幅後に測定）を、AlphaOmega Microguide（登録商標）記録システムを使用して捕捉した。記録は、４８ｋＨｚのサンプリング・レートを使用し、期間が９〜１００秒間であり、刺激は、７〜８０秒間、持続した。刺激は、同じ脳核中に２ｍｍだけ離して置かれた微小電極を介し、周波数が１３０Ｈｚ（すなわち、期間が、ほぼ３７０個のサンプリングきざみ）で、パルス幅が６０μｓで、周期的に送られた。パルス電流強さは、１０００μＡ〜６０００μＡの範囲内であった。
アーチファクト振幅に依存したアルゴリズムの振る舞いを詳しく調査するために、ＡＳＲが事前に定められた人為的な信号のセットを構築した。神経細胞活性がない場合のアーチファクトを含む記録（６秒セグメント）には、単一神経細胞のスパイクを有し刺激がない場合の他の信号を加えた。適切な定数によって、ＡＳＲ所望の値を与える。
X Synthetic(t)=X Artefacts(t)+KX Single-Neuron(t)
ＡＳＲが２〜２０の範囲にある、そのような信号のコレクションを生成した。
式xSt(t+kT)=x0St(t)+ξk(t)、0<t≦Tによる加法ノイズを伴う標準モデルと、式
【数２３】
2010524612.tif 000025
によるモデルによる信号記述の性質を比較するために、残渣のしかるべき位相依存の標準偏差STD(t)、0<t≦Tを図１６に示すように計算し、図１６では、破線が、
式xSt(t+kT)=x0St(t)+ξk(t)、0<t≦T
による加法ノイズを伴う標準モデルによる信号を表し、実線が、式
【数２４】
2010524612.tif 000026
による非線形振動モデルによる信号の表現に基づいた、位相依存のアーチファクトのＳＴＤを表す。
式xSt(t+kT)=x0St(t)+ξk(t)、0<t≦T
によるモデル中の加法ノイズξk(t)の標準偏差σζ(t)、0<t≦T、および式
【数２５】
2010524612.tif 000027
による垂直偏差のベクトルの第１の座標n1(t)の標準偏差σn1(θ)、0<θ≦Tを見積もった。値n1(t)は、式
【数２６】
2010524612.tif 000028
によるモデルについて、時間領域中のアーチファクトの残渣を表す。
８つの実際のデータ記録は、上記に述べたように、STD(t)を見積もるために、使用された。式
【数２７】
2010524612.tif 000029
による位相モデルは、特にピークの近くで、より一層小さい分散を生成する（図１７参照）。位相モデルの場合、残渣全体にわたるＳＴＤの和は、平均して２．１３倍小さく、ＳＴＤの最大値は、１．５〜２．１（平均して１．８）倍小さかった。
【００４５】
このテストは、刺激があり洗浄手順中に、スパイク列およびスパイク自体に何が起きるかを明らかにすることを意図している。アーチファクトの混入の前とすべての処理の終了後に神経細胞活性のオリジナル信号X Single-Neuron(t)から得られたスパイク列を比較した。スパイク列の抽出のために、単純な閾値フィルタを使用した。処理された信号中でそれらの正しい位置で見つかったスパイクのパーセンテージを測定した。係数Ｋを用いた信号の合計によってスパイク・ノイズ比が変わることに注意されたい。図１８および図２０に、ＡＳＲに依存しＳＮＲに対応するスパイク検出のエラーのグラフを表す。図１９は、ＡＳＲに依存したスパイクの２つの第１の主な成分の分散を表す。スパイクは、約１０〜１２のＡＳＲまで事実上変わらないままであることを見ることができる。さらに、ＡＳＲ＝０が、刺激を加えられていない場合のオリジナルの神経細胞活性に対応することに注意されたい。
【００４６】
上記に述べた記録は、実際のデータに対してアルゴリズムをテストするために、使用された。信号は、ASR∈[1.3、14]、ANR∈[3、160]およびSNR∈[1.4、5.1]によって特徴付けられる。図３〜６に、オリジナル信号およびフィルタリングされた信号の例を示す。図７〜１０は、刺激期間のウィンドウのセットにスライスされた、フィルタリングの前およびその後の信号を示す。これらの図は、オリジナルのアーチファクトおよび処理が完全に終わった後に残された残渣の形状を表す。神経細胞のスパイクも見ることができる。図１２および１３に、しかるべきアーチファクトのセットの相図を３Ｄで示す。
【００４７】
図２１に、位相空間および時間領域の両方中での信号フィルタリングの前とその後のアーチファクト・スパイク比（ＡＳＲ）の間の依存関係を示す。比較のために、合成データについて、そのような結果も示す。８つの実際の記録および合成信号についての結果は、良好にフィットしている。位相空間フィルタリングの使用は、実際のデータおよび合成データに関するすべての実験で、アーチファクト残渣をより少なく生成する。すなわち、ASR Time/ASR Phase Space>1である。ここで、ASR Phase SpaceおよびASR Timeは、位相空間中および時間領域中での処理の後のＡＳＲである。実際のデータ・セット（初期ASR∈[1.3、14]）の場合、平均してASR Phase Space=0.239およびASR Time=0.537であり、合成データ（初期ASR∈[2、20]）の場合、平均してASR Phase Space=0.418およびASR Time=1.151である。位相空間フィルタリングの効果は、信号の初期ＡＳＲが増加によって高まる。実際のデータおよび合成データの両方に関する実験で、ＡＳＲがフィルタリング前に４を超える場合、アーチファクト残渣の比が、２を超える、すなわちASR Time/ASR Phase Space>2である。すなわち、実際のデータ・セットについて、ASR>4の場合、平均でASR Phase Space/ASR Time=2.87であり、合成データについて、平均でASR Phase Space/ASR Time=2.81である。
【００４８】
したがって、脳‐コンピュータ・インターフェース中で神経細胞信号を処理するための方法およびシステムを上記に述べてきた。ここで述べた個々のシステム構成要素は、本発明の様々なプロセスおよびオペレーションを実施するために、具体的に構築することができ、あるいは、それらは、必要な機能を備えるためにプログラム・コードによって選択的に起動される、または再構成される汎用コンピュータまたは計算プラットフォームを含むことができる。本明細書に開示されたプロセスは、本質的に、いかなる特定のコンピュータまたは装置にも関連付けられなく、ハードウェア、ソフトウェアおよび／またはファームウェアのいずれの適切な組み合わせによって実施することができる。たとえば、様々な多目的のマシーンは、本発明の教示によって書かれたプログラムとともに使用することができ、あるいは、専門化された装置またはシステムを構築し、それによって必要な方法または技術を実施することは、さらに都合がよいはずである。
【００４９】
本発明の他の実施形態は、本明細書を検討し、ここに述べられた本発明を実施することによって当業者に明らかになる。本明細書および実施例は、単に例示するものとして見なし、本発明の真の範囲および精神は、以下のクレームによって指示されるべきと意図するものである。
【図１】
2010524612.tif 000030
【図２】
2010524612.tif 000031
【図３】
2010524612.tif 000032
【図４】
2010524612.tif 000033
【図５】
2010524612.tif 000034
【図６】
2010524612.tif 000035
【図７】
2010524612.tif 000036
【図８】
2010524612.tif 000037
【図９】
2010524612.tif 000038
【図１０】
2010524612.tif 000039
【図１１】
2010524612.tif 000040
【図１２】
2010524612.tif 000041
【図１３】
2010524612.tif 000042
【図１４】
2010524612.tif 000043
【図１５】
2010524612.tif 000044
【図１６】
2010524612.tif 000045
【図１７】
2010524612.tif 000046
【図１８】
2010524612.tif 000047
【図１９】
2010524612.tif 000048
【図２０】
2010524612.tif 000049
【図２１】
2010524612.tif 000050
【国際調査報告】
2010524612.tif 000051
2010524612.tif 000052
2010524612.tif 000053