本発明は、神経細胞から測定される生体電気信号を処理するスパイクソーティング装置、スパイクソーティング方法、スパイクソーティングプログラム、および記録媒体に関する。

神経科学分野における強力な計測法として、マルチニューロン活動記録法が幅広く使用されている。 この方法では、近接する多数の神経細胞（ニューロン）から、スパイク波形を同時記録する。 具体的には、多数の記録電極から構成される１本の多点電極を用意し、この多点電極を用いて、狭い範囲に存在するほぼ全ての単一ニューロン活動波形を同時に検出する。

現在、自作可能で最も広く活用されている多点電極に、４本の記録電極を束ねたテトロード（ｔｅｔｒｏｄｅ）がある。 テトロードなどの多点電極を用いてスパイク波形を記録する場合、スパイク波形を表すスパイク波形データを処理することによって、スパイク波形から単一ニューロン活動波形を分離し、抽出することが求められる。 これを一般に、スパイクソーティングと呼ぶ。 スパイクソーティングの段階において、スパイク波形データの取り扱いを誤ると、スパイク波形データに秘められている機能的情報を的確に読み取ることができない。 そればかりか、事実とは異なる誤った結論に到達する可能性もある。

近接する複数の神経細胞は、時に、同時に発火する。 これにより、発生したスパイク波形にオーバーラップが生じる。 オーバーラップが生じると、単一ニューロン活動波形を正確に記録できない。 しかし、既存のスパイクソーティング方法は、このようなスパイク波形のオーバーラップに対処できない。 さらには、スパイク波形の変動などの問題にも対処できない。 したがって、近接した神経細胞間の詳細な相互作用、すなわち、数ミリ秒の時間内に神経細胞間に作用している関係性を検出できない。

このような問題に対処するため、非特許文献１には、マルチニューロン活動波形を適当な数のクラスタに分割する手法であるｋ−ｍｅａｎｓクラスタリングと、観測された信号から原波形を抽出する独立成分分析法とを組み合わせたスパイクソーティング方法が開示されている。

非特許文献１によれば、生波形の記録に用いる記録電極の本数よりも多い数の単一ニューロン活動波形を分離できる。 したがって、スパイク波形のオーバーラップや変動などの、既存のスパイクソーティングでは対処できない問題に対処できる。 なお、この方法は、非特許文献２〜４にも開示されている。
独立成分分析を応用したマルチニューロン活動の解析、脳の科学、２５巻、５３〜６０頁、２００３年 Ｃｌａｓｓｉｆｉｃａｔｉｏｎ ｏｆ ｎｅｕｒｏｎａｌ ａｃｔｉｖｉｔｉｅｓ ｆｒｏｍ ｔｅｔｒｏｄｅ ｒｅｃｏｒｄｉｎｇｓ ｕｓｉｎｇ ｉｎｄｅｐｅｎｄｅｎｔ ｃｏｍｐｏｎｅｎｔ ａｎａｌｙｓｉｓ、４９巻、２８９〜２９８頁、２００２年 Ａ ｎｅｗ ａｐｐｒｏａｃｈ ｔｏ ｓｐｉｋｅ ｓｏｒｔｉｎｇ ｆｏｒ ｍｕｌｔｉ−ｎｅｕｒｏｎａｌ ａｃｔｉｖｉｔｉｅｓ ｒｅｃｏｒｄｅｄ ｗｉｔｈ ａ ｔｅｔｒｏｄ − ｈｏｗ ＩＣＡ ｃａｎ ｂｅ ｐｒａｃｔｉｃａｌ、Ｎｅｕｒｏｓｃｉｅｎｃｅ Ｒｅｓｅａｒｃｈ、４６巻、２６５〜２７２頁、２００３年 Ａｕｔｏｍａｔｉｃ Ｓｏｒｔｉｎｇ ｆｏｒ Ｍｕｌｔｉ−Ｎｅｕｒｏｎａｌ Ａｃｔｉｖｉｔｙ Ｒｅｃｏｒｄｅｄ Ｗｉｔｈ Ｔｅｔｒｏｄｅｓ ｉｎ ｔｈｅ Ｐｒｅｓｅｎｓｅ ｏｆ Ｏｖｅｒｌａｐｐｉｎｇ Ｓｐｉｋｅｓ、Ｊ Ｎｅｕｒｏｐｈｙｓｉｏｌ、８９巻、ｐ２２４５〜２２５８頁、２００３年

しかし、上述した従来技術には、処理速度が遅く、分離精度が低いという問題がある。 これは、ｋ−ｍｅａｎｓクラスタリングを用いて、スパイク波形のクラスタリングを行うからである。

ｋ−ｍｅａｎｓクラスタリングは、神経細胞における単一ニューロン活動と異なる他の生体電気現象など、一般的なデータをクラスタリングするための汎用的な手法である。 ニューロン活動に対して、最適化されているわけではない。 このため、ニューロン活動の特性を考慮せず一連の処理を数十回繰り返してしまう。 このことから、単一ニューロン活動波形の分離に、大幅に長い時間を要する。 また、正規分布にしたがってクラスタ境界を設定するため、正規分布にしたがわない特性のパラメータを有するスパイク波形を誤って単一ニューロン活動波形に分離してしまう。 このことから、単一ニューロン活動波形の分離精度が低くなる。

本発明は上記の課題を解決するためになされたものであり、その目的は、より高い分離精度で、かつ、より短い時間で単一ニューロン活動波形を分離抽出できるスパイクソーティング装置、スパイクソーティング方法、スパイクソーティングプログラムおよび記録媒体を提供することにある。

本発明に係るスパイクソーティング装置は、上記の課題を解決するために、複数の記録電極から構成される多点電極によって神経細胞から測定された複数の生波形から、単一ニューロン活動波形を分離するスパイクソーティング装置であって、上記複数の生波形から、スパイク波形を抽出するスパイク波形抽出手段と、上記スパイク波形から、期待値最大化法に基づき、複数のクラスタ波形を生成するクラスタ波形生成手段と、上記複数のクラスタ波形から、それぞれ、独立成分分析法に基づき、独立成分波形を生成する独立成分波形生成手段と、上記複数の独立成分波形のうち、互いに類似する波形を、所定の条件が成立する場合に連結することによって、単一ニューロン活動波形を抽出する単一ニューロン活動波形抽出手段とを備えていることを特徴としている。

上記の構成によれば、本装置は、テトロード等の多点電極によって神経細胞から測定される複数の生波形から、単一ニューロン活動波形を分離する。 このとき、本装置は、生波形からスパイク波形を抽出し、抽出したスパイク波形から、期待値最大化法に基づき、クラスタ波形を生成する。 この期待値最大化法として、たとえば、ＫｌｕｓｔａＫｗｉｋのアルゴリズムを使用する。 また、本装置は、複数のクラスタ波形から、独立成分分析法に基づき、独立成分波形を生成する。 この独立成分分析法として、たとえば、ＦａｓｔＩＣＡのアルゴリズムを使用する。 さらに、上記複数の独立成分波形のうち、互いに類似する波形を、所定の条件、たとえば、連結後の波形に不応期が存在する場合などに、互いに連結する。

以上のように、本装置は、マルチニューロン活動波形から、単一ニューロン活動波形を分離抽出するための方法を利用している。 そのため、クラスター境界を、正規分布に厳密にしたがって決定するわけではない。 ニューロン活動の場合、距離の非常に近いクラスターが、同じ神経細胞から発生していると経験的に仮定されている。 そのため、本装置では、それらを同一クラスターとして連結する。

このように、本装置は、ニューロン活動の特性を考慮にいれてクラスタリングする。 これにより、ｋ−ｍｅａｎｓクラスタリングを利用している従来技術に比べて、単一ニューロン活動波形の分離精度をより向上できる効果を奏する。 また、本装置は、ＡＩＣ基準（Ａｋａｉｋｅ'ｓ Ｉｎｆｏｒｍａｔｉｏｎ Ｃｒｉｔｅｒｉａ；赤池の情報量基準）などを利用して、クラスタを生成する際の演算を無駄に繰り返さない。 そのため、ｋ−ｍｅａｎｓクラスタリングを利用している従来技術に比べて、単一ニューロン活動波形をより高速に分離できる効果を奏する。

本発明に係るスパイクソーティング方法は、上記の課題を解決するために、複数の記録電極から構成される多点電極によって神経細胞から測定された複数の生波形を処理するスパイクソーティング装置において、単一ニューロン活動波形を分離するスパイクソーティング方法であって、上記複数の生波形から、スパイク波形を抽出するスパイク波形抽出ステップと、上記スパイク波形から、期待値最大化法に基づき、複数のクラスタ波形を生成するクラスタ波形生成ステップと、上記複数のクラスタ波形から、それぞれ、独立成分分析法に基づき、複数の独立成分波形を抽出する独立成分波形抽出ステップと、上記複数の独立成分波形のうち、互いに類似する波形を、所定の条件が成立する場合に連結することによって、上記単一ニューロン活動波形を抽出する単一ニューロン活動波形抽出ステップとを含んでいることを特徴としている。

この方法によれば、上述したスパイクソーティング装置と同様の作用効果を奏する。

また、本発明に係るスパイクソーティング装置では、さらに、上記所定の条件は、上記複数の独立成分波形のうち２つを連結した場合に、連結後の波形に所定数の不応期が存在することであることが好ましい。 この構成によれば、スパイクソーティング装置は、連結後の独立成分波形に所定数の不応期が存在する場合に、独立成分波形から単一ニューロン活動波形を抽出する。 そのため、不応期が存在しない波形を、すなわち、ノイズとして取り扱うべき波形を、あやまって単一ニューロン活動波形として抽出する可能性が減少する。 これにより、本構成のスパイクソーティング装置は、単一ニューロン活動波形の分離精度をよりいっそう高めることができる効果を奏する。

なお、上記スパイクソーティング装置は、コンピュータによって実現してもよい。 この場合、コンピュータを上記各手段として動作させることにより上記スパイクソーティング装置をコンピュータにおいて実現するスパイクソーティングプログラム、およびそのスパイクソーティングプログラムを記録したコンピュータ読み取り可能な記録媒体も、本発明の範疇に入る。

以上のように、本発明に係るスパイクソーティング装置は、スパイク波形から、期待値最大化法に基づき、複数のクラスタ波形を生成するクラスタ波形生成手段と、上記複数のクラスタ波形から、それぞれ、独立成分分析法に基づき、複数の独立成分波形を生成する独立成分波形生成手段と、上記複数の独立成分波形のうち、互いに類似する波形を、所定の条件が成立する場合に連結することによって、単一ニューロン活動波形を分離する単一ニューロン活動波形抽出手段とを備えているため、より高い分離精度で、かつ、より短い時間で単一ニューロン活動波形を分離抽出できる効果を奏する。

本発明の一実施形態について、図１〜図１１を参照して以下に説明する。

本発明のスパイクソーティング装置１は、ｎ（ｎは２以上の整数）個の記録電極を含む多点電極を用いて、神経細胞から測定されるｎ個の生波形データを処理する装置である。 ユーザは、テトロードなどの多点電極を用いて、神経細胞からマルチニューロン活動を記録する。 そして、記録したマルチニューロン活動を、スパイクソーティング装置１に入力する。 これにより、スパイクソーティング装置１は、入力されたマルチニューロン活動の生波形データを処理することによって、生波形から単一ニューロン活動波形を抽出する。

まず、スパイクソーティング装置１が処理する生波形データの取得方法を簡単に説明する。 本実施形態では、４本のタングステンワイヤからなるテトロードを使用する。 ユーザは、テトロードを、可動式慢性埋め込み用のマイクロドライブに取り付ける。 このマイクロドライブを操作して、ユーザは、テトロードをたとえばサルの脳内に挿入する。 挿入後、１時間に渡り、サルの神経細胞から、テトロードを通じてマルチニューロン活動を記録する。 ４本のタングステンワイヤのそれぞれから、生波形データを１つずつ記録する。 したがって、１本のテトロードから、各タングステンワイヤに由来する４つの生波形データを記録できる。

タングステンワイヤのそれぞれから記録された生波形信号を、増幅器に入力する。 増幅器は、入力された生波形信号を、２０００倍に増幅し、増幅後の生波形信号をＡＤ変換器に入力する。 ＡＤ変換器は、入力された生波形信号を、２５ｋＨｚのサンプリングレートで、デジタルデータに変換する。 さらに、ＡＤ変換器は、デジタル化により生成した生波形データを、スパイクソーティング装置１に入力する。 スパイクソーティング装置１は図示しない記憶部を備えており、入力された生波形データを、この記憶部にいったん保存する。

スパイクソーティング装置１は、入力された生波形データを処理することによって、生波形から単一ニューロン活動波形を分離する。 なお、本実施形態では、４本のタングステンワイヤを通じて、４つの生波形データが記録される。 そこで、これらの４つの生波形データを、それぞれ、第１〜第４チャンネルの生波形データと表現する。 すなわち、第１記録電極を通じて取得した生波形データを、第１チャンネルの生波形データ（第１生波形データ）と表現する。 また、第２記録電極を通じて取得した生波形データを、第２チャンネルの生波形データ（第２生波形データ）と表現する。 また、第３記録電極を通じて取得した生波形データを、第３チャンネルの生波形データ（第３生波形データ）と表現する。 また、第４記録電極を通じて取得した生波形データを、第４チャンネルの生波形データ（第４生波形データ）と表現する。

本実施形態におけるスパイクソーティング装置１の構成について、図１を参照して以下に説明する。 図１は、本発明の一実施形態に係るスパイクソーティング装置１の構成を示すブロック図である。 この図に示すように、スパイクソーティング装置１は、スパイク波形抽出部１０（スパイク波形抽出手段）、スパイク特徴抽出部１１（スパイク特徴抽出手段）、クラスタ波形生成部１２（クラスタ波形生成手段）、独立成分波形抽出部１３（独立成分波形抽出手段）、再構成波形生成部１４（再構成波形生成手段）、および単一ニューロン活動波形抽出部１５（単一ニューロン活動波形抽出手段）を備えている。

スパイク波形抽出部１０は、生波形からスパイク波形を抽出する。 このとき、スパイク波形抽出部１０は、生波形においてスパイク波形が出現する時間を算出することによって、この時間をタイムスタンプとする。

スパイク特徴抽出部１１は、スパイク波形から、スパイク波形の特徴を表すスパイク特徴値を抽出する。 たとえば、スパイク特徴抽出部１１は、スパイク波形の正ピーク値および負ピーク値を、スパイク特徴値として抽出する。

クラスタ波形生成部１２は、スパイク特徴値に基づき、スパイク波形をクラスタリングする。 このとき、クラスタ波形生成部１２は、ある同一のクラスタに属すると判定した複数のスパイク波形を連結することによって、連結クラスタ波形群を生成する。 なお、クラスタ波形生成部１２は、期待値最大化法に基づき、スパイク波形をクラスタリングする。

再構成波形生成部１４は、連結クラスタ波形群から独立成分クラスタ波形を生成する。 このとき、再構成波形生成部１４は、独立成分分析法に基づき、独立成分クラスタ波形を生成する。

単一ニューロン活動波形抽出部１５は、独立成分クラスタ波形から単一ニューロン活動波形を抽出する。 このとき、単一ニューロン活動波形抽出部１５は、２つの独立成分クラスタ波形を連結した波形に、不応期が存在する場合に、単一ニューロン活動波形候補を生成する。 さらに、単一ニューロン活動波形抽出部１５は、不応期が存在すると判定した単一ニューロン活動波形候補を、単一ニューロン活動波形とする。

スパイクソーティング装置１が、入力された生波形データを処理することによって、単一ニューロン活動波形を表す単一ニューロン活動波形データを生成する際のデータ処理の流れについて、以下に詳細に説明する。 図１に示すように、スパイクソーティング装置１のスパイク波形抽出部１０に、生波形を表す生波形データが入力される。 スパイクソーティング装置１が処理する生波形データは、記録された時間順に電圧値（観測信号値）が並んでいる配列データである。 ここでいう「時間」とは、生波形データの記録を開始してからの経過時間を意味する。

スパイク波形抽出部１０は、入力された生波形データに基づき、所定の演算を実行して、生波形から所定の閾値を算出する。 本実施形態では、スパイク波形抽出部１０は、まず、生波形の平均２乗誤差を算出する。 具体的には、スパイク波形抽出部１０は、生波形データに含まれている所定数の連続した電圧値を用いて、生波形の平均２乗誤差を算出する。 たとえばスパイク波形抽出部１０は、生波形データの先頭から１０分間分の電圧値を使用し、平均２乗誤差を算出する。 なお、平均２乗誤差の算出方法は公知であるため、詳細な説明を省略する。

スパイク波形抽出部１０は、４チャンネルの各生波形データから、それぞれ、対応する平均２乗誤差を算出する。 たとえば、第１生波形データ（第１チャンネルの生波形データ）から、第１平均２乗誤差を算出する。 また、第２生波形データ（第２チャンネルの生波形データ）から、第２平均２乗誤差を算出する。 また、第３生波形データ（第３チャンネルの生波形データ）から、第３平均２乗誤差を算出する。 また、第４生波形データ（第４チャンネルの生波形データ）から、第４平均２乗誤差を算出する。 スパイク波形抽出部１０は、算出した平均２乗誤差を、それぞれ対応する生波形データに関連づける。 さらにスパイク波形抽出部１０は、算出した平均２乗誤差をいったん記憶部に保存する。

スパイク波形抽出部１０は、算出した平均２乗誤差に５を乗算することによって、それぞれをいずれも５倍の値に変更する。 これにより、スパイク波形抽出部１０は、第１〜第４生波形データから、それぞれに対応する第１〜第４閾値を算出する。 これらの閾値を算出したあと、スパイク波形抽出部１０は、生波形データに含まれる各電圧値を、先頭から順に、それぞれ対応する閾値と比較する。

スパイク波形抽出部１０は、第１〜第４生波形データのうち、少なくともいずれかの生波形データに含まれるある電圧値が、比較した閾値を超えるか否かを判定する。 ここで、判定結果が真である場合、スパイク波形抽出部１０は、閾値を超えると判定した電圧値が生波形において出現する時間を、タイムスタンプとして記憶部に保存する。 たとえば、第１〜第４生波形データのうち、第４生波形データに含まれる電圧値が、第４閾値を超えたと判定する場合、スパイク波形抽出部１０は、その電圧値が第４生波形データにおいて出現する時間を、第１〜第４生波形データに共通するタイムスタンプとして、記憶部に保存する。

スパイク波形抽出部１０は、算出したタイムスタンプに基づき、第１〜第４生波形データの全てから、それぞれ、スパイク波形を抽出する。 この例を、図２に示す。 図２は、スパイクソーティング装置１が処理した生波形、および、生波形から抽出したスパイク波形の一例を示す図である。

図２に示すように、スパイク波形抽出部１０は、４チャンネルの生波形データから、スパイク波形群ｆ１を抽出する。 スパイク波形群ｆ１の場合、第４チャンネルの生波形データにおける電圧値が第４閾値を超えている。 そこでスパイク波形抽出部１０は、この電圧値が第４生波形において現れる時間を、タイムスタンプとして算出する。 つぎにスパイク波形抽出部１０は、算出したタイムスタンプを基準にして、第１〜第４生波形データの全てから、それぞれ、スパイク波形を抽出する。 すなわち、スパイク波形抽出部１０は、第１〜第４生波形データにおいていずれも同じ時間に出現するスパイク波形を抽出する。

このときスパイク波形抽出部１０は、図３に示すように、抽出したスパイク波形の１６ポイント目がスパイク波形の負ピーク値になるようにする。 たとえば、スパイク波形抽出部１０は、各生波形データから、負ピーク値を中心に、前後１６ポイント分のデータを取り出す。 ここでいう「ポイント」とは、スパイク波形を構成する電圧値の数を意味する。

この操作を、スパイク波形抽出部１０は、４チャンネルの全ての生波形データに対して適用する。 したがって、スパイク波形抽出部１０は、算出したタイムスタンプに基づき、４つの生波形データから、合計して４つのスパイク波形を抽出する。 スパイク波形抽出部１０は、抽出される４つのスパイク波形と、これらのスパイク波形の抽出に用いたタイムスタンプとを関連づけて、スパイク波形群ｆ１を生成する。 生成したスパイク波形群ｆ１を、スパイク波形抽出部１０は、記憶部に保存する。

スパイク波形抽出部１０は、最初のスパイク波形群ｆ１の分離処理が終了してから、あらたなスパイク波形の分離処理を開始する。 すなわち、抽出したスパイク波形に含まれる最後の電圧値の次の電圧値から、再び、閾値との比較を再開する。 もし、スパイク波形群ｆ１を抽出したデータと同様に、４チャンネルのうち少なくともいずれかの生波形データに含まれる電圧値が、閾値を超えると判定する場合、その電圧値が出現するタイムスタンプに基づき、あらたなスパイク波形を抽出する。 図２の例では、スパイク波形抽出部１０は、スパイク波形群ｆ１を抽出したあと、スパイク波形群ｆ２を抽出する。 このとき、スパイク波形群ｆ２を抽出するためのタイムスタンプを算出し、タイムスタンプに基づき各生波形データからスパイク波形を抽出して、スパイク波形群ｆ２を生成する。 このようにして、スパイク波形抽出部１０は、図２に示す生波形データから、図３に示すように、スパイク波形群ｆ１およびスパイク波形群ｆ２を生成する。

スパイク波形抽出部１０は、生波形データの最初から最後まで、このようなスパイク波形の抽出処理を続ける。 たとえば、１時間にわたって神経細胞から測定した生波形データを処理した場合、スパイク波形抽出部１０は、４チャンネルの生波形データから、それぞれ、数万個のスパイク波形を抽出する。 スパイク波形抽出部１０は、抽出したスパイク波形を、いずれも、スパイク波形の抽出に用いたタイムスタンプと関連づけることによって、スパイク波形群を生成する。

以上の処理によって、スパイク波形抽出部１０は、生波形から複数のスパイク波形群を抽出する。 スパイク波形抽出部１０は、スパイク波形群を表すスパイク波形群データをスパイク特徴抽出部１１に出力する。

スパイク特徴抽出部１１は、入力されたスパイク波形群データを処理することによって、スパイク波形群からスパイク特徴値を抽出する。 この例について、図４を参照して以下に説明する。 図４は、スパイク波形から抽出されるスパイク特徴値の一例を示す図である。

図４に示すように、スパイク特徴抽出部１１は、スパイク特徴値として、スパイク波形群に含まれる各スパイク波形から、それぞれ、正ピーク値および負ピーク値を算出する。 このとき、スパイク波形から正ピーク値および負ピーク値を抽出するアルゴリズムは任意のものでよい。 図４に示す例では、スパイク特徴抽出部１１は、スパイク波形群ｆ１の第１チャンネルのスパイク波形から、正ピーク値として１１を算出し、負ピーク値として−２２を算出する。 また、第２チャンネルのスパイク波形から、正ピーク値として１０を算出し、負ピーク値として−２２を算出する。 また、第３チャンネルのスパイク波形から、正ピーク値として１０を算出し、負ピーク値として−１８を算出する。 また、第４チャンネルのスパイク波形から、正ピーク値として１５を算出し、負ピーク値として−２４を算出する。

このように、スパイク特徴抽出部１１は、１つのスパイク波形群に含まれる４つのスパイク波形から、４つの正ピーク値と、４つの負ピーク値とを算出する。 つぎにスパイク特徴抽出部１１は、算出したスパイク特徴値を、一つのスパイク特徴値ベクトルにまとめる。 たとえば、スパイク波形群ｆ１から、（１１、−２２、１０、−２２、１０、−１８、１５、−２４）のように、８次元のスパイク特徴値ベクトルｖ１を生成する。 スパイク特徴抽出部１１は、この８次元のスパイク特徴値ベクトルｖ１を表すスパイク特徴値ベクトルデータを生成し、スパイク波形群ｆ１と関連づけ、記憶部に保存する。 また、スパイク特徴抽出部１１は、このスパイク特徴値ベクトルデータをクラスタ波形生成部１２に出力する。

スパイク特徴抽出部１１は、図４に示すように、スパイク波形群ｆ１の場合と同様に、スパイク波形群ｆ２の各チャンネルのスパイク波形からも正ピーク値および負ピーク値を算出する。 これにより、スパイク特徴抽出部１１は、図５に示すように、スパイク波形群ｆ２から、（２０、−２２、２９、−３３、１２、−１４、１５、−１９）のように、８次元のスパイク特徴値ベクトルｖ２を生成する。 スパイク特徴抽出部１１は、抽出したスパイク特徴値ベクトルｖ２を表すスパイク特徴値ベクトルデータを生成し、スパイク波形群ｆ２に関連づける。 さらに、このスパイク特徴値ベクトルデータをクラスタ波形生成部１２に出力する。

スパイク特徴抽出部１１は、入力された全てのスパイク波形群データから、８次元のスパイク特徴値ベクトルを算出する。 これにより、スパイク特徴抽出部１１は、入力されたスパイク波形群データの数と同一の数のスパイク特徴値ベクトルデータを生成し、これらをいずれも、クラスタ波形生成部１２に出力する。 たとえば、スパイク波形群が１０，０００個ある場合、スパイク特徴抽出部１１は、１０，０００個のスパイク特徴値ベクトルデータを生成し、クラスタ波形生成部１２に出力する。 この場合、スパイク特徴抽出部１１は、１０，０００個のスパイク特徴値ベクトルデータのそれぞれを、いずれも、スパイク特徴値ベクトルデータの生成に用いたスパイク波形群データに関連づける。

クラスタ波形生成部１２は、スパイク特徴値ベクトルに基づき、スパイク波形群をクラスタリングする。 具体的には、入力されたスパイク特徴値ベクトルデータを用いて、所定の演算を実行して、スパイク特徴値ベクトルデータをクラスタリングする。 実際には、スパイク特徴値ベクトルデータに関連づけられているスパイク波形群データのクラスタ番号を算出することによって、このクラスタ番号を、スパイク特徴値ベクトルデータに関連づける。 クラスタ波形生成部１２は、スパイク特徴値ベクトルデータに関連づけられているクラスタ番号に基づき、スパイク特徴値ベクトルデータに関連づけられているスパイク波形群データをクラスタリングする。 すなわち、クラスタ波形生成部１２は、スパイク波形群データのクラスタを生成する。

このとき、クラスタ波形生成部１２は、予め定められた初期値に基づき、生成するクラスタの数を決定する。 初期値として、たとえば、上限として３０個、下限として１０個が設定されていればよい。 これらの初期値は、図示しない記憶部に予め格納されている。 そのため、クラスタ波形生成部１２は、スパイク波形群データをクラスタリングするとき、記憶部にアクセスし、初期値を読み出す。

なお、ユーザが、この初期値を、スパイクソーティング装置１の図示しない入力部を介して入力することによって、スパイクソーティング装置１に、クラスタの上限および下限を指示してもよい。 このとき、クラスタ波形生成部１２は、ユーザによって入力された値を、クラスタ個数の初期値として使用する。

クラスタ波形生成部１２は、期待値最大化法に基づき、スパイク波形群データをクラスタリングする。 たとえば、クラスタ波形生成部１２は、各スパイク特徴値ベクトルが、８次元空間において互いに近接した位置に存在するか否かを判定する。 このとき、クラスタ波形生成部１２は、あるスパイク特徴値ベクトルと、他のスパイク特徴値ベクトルとが、８次元空間において、予め定められた所定の距離をおいて位置するか否かを判定する。 ２つのスパイク特徴値ベクトルが、所定の距離よりも近い距離に位置すると判定する場合、クラスタ波形生成部１２は、これらのスパイク特徴値ベクトルに同一のクラスタ番号を設定する。

クラスタ波形生成部１２は、入力される全てのスパイク特徴値ベクトルの組み合わせに対して、上述した判定を行い、スパイク波形群データをクラスタリングする。 これにより、クラスタ波形生成部１２は、入力されたスパイク特徴値ベクトルの全てにクラスタ番号を設定する。 すなわち、各スパイク特徴値ベクトルと、このスパイク特徴値ベクトルが属するクラスタの番号とを関連づけて、記憶部に保存する。 これにより、スパイク特徴値ベクトルに関連づけられているスパイク波形群に、クラスタ番号を関連づける。

なお、クラスタ波形生成部１２は、ＫｌｕｓｔａＫｗｉｋにおいて採用されているアルゴリズムに基づき、スパイク波形群データをクラスタリングする。 ＫｌｕｓｔａＫｗｉｋは、フリーソフトウェアとして広く流通している（ＵＲＬ：ｈｔｔｐ：／／ｋｌｕｓｔａｋｗｉｋ．ｓｏｕｒｃｅｆｏｒｇｅ．ｎｅｔ／を参照）ものであり、そのアルゴリズ自体は公知であるため、詳細な説明を省略する。

クラスタ波形生成部１２は、スパイク波形群に設定したクラスタ番号に基づき、同一のクラスタ番号のスパイク波形群を、互いに連結する。 たとえば、図６に示すように、２つのスパイク波形群を連結する。 実際には、クラスタ波形生成部１２は、スパイク波形群を表すスパイク波形群データを、互いに関連づけることによって、スパイク波形群を連結する。 このとき、クラスタ波形生成部１２は、タイムスタンプの出現順（すなわち、タイムスタンプの値が小さい順）にスパイク波形群を連結する。 つまり、タイムスタンプがより大きいスパイク波形群を、タイムスタンプがより小さいスパイク波形群の後ろに連結する。 また、クラスタ波形生成部１２は、チャンネル別に、同一のチャンネルのスパイク波形同士を互いに連結する。 たとえば、クラスタ波形生成部１２は、第１チャンネルのスパイク波形同士を連結し、第２チャンネルのスパイク波形同士を連結し、第３チャンネルのスパイク波形同士を連結し、第４チャンネルのスパイク波形同士を連結する。

このようにして、クラスタ波形生成部１２は、複数のスパイク波形群が連結した連結クラスタ波形群を生成する。 また、クラスタ波形生成部１２は、連結した各スパイク波形のタイムスタンプを、互いに連結することによって、一つのタイムスタンプベクトルを生成する。 たとえば、図６に示すように、タイムスタンプｔ１およびｔ２を連結して、ｔ１およびｔ２からなるタイムスタンプベクトルを生成する。 他の例を挙げれば、１０個のスパイク波形群を連結して１つの連結クラスタ波形群を生成する場合は、１０個のタイムスタンプを連結した１つのタイムスタンプベクトルを生成する。 クラスタ波形生成部１２は、生成したタイムスタンプベクトルを、対応する連結クラスタ波形群に関連づける。 さらにクラスタ波形生成部１２は、生成した連結クラスタ波形群を独立成分波形抽出部１３に出力する。

独立成分波形抽出部１３は、入力された連結クラスタ波形群から、独立成分波形を抽出する。 このとき、独立成分波形抽出部１３は、独立成分分析法のアルゴリズムに基づき、独立成分波形を抽出する。 独立成分波形抽出部１３は、１つの連結クラスタ波形群から、４つの独立成分波形を抽出する。 具体的には、連結クラスタ波形群に含まれる４つのスパイク波形から、４つの独立成分波形を抽出する。 すなわち、生波形データの測定に用いたタングステンワイヤの数と同一の数の独立成分波形を抽出する。 独立成分波形抽出部１３が抽出した独立成分波形の一例を、図７に示す。 図７に示す独立成分波形は、図６に示す連結クラスタ波形群から抽出されるものである。

独立成分波形抽出部１３はさらに、連結クラスタ波形群から、１つの分離行列Ｗを算出する。 この分離行列Ｗとして、４行４列のベクトルを算出する。 独立成分波形抽出部１３は、算出した分離行列Ｗを、この分離行列Ｗの生成に使用した連結クラスタ波形群に関連づける。 同時に、算出した分離行列Ｗから、その逆行列である逆行列Ｗ ^−１を算出する。 このとき、算出した逆行列Ｗ ^−１から、独立成分波形ごとに４行１列の基底ベクトルを算出する。 なお、この基底ベクトルは、逆行列Ｗ ^−１のある一行を抜き出した縦ベクトルである。 逆行列Ｗ ^−１は４行４列のベクトルであるため、独立成分波形抽出部１３は、逆行列Ｗ ^−１から、行ごとに合計して４つの基底ベクトルを算出する。

独立成分波形抽出部１３は、生成した分離行列Ｗ、逆行列Ｗ ^−１ 、および４つの基底ベクトルを、いずれも、４つの独立成分波形からなる独立成分波形群に関連づける。 独立成分波形抽出部１３は、分離行列Ｗ等が関連づけられている独立成分波形群を、再構成波形生成部１４に出力する。

なお、独立成分波形抽出部１３は、ＦａｓｔＩＣＡにおいて採用されているアルゴリズムに基づき、スパイク波形をクラスタリングする。 ＦａｓｔＩＣＡは、フリーソフトウェアとして広く流通している（ＵＲＬ：ｈｔｔｐ：／／ｗｗｗ．ｃｉｓ．ｈｕｔ．ｆｉ／ｐｒｏｊｅｃｔｓ／ｉｃａ／ｆａｓｔｉｃａ／を参照）ものであり、そのアルゴリズ自体は公知であるため、詳細な説明を省略する。

再構成波形生成部１４は、つぎに示す式（１）を用いて、独立成分波形群から再構成波形群を生成する。

式（１）において、行列Ｗ' ^−１は、逆行列Ｗ ^−１の中で、指定行以外の要素を全てゼロにした行列である。 すなわち、この行列Ｗ' ^−１は、４行４列の行列である。 また、ｒは、再構成波形生成部１４が生成した再構成波形である。

再構成波形生成部１４は、式（２）を用いた処理によって、１つの独立成分波形群から、４つの再構成波形群を生成する。 すなわち、４つの独立成分波形を含む１つの独立成分波形群から、それぞれ、４つの再構成波形を含む再構成波形群を４つ、生成する。 具体的には、独立成分波形群に、逆行列Ｗ ^−１の１行目以外の要素をゼロにした行列Ｗ' ^−１を適用することによって、第１再構成波形群を生成する。 また、独立成分波形群に、逆行列Ｗ ^−１の２行目以外の要素をゼロにした行列Ｗ' ^−１を適用することによって、第２再構成波形群を生成する。 また、独立成分波形群に、逆行列Ｗ ^−１の３行目以外の要素をゼロにした行列Ｗ' ^−１を適用することによって、第３再構成波形群を生成する。 また、独立成分波形群に、逆行列Ｗ ^−１の４行目以外の要素をゼロにした行列Ｗ' ^−１を適用することによって、第４再構成波形群を生成する。 このとき生成した再構成波形群の一例を、図８（ａ）に示す、図８（ａ）は、第１再構成波形群を示す。

再構成波形生成部１４は、再構成波形群を表す再構成波形群データを生成し、記憶部に保存する。 つぎに再構成波形生成部１４は、再構成波形群に含まれる各チャンネルの再構成波形から、再構成スパイク波形を抽出する。 このとき、再構成波形生成部１４は、各チャンネルの再構成波形の平均２乗誤差を算出する。 再構成波形生成部１４は、算出した平均２乗誤差を５倍し、閾値を算出する。

再構成波形生成部１４は、再構成波形群に含まれている各チャンネルの再構成波形から、先頭から順に３２ポイント分の信号値を取り出す。 取り出した再構成スパイク波形を構成する信号値を、再構成波形生成部１４は、順次、閾値と比較する。 これにより、再構成スパイク波形を構成する信号値が、閾値を超えるか否かを判定する。 閾値を超えると判定する場合、再構成波形生成部１４は、その再構成スパイク波形を取り出し、記憶部に保存する。 閾値を超えないと判定する場合、再構成波形生成部１４は、その再構成スパイク波形を破棄する。

再構成波形生成部１４は、再構成波形群に含まれるいずれかのチャンネルの再構成波形において、信号値が閾値を超えると判定する場合に、全てのチャンネルの再構成波形から、それぞれ、再構成スパイク波形を抽出する。

このようにして、再構成波形生成部１４は、再構成波形群に含まれている全ての再構成スパイク波形の電圧値を、閾値と比較する。 再構成波形生成部１４は、閾値を超える信号値を含む再構成スパイク波形を取り出し、記憶部に保存する。 また、取り出した再構成スパイク波形を、タイムスタンプの出現順に、互いに連結する。 このとき、再構成波形生成部１４は、チャンネル別に、再構成スパイク波形を連結する。 具体的には、ある再構成波形群の第１チャンネルの再構成スパイク波形と、別の再構成波形群の第１チャンネルの再構成スパイク波形とを、互いに連結する。 他のチャンネルについても、同様である。 この例を、図８（ｂ）に示す。 図８（ｂ）の例では、再構成波形生成部１４は、図８（ａ）の再構成波形群から、１、２、３、５、８、９、および１０番目の再構成スパイク波形を抽出する。 これらの抽出した再構成スパイク波形を、再構成波形生成部１４は、チャンネル別に連結することによって、４つの独立成分クラスタ波形を含む独立成分クラスタ波形群を生成する。 実際には、独立成分クラスタ波形群を表す独立成分クラスタ波形群データを生成する。

上述したように、再構成波形群には、タイムスタンプが関連づけられている。 スパイク波形抽出部１０は、このタイムスタンプを、生波形からスパイク波形を抽出するときに算出する。 しかし再構成波形生成部１４は、再構成スパイク波形を、生波形ではなく、独立成分波形から抽出する。 そのため再構成スパイク波形の形状は、スパイク波形のそれとが異なっている可能性が高い。 これにより、再構成スパイク波形のタイムスタンプは、スパイク波形のタイムスタンプと異なる値になる可能性も高い。 したがって、再構成波形群に関連づけられているタイムスタンプを、そのまま再構成スパイク波形のタイムスタンプとして使用する場合、タイムスタンプにずれが生じるおそれがある。 タイムスタンプにずれが生じていると、後述する不応期の判定などに支障が生ずる。

そこで再構成波形生成部１４は、抽出した再構成スパイク波形にあわせて、タイムスタンプを補正する。 具体的には、まず、再構成波形生成部１４は、取り出した各再構成スパイク波形から、再構成スパイク波形の負ピーク値を算出する。 このとき、再構成波形生成部１４は、各チャンネルの再構成スパイク波形から、それぞれ、対応する負ピーク値を算出する。 したがって、４つの再構成スパイク波形から、４つの負ピーク値を算出する。 再構成波形生成部１４は、算出した負ピーク値のうち、もっとも値が大きい負ピーク値が出現する時間を算出する。 これにより、再構成波形生成部１４は、この時間から１６を引いた値を、補正前のタイムスタンプに加算することによって、タイムスタンプを補正する。 この例を、図９に示す。 再構成波形生成部１４は、図９（ａ）にタイムスタンプｔ１〜ｔｈを、図９（ｂ）に示す補正後タイムスタンプｔ'１〜ｔ'ｈに補正する。

再構成波形生成部１４は、独立成分クラスタ波形群に含まれる各再構成スパイク波形のタイムスタンプを、このようにして補正する。 また、補正後のタイムスタンプを、再構成スパイク波形ごとに抽出し、一つにまとめて、補正後タイムスタンプベクトルを生成する。 したがって、補正後タイムスタンプベクトルには、再構成スパイク波形の数と同一数の補正後タイムスタンプが含まれている。

再構成波形生成部１４は、生成した補正後タイムスタンプを、独立成分クラスタ波形群に関連づける。 ここで、上述したように、再構成波形群には、基底ベクトルが関連づけられている。 すなわち、再構成波形生成部１４は、独立成分クラスタ波形群に、補正後タイムスタンプベクトルと、独立成分クラスタ波形群の生成に使用した再構成波形群の基底ベクトルとを、関連づける。 再構成波形生成部１４は、このようにして生成した独立成分クラスタ波形群を、単一ニューロン活動波形抽出部１５に出力する。

単一ニューロン活動波形抽出部１５は、入力された独立成分クラスタ波形群から単一ニューロン活動波形を抽出する。 このとき、単一ニューロン活動波形抽出部１５は、独立成分クラスタ波形群に関連づけられている基底ベクトルに基づき、単一ニューロン活動波形を抽出する。 たとえば、単一ニューロン活動波形抽出部１５は、ある独立成分クラスタ波形群の基底ベクトルと、別の独立成分クラスタ波形群の基底ベクトルとを比較する。 このとき、単一ニューロン活動波形抽出部１５は、次の式（２）を用いて２つの基底ベクトルの角度を算出することによって、２つの基底ベクトルを比較する。

式（２）において、ａおよびｂは、比較対象である２つの基底ベクトルである。 ＜ａ，ｂ＞は、ａおよびｂの内積である。 ｜ａ｜は、ａの絶対値である。 ｜ａ｜は、ｂの絶対値である。 ｜｜ａ｜｜は、ａのノルムである。 ｜｜ｂ｜｜は、ｂのノルムである。 Ｍ _ａｂは、ａおよびｂの角度である。

単一ニューロン活動波形抽出部１５は、算出した角度が、予め定められている所定の角度よりも大きいか否かを判定する。 この判定結果が真である場合、単一ニューロン活動波形抽出部１５は、２つの基底ベクトルに関連づけられている補正後タイムスタンプベクトルを互いに連結する。 さらに単一ニューロン活動波形抽出部１５は、連結した補正後タイムスタンプベクトルに含まれる各補正後タイムスタンプを、値が小さい順に並び替える。

単一ニューロン活動波形抽出部１５は、算出したタイムスタンプベクトルに基づき、連結した独立成分クラスタ波形に、一定数の不応期が存在するか否かを判定する。 具体的には、単一ニューロン活動波形抽出部１５は、まず、補正後タイムスタンプベクトルに含まれる要素（補正後タイムスタンプ）の数を算出する。 つぎに、単一ニューロン活動波形抽出部１５は、補正後タイムスタンプベクトルに含まれる連続した２つのタイムスタンプの間のスパイク間インターバルを算出する。 たとえば、タイムスタンプベクトルに含まれる２番目の補正後タイムスタンプから、最初の補正後タイムスタンプを減算し、第１のスパイク間インターバルを算出する。 さらに、タイムスタンプベクトルに含まれる３番目の補正後タイムスタンプから、２番目の補正後タイムスタンプを減算し、第２のスパイク間インターバルを算出する。 このようにして、単一ニューロン活動波形抽出部１５は、ｎ個のタイムスタンプを要素として持つ補正後タイムスタンプベクトルから、ｎ−１個のスパイク間インターバルを算出する。

つぎに、単一ニューロン活動波形抽出部１５は、ｎ−１個のスパイク間インターバルのうち、所定の値以下のものを算出する。 たとえば、スパイク間インターバルのうち、２ミリ秒以下のものの数を算出する。 さらに、単一ニューロン活動波形抽出部１５は、タイムスタンプの要素の数に、所定の設定値（たとえば０．００１）を乗算する。

単一ニューロン活動波形抽出部１５は、算出した２つの値を比較する。 すなわち、単一ニューロン活動波形抽出部１５は、スパイク間インターバルのうち２ミリ秒以下のものの数が、補正後タイムスタンプベクトルの要素の数に０．００１を乗算した値よりも小さいか否かを判定する。 ここで、小さいと判定する場合、単一ニューロン活動波形抽出部１５は、この補正後タイムスタンプベクトルに関連づけられている２つの独立成分クラスタ波形群を連結した波形に、不応期が存在するとみなす。 これにより、単一ニューロン活動波形抽出部１５は、２つの独立成分クラスタ波形群を互いに連結する。 このとき、単一ニューロン活動波形抽出部１５は、チャンネル別に、２つの独立成分クラスタ波形を互いに連結する。 一方、大きいと判定する場合、単一ニューロン活動波形抽出部１５は、この補正後タイムスタンプベクトルに関連づけられている２つの独立成分クラスタ波形群を連結した波形に、不応期が存在しないとみなす。 したがって、単一ニューロン活動波形抽出部１５は、これらの２つの独立成分クラスタ波形を互いに結合しない。

単一ニューロン活動波形抽出部１５は、独立成分クラスタ波形群から選択される任意の２つの組み合わせの全てについて、上述した不応期の判定を行う。 これにより、単一ニューロン活動波形抽出部１５は、複数の独立成分クラスタ波形群が互いに連結された単一ニューロン活動波形候補を生成し、記憶部に保存する。 さらに、単一ニューロン活動波形抽出部１５は、独立成分クラスタ波形群の補正後タイムスタンプを、互いに連結し、出現順に並び替えることによって、補正後タイムスタンプベクトルを生成し、単一ニューロン活動波形候補に関連づける。 たとえば、単一ニューロン活動波形抽出部１５は、１２個の独立成分クラスタ波形群が互いに連結された単一ニューロン活動波形候補を生成する。 他にも、単一ニューロン活動波形抽出部１５は、４個の独立成分クラスタ波形群が互いに連結した単一ニューロン活動波形候補を生成する。 なお、単一ニューロン活動波形抽出部１５は、他の独立成分クラスタ波形群とまったく連結されなかった独立成分クラスタ波形群も、取り敢えず単一ニューロン活動波形候補として、記憶部に保存する つぎに単一ニューロン活動波形抽出部１５は、単一ニューロン活動波形候補に不応期が存在するか否かを判定する。 すなわち、単一ニューロン活動波形抽出部１５は、単一ニューロン活動波形候補に関連づけられている補正後タイムスタンプベクトルに基づき、上述した不応期の判定処理と同一の判定を行う。 単一ニューロン活動波形抽出部１５は、不応期が存在すると判定した単一ニューロン活動波形候補を、単一ニューロン活動波形として記憶部に保存する。 実際には、単一ニューロン活動波形を表すデータを、記憶部に保存する。 一方、単一ニューロン活動波形抽出部１５は、不応期が存在しないと判定した単一ニューロン活動波形候補を、単一ニューロン活動波形として記憶部に保存しない。 すなわち、単一ニューロン活動波形抽出部１５は、この単一ニューロン活動波形候補を表す単一ニューロン活動波形候補データを、ノイズとして記憶部に保存する。

単一ニューロン活動波形抽出部１５は、入力される全ての単一ニューロン活動波形候補に対し、不応期が存在するか否かを判定する。 これにより、単一ニューロン活動波形抽出部１５は、単一ニューロン活動波形候補から、ノイズ以外の単一ニューロン活動波形のみを選別し、抽出する。

以上のように、スパイクソーティング装置１は、マルチニューロン活動記録法により計測されたデータを処理する。 マルチニューロン活動記録法は、細胞外記録法の一つである。 すなわち、神経細胞における単一ニューロン活動を、直接に同定できるものではない。 そこで、単一ニューロン活動を正確に同定できる細胞内記録によって計測されるスパイク波形と、同時に、細胞外から計測されたマルチニューロン活動とを比較することによって、スパイクソーティング装置１における単一ニューロン活動波形分離の精度を算出した。

なお、単一ニューロン活動波形分離の際のエラー率を、以下のように第一種（Ｍ−Ｘ）／Ｎと第二種（Ｎ−Ｘ）／Ｎの過誤を加算して求めた。

この式において、Ｎは、細胞内記録法により計測された単一ニューロン活動波形の数を表す。 Ｍは、マルチニューロン活動波形に基づき、スパイクソーティング装置１（または従来技術に基づく装置）が分離抽出した単一ニューロン活動波形の数を表す。 Ｘは、細胞内記録法により計測された単一ニューロン活動波形と、スパイクソーティング装置１（または従来技術に基づく装置）が分離抽出した単一ニューロン活動波形とが、同時刻に発生した数を表す。

このようにして算出したノイズ率を図１１図に示す。 図１１は、スパイクソーティング装置１のノイズ率と、従来技術のノイズ率とを比較して示した図である。 この図に示すように、従来技術では、ノイズ率が２９．６〜４３．８％である。 一方、スパイクソーティング装置１では、ノイズ率が０．５〜７．１％である。 このことから、スパイクソーティング装置１では、従来技術に比べて、分離結果に含まれるノイズの数が優位に少ないことがわかる。 すなわち、スパイクソーティング装置１は、従来技術に比べて、より精度よく単一ニューロン活動波形を分離できる。

さらに、単一ニューロン活動波形の分離時間の例を、図１２に示す。 図１２は、スパイクソーティング装置１における単一ニューロン活動波形の分離時間と、従来技術における単一ニューロン活動波形の分離時間とを比較して示した図である。 図１２に示す結果は、神経細胞から１時間に渡って記録した生波形データを、それぞれ、スパイクソーティング装置１および従来技術において処理したときに要した時間である。 この図に示すように、従来技術では、単一ニューロン活動波形の分離に３９時間を要した。 一方、スパイクソーティング装置１では、２２時間の時間を要した。 したがって、スパイクソーティング装置１は、従来技術に比べて約１．８倍の速度で単一ニューロン活動波形を分離できたといえる。 このように、スパイクソーティング装置１における分離速度は、従来技術に比べて優位に早いことが実証された。

以上のように、スパイクソーティング装置１は、テトロードによって神経細胞から測定される複数の生波形から、単一ニューロン活動波形を分離する。 このとき、スパイクソーティング装置１は、生波形からスパイク波形を抽出し、抽出したスパイク波形から、期待値最大化法に基づき、クラスタ波形を生成する。 この期待値最大化法として、たとえば、ＫｌｕｓｔａＫｗｉｋのアルゴリズムを使用する。 また、スパイクソーティング装置１は、複数のクラスタ波形から、独立成分分析法に基づき、独立成分波形を生成する。 この独立成分分析法として、たとえば、ＦａｓｔＩＣＡのアルゴリズムを使用する。 さらに、互いに類似する独立成分波形を、所定の条件、たとえば、不応期が存在する場合などに、互いに連結する。

すなわち、スパイクソーティング装置１は、マルチニューロン活動波形から、単一ニューロン活動波形を分離抽出するための方法を利用している。 そのため、クラスター境界を、正規分布に厳密にしたがって決定するわけではない。 ニューロン活動の場合、距離の非常に近いクラスターが、同じ神経細胞から発生していると経験的に仮定されている。 そのため、スパイクソーティング装置１では、それらを同一クラスターとして連結する。

このように、スパイクソーティング装置１は、ニューロン活動の特性を考慮にいれてクラスタリングする。 これにより、ｋ−ｍｅａｎｓクラスタリングを利用している従来技術に比べて、単一ニューロン活動波形の分離精度をより向上できる。 また、スパイクソーティング装置１は、ＡＩＣ基準（Ａｋａｉｋｅ'ｓ Ｉｎｆｏｒｍａｔｉｏｎ Ｃｒｉｔｅｒｉａ；赤池の情報量基準）などを利用して、クラスタを生成する際の演算を無駄に繰り返さない。 そのため、ｋ−ｍｅａｎｓクラスタリングを利用している従来技術に比べて、単一ニューロン活動波形をより高速に分離できる。

なお、本発明は上述した実施形態に限定されるものではなく、請求項に示した範囲で種々の変更が可能である。 すなわち、請求項に示した範囲で適宜変更した技術的手段を組み合わせて得られる実施形態についても、本発明の技術的範囲に含まれる。

たとえば、スパイク波形抽出部１０は、平均２乗誤差以外の他の閾値を使用してもよい。 すなわち、抽出されるスパイク波形の特性に応じた、任意のパラメータを閾値として使用してもよい。

また、スパイク特徴抽出部１１は、正ピーク値および負ピーク値とは異なる、他のスパイク特徴値をスパイク波形から抽出してもよい。 たとえば、スパイク波形の幅、スパイク波形の面積、生波形データに対して主成分分析を行って得られる固有値を、スパイク特徴値として抽出できる。

また、スパイクソーティング装置１が処理する生波形データは、サルから記録した生波形データに限定されない。 テトロードなどの多点電極によって神経細胞から記録した生波形データは、全て、スパイクソーティング装置１によって処理できる。 たとえば、ラットやマウスなどの小動物、サルなどの大型動物も対象にしてよい。 さらに、倫理的問題が解決さえすれば、人間を対象にすることも可能である。

生波形データの記録に用いるテトロードのインピーダンスは、１００ｋΩ〜１．５ＭΩの範囲内にあることが好ましい。 タングステンワイヤの太さは、２０ミクロン以下であることが求められる。 タングステンワイヤを束ねたテトロードの太さは、１００ミクロン以下であることが求められる。 なお、独立成分分析は、電極間の波形分布の差を、情報としている。 そのため、同一の神経細胞から発生しているスパイク波形を計測しなければ、単一ニューロン活動波形を分離するために充分な情報を取得できない。 そこで、テトロードのサイズを、計測範囲が重なる１００ミクロン以下とすればよい。

上述した本実施形態では、４本のタングステンワイヤからなるテトロードを例に、説明した。 しかし、テトロードに含まれるタングステンワイヤの本数は、４本に限定されない。 上述した１００ミクロンの範囲内に収まるのであれば、任意の本数が可能である。 したがって、スパイクソーティング装置１では、４以上のチャンネルの生波形データを処理できる。

この観点から言えば、スパイクソーティング装置１は、ｎ個（ｎは２以上の整数）の記録電極を含むテトロードを用いて神経細胞から測定されたｎ個の第１，第２，・・・第ｎ生波形から、ｎ個の生波形に含まれている少なくとも１つの単一ニューロン活動波形を分離する装置である。

また、スパイク波形抽出部１０は、閾値算出手段として、第１，第２，・・・第ｎ生波形のそれぞれに、所定の演算を施すことによって、第１，第２，・・・第ｎ生波形のそれぞれに対応する所定の第１，第２，・・・第ｎ閾値を算出する。

さらに、スパイク波形抽出部１０は、閾値判定手段として、第１，第２，・・・第ｎ生波形に含まれる、ある同一時間における第１，第２，・・・第ｎ観測信号が、それぞれ、第１，第２，・・・第ｎ閾値を超えるか否かを判定する。

さらに、スパイク波形抽出部１０は、タイムスタンプ算出手段として、第１，第２，・・・第ｎ生波形の少なくともいずれかに対して、閾値超過判定手段による判定結果が真である場合に、閾値を超えると判定された観測信号が出現する時間を、スパイク出現時間として算出する。

さらに、スパイク波形抽出部１０は、スパイク波形抽出手段として、スパイク出現時間において出現する観測信号を含む連続した所定の長さの波形を、第１，第２，・・・第ｎ生波形から、それぞれ、ｎ個の第１，第２，・・・第ｎスパイク波形として抽出する。

さらに、スパイク波形抽出部１０は、スパイク波形群生成手段として、抽出されたｎ個の第１，第２，・・・第ｎスパイク波形と、これらのスパイク波形の抽出に参照されたスパイク出現時間とを関連づけることによって、スパイク波形群を生成する。

さらに、スパイク特徴抽出部１１は、スパイク特徴抽出手段として、同一のスパイク波形群に含まれる第１，第２，・・・第ｎスパイク波形に対して所定の演算を行うことによって、第１，第２，・・・第ｎスパイク波形の特徴をそれぞれ表すｎ個の第１，第２，・・・第ｎスパイク特徴を抽出する。

さらに、スパイク特徴抽出部１１は、スパイク特徴ベクトル生成手段として、同一のスパイク波形群に含まれるスパイク波形から抽出されたｎ個のスパイク特徴を互いに関連づけることによって、スパイク波形群に対応するスパイク特徴ベクトルを生成する。

さらに、クラスタ波形生成部１２は、スパイク波形群類似性判定手段として、スパイク特徴ベクトル同士を互いに比較することによって、あるスパイク波形群と、他のスパイク波形群とが類似しているか否かを判定する。

さらに、クラスタ波形生成部１２は、クラスタ生成手段として、スパイク波形類似性判定手段による判定結果に基づき、あるスパイク波形群に含まれる第１，第２，・・・第ｎスパイク波形を、それぞれ、このスパイク波形群に類似すると判定されたスパイク波形群に含まれる第１，第２，・・・第ｎスパイク波形と関連づけることによって、互いに類似関係にある複数のスパイク波形群からなるスパイク波形クラスタを生成する。

さらに、独立成分波形抽出部１３は、独立成分波形抽出手段として、ｎ個の第１，第２，・・・第ｎ連結スパイク波形を含むスパイク波形クラスタから、１つの分離行列、および、記録電極の本数と同一数であるｎ個の第１，第２，・・・第ｎ独立成分波形を抽出する。

さらに、独立成分波形抽出部１３は、逆行列算出手段として、分離行列から、ｎ個の逆行列を算出する。

さらに、独立成分波形抽出部１３は、基底ベクトル生成手段として、ｎ個の独立成分波形から、それぞれ、ｎ個の第１，第２，・・・第ｎ基底ベクトルを生成する。

さらに、再構成波形生成部１４は、再構成波形生成手段として、スパイク波形クラスタに含まれる第１，第２，・・・第ｎ連結スパイク波形の全てに対して、第１，第２，・・・第ｎ逆行列をそれぞれ適用することによって、ｎ×ｎ個の再構成波形を生成する。

さらに、再構成波形生成部１４は、第２閾値算出手段として、ｎ×ｎ個の再構成波形から、それぞれに対応するｎ×ｎ個の第２閾値を算出する。

さらに、再構成波形生成部１４は、第２閾値超過判定手段として、再構成スパイク波形のいずれかが、第２閾値を超えるか否かを判定する。

さらに、再構成波形生成部１４は、第２連結手段として、第２閾値を超えると判定される再構成スパイク波形を互いに連結し、さらに、基底ベクトルを関連づけることによって、連結再構成スパイクを生成する。

さらに、単一ニューロン活動波形抽出部１５は、独立成分クラスタ類似性判定手段として、独立成分クラスタの全ての組み合わせに対して、独立成分クラスタに関連づけられている基底ベクトルに基づき、独立成分クラスタ同士が類似するか否かを判定する。

さらに、単一ニューロン活動波形抽出部１５は、単一ニューロン活動波形算出手段として、互いに独立すると判定された独立成分クラスタ同士を関連づけることによって、単一ニューロン活動波形を抽出する。

なお、上述した各部材は、いずれも機能ブロックである。 したがって、これらの部材は、ＣＰＵなどの演算手段が、図示しない記憶部に格納されたスパイクソーティングプログラムを実行し、図示しない入出力回路などの周辺回路を制御することによって、実現される。

したがって、本発明の目的は、上述した機能を実現するソフトウェアであるスパイクソーティングプログラムのプログラムコード（実行形式プログラム、中間コードプログラム、ソースプログラム）をコンピュータによって読み取り可能に記録している記録媒体を、スパイクソーティング装置に供給し、スパイクソーティング装置に備えられるコンピュータ（またはＣＰＵやＭＰＵ、ＤＳＰ）が、記録媒体に記録されているプログラムコードを読み出し実行することによって、達成可能である。

この場合、記録媒体から読み出されたプログラムコード自体が、上述した機能を実現する。 そのため、そのプログラムコードを記録している記録媒体は、本発明を構成することになる。

一方で、上述した各部材は、上述したソフトウェアと同様の処理を行うハードウェアとして実現してもよい。 この場合、本発明の目的は、ハードウェアとしてのスパイクソーティング装置によって達成されることになる。

ここで、プログラムコードを読み出し実行する演算手段は、単体の構成であればよい。 または、スパイクソーティング装置内部のバスや各種の通信路を介して接続されている複数の演算手段が、プログラムコードを協同して実行する構成であってもよい。

演算手段によって直接的に実行可能なプログラムコードを、このプログラムコードを格納しているコンピュータ読み取り可能な記録媒体を通じて、スパイクソーティング装置に配布すればよい。 また、プログラムコードを、後述する解凍などの処理によってプログラムコードを生成可能なデータとして、当該データを格納しているコンピュータ読み取り可能な記録媒体に通じて、スパイクソーティング装置に配布してもよい。 あるいは、これらのプログラムコードまたはデータを、有線または無線の通信路を介してデータを伝送する通信ネットワークを通じて、スパイクソーティング装置に配布または送信してもよい。 いずれの手段によって配布または送信されても、プログラムコードは、スパイクソーティング装置に備えられる演算手段によって実行される。

このとき、特定のものに限定されない各種の通信ネットワークを通じて、プログラムコードまたはデータを伝送できる。 このような通信ネットワークの具体例を挙げると、インターネット、イントラネット、エキストラネット、ＬＡＮ、ＩＳＤＮ、ＶＡＮ、ＣＡＴＶ通信網、仮想専用網（Ｖｉｒｔｕａｌ Ｐｒｉｖａｔｅ Ｎｅｔｗｏｒｋ）、電話回線網、移動体通信網、衛星通信網等がある。 また、通信ネットワークを構成する伝送媒体（通信路）も、特に限定されない。 具体的には、ＩＥＥＥ１３９４規格による回線、ＵＳＢ回線、電力線、ケーブルＴＶ回線、電話線、およびＡＤＳＬ回線等の有線を、伝送媒体として利用できる。 さらに、ＩｒＤＡやリモコンに用いられている赤外線を利用した無線、Ｂｌｕｅｔｏｏｔｈ規格またはＩＥＥＥ８０２．１１無線規格に規定されている無線、ＨＤＲ、携帯電話網、衛星回線、地上波デジタル網等を利用した無線も、伝送媒体として利用できる。

なお、プログラムコードをスパイクソーティング装置に配布するための記録媒体は、プログラムコードの配布前には、取り外し可能になっていることが好ましい。 しかし、プログラムコードを配布した後には、スパイクソーティング装置装置から取り外し可能になっていてもよく、スパイクソーティング装置と一体化されて取り外し不可能になっていてもよい。

また、記録媒体は、プログラムコードが記録されてさえいれば、書き換え（書き込み）可能であってもよく、不可能であってもよい。 また、揮発性であってもよく、非揮発性であってもよい。 さらに、記録媒体へのプログラムコードの記録方法、および記録媒体の形状も、任意のものでよい。

このような条件を満たす記録媒体を例示すると、磁気テープやカセットテープなどのテープ、フロッピー（登録商標）ディスクやハードディスクなどの磁気ディスク、ＣＤ−ＲＯＭや光磁気ディスク（ＭＯ）、ミニディスク（ＭＤ）、デジタルビデオディスク（ＤＶＤ）などのディスクがある。 さらに、ＩＣカードや光カードのようなカード型メモリ、あるいは、マスクＲＯＭやＥＰＲＯＭ、ＥＥＰＲＯＭまたはフラッシュＲＯＭなどの半導体メモリも該当する。 さらに、ＣＰＵなどの演算手段内に形成されているメモリも該当する。

なお、プログラムコードを記録媒体から読み出して主記憶に格納するためのプログラムは、あらかじめ、スパイクソーティング装置内に、コンピュータによって実行可能に格納されている。 また、プログラムコードを通信ネットワークを通じてスパイクソーティング装置に配布する場合、通信ネットワークからプログラムコードをダウンロードするプログラムは、あらかじめ、スパイクソーティング装置内に、コンピュータによって実行可能に格納されている。

また、プログラムコードは、上述した各処理の全手段を演算手段へ指示するコードであればよい。 なお、コンピュータには、プログラムコードによる各処理の一部または全部を所定の手順で呼び出すことによって実行可能な基本プログラム（たとえば、オペレーティングシステムやライブラリなど）がすでに存在している場合がある。 この場合、プログラムコードにおける全手順の一部または全部を、この基本プログラムの呼び出しを演算手段へ指示するコードやポインタなどに置き換えたプログラムコードものを、スパイクソーティングプログラムのプログラムコードとしてもよい。

また、記録媒体に、実メモリにプログラムコードを配置した状態のように、スパイクソーティングプログラムを格納すればよい。 具体的には、演算手段が記録媒体にアクセスしてプログラムコードを実行できる形式によって、スパイクソーティングプログラムを記録媒体に格納すればよい。 または、実メモリにプログラムコードを配置する前であり、かつ、演算手段が常時アクセス可能なローカルな記録媒体（たとえばハードディスクなど）にインストールした後の格納形式によって、スパイクソーティングプログラムを記録媒体に格納してもよい。 あるいは、通信ネットワークや搬送可能な記録媒体などからローカルな記録媒体にインストールする前の格納形式によって、スパイクソーティングプログラムを記録媒体に格納してもよい。

スパイクソーティングプログラムは、コンパイルされた後のオブジェクトコードに限られない。 たとえば、スパイクソーティングプログラムは、ソースコードとして記録媒体に格納されていてもよい。 あるいは、インタプリトまたはコンパイルの途中において生成される中間コードとして、記録媒体に格納されていてもよい。 上述したいずれの場合であっても、記録媒体に格納されているプログラムコード（中間コード）は、演算手段が実行可能な形式に変換可能なものであればよい。

すなわち、プログラムコード（中間コード）は、所定の形式変換プログラムが、圧縮されたプログラムコードを解凍したり、符号化されたプログラムコードを復元したり、ソースコードをインタプリト、コンパイル、リンク、または、実メモリへ配置したりすることによって、あるいはこれらの処理を組み合わせて実行することによって、演算手段が実行可能な形式に変換されるものであればよい。 これにより、スパイクソーティングプログラムを記録媒体に格納する際の格納形式にかかわらず、同様の効果を得ることができる。

本発明は、脳神経活動だけで機械を操作する最新技術である脳−機械インターフェース（Ｂｒａｉｎ Ｍａｃｈｉｎｅ Ｉｎｔｅｒｆａｃｅ）作成の基礎技術となる。 したがって、脊椎損傷によって四肢が不自由な患者などを対象とした、医療分野や介護分野に広く応用できる。 また、脳の神経活動を解析するマルチニューロン活動記録の解析装置の製造にも応用できる。

本発明の一実施形態に係るスパイクソーティング装置の構成を示すブロック図である。

スパイクソーティング装置が処理した生波形、および、生波形から抽出したスパイク波形の一例を示す図である。

生波形から抽出されるスパイク波形の一例を示す図である。

スパイク波形から抽出されるスパイク特徴値の一例を示す図である。

スパイク波形、タイムスタンプ、およびスパイク特徴値ベクトルが関連づけられたスパイク波形群の一例を示す図である。

２つのスパイク波形データを連結することによって、１つの連結クラスタ波形群を生成する例を示す図である。

スパイクソーティング装置によって抽出した独立成分波形の一例を示す図である。

（ａ）は、再構成波形群の一例を示す図であり、（ｂ）は、再構成波形群から抽出した再構成スパイク波形の例を示す図である。

（ａ）は、補正前のタイムスタンプが関連づけられている再構成波形群の一例を示す図であり、（ｂ）は、補正後のタイムスタンプが関連づけられている再構成波形群の一例を示す図である。

スパイクソーティング装置における単一ニューロン活動波形の分離精度と、従来技術における単一ニューロン活動波形の分離精度とを比較して示した図である。

スパイクソーティング装置のノイズ率と、従来技術のノイズ率とを比較して示した図である。

スパイクソーティング装置における単一ニューロン活動波形の分離時間と、従来技術における単一ニューロン活動波形の分離時間とを比較して示した図である。

符号の説明

１ スパイクソーティング装置１０ スパイク波形抽出部（スパイク波形抽出手段）
１１ スパイク特徴抽出部（スパイク特徴抽出手段）
１２ クラスタ波形生成部（クラスタ波形生成手段）
１３ 独立成分波形抽出部（独立成分波形抽出手段）
１４ 再構成波形生成部（再構成波形生成手段）
１５ 単一ニューロン活動波形抽出部（単一ニューロン活動波形抽出手段）