JP2021531096A - 神経細胞信号を監視するためのシステム及び方法 - Google Patents

Abstract

本開示の一態様では、患者の脳状態における少なくとも一つの振動成分を推定する方法が提供される。この方法は、脳波記録（ＥＥＧ）データを受け取ることと、少なくとも一つの振動成分を含む状態空間モデルをＥＥＧデータに当て嵌めることと、状態空間モデルに基づいてＥＥＧデータの交差周波数結合の少なくとも一つの値を推定することと、交差周波数結合の少なくとも一つの値に基づいて報告を生成することと、表示部とメモリの少なくとも一方に報告を出力することと、を含む。【選択図】図２３

Description

＜関連出願の相互参照＞
本願は、「神経細胞の振動を監視するためのシステム及び方法」と題されて２０１８年７月１６日付けで出願された米国仮特許出願第６２／６９８，４３５号に基づき、その優先権を主張するとともに、その全文を参照により本願に援用する。

＜連邦政府による資金提供を受けた研究の記載＞
本発明は、米国立衛生研究所による助成金第Ｒ０１ＡＧ０５６０１５−０１号、第Ｒ０１ＡＧ０５４０８１−０１Ａ１号及び第ＧＭ１１８２６９−０１Ａ１号に基づく資金提供を受けてなされたものであり、連邦政府は本発明に対する特定の権利を有する。

神経細胞信号の振動活動は普遍的であり、広範な時間的且つ空間的スケールにわたる神経細胞集団の協調活動を反映する。神経細胞の振動は、注意力や認識力や感覚処理及び意識状態をはじめとする多くの脳機能要素において重要な媒介的役割を果たすと考えられている。しかるべく、神経認知障害や精神医学的障害ならびに意識変容状態も、脳振動の変容又は撹乱に関連するとされている。神経細胞の振動活動は、動物モデルとヒトモデルのいずれにおいても、個々の神経細胞膜電位や神経細胞回路、更には非侵襲的頭皮電磁界に及ぶスケールにて記録できる。したがって、神経細胞の振動により、異なるスケールやモデルにわたり観測結果をリンクするための貴重な機構的土台を得ることができる。

あいにくのところ、神経細胞信号は、他の多くの生理学的信号と比べると極めて微妙なものであり、その全て又は一部が多種多様な影響因子によって識別が困難になる可能性がある。そのため、たいていの場合、神経信号を取得且つ又は分析する前にロバストな信号処理且つ又は信号識別が先行する。

例えば、ノンパラメトリックスペクトル分析などの周波数領域法や線形帯域通過フィルタ処理に基づく時間領域法などが脳振動分析によく使用される。典型的な分析法では、特定の周波数帯域における振動（例：周波数８〜１２Ｈｚであるアルファ帯域での振動）におけるパワーの推定値を求める場合がある。大抵の場合、８〜１２Ｈｚのパワーは、周波数領域においてパワースペクトルを用いて計算されたり、又は帯域通過フィルタ処理した信号のパワー又は分散を推定することにより時間領域において計算されたりする。

この種の処理は、様々な推定を前提としており、かかる推定は誤差を伴うものもあれば、意識的なトレードオフ（妥協により両立させるための手段）としてなされるものもある。例えば、上記周波数の選択とフィルタ処理では本質的に「ノイズ」の除去を目的として一部の「信号」を差別する。但し、「ノイズ」を取り消す試みの一環として失われた「信号」の一部が想定以上に貴重である場合もある。更に、バンドパス周波数フィルタ処理などの信号処理手法の多くは容易に実施可能であるが、「信号」を最も良く捉えるための信号処理手法としては理想的と言えない場合もある。

したがって、より簡単に神経細胞信号の獲得と分析を図れるシステムと方法を備えていることが望ましい。

本開示は、脳監視を超越した用途において患者の監視や信号の復調などのために使用される、振動状態に応じて神経細胞信号を処理するためのシステム及び方法を提供する。本開示の非限定的な一態様では、状態空間モデルにおける線形振動表現の組み合わせとして神経細胞信号を捕獲するためのシステムと方法を提供する。これらのモデルのパラメータとしては、データ中に存在する振動を特定するのに適したモデルを選択し、それによって所望の神経細胞信号を抽出するために期待値最大化アルゴリズムその他の方法を使用できる。本開示は、これらのシステムと方法により適切な抽出を成し遂げることにより、本来対象となるべき振動と、発生する可能性のある二次的な生理学的信号とを区別又は分離できる（かかる区別又は分離なくしては、二次的な生理学的信号がそれ以降の分析に混同して困惑を招き得る）。このように、本開示のシステムと方法により、位相振幅相互相関や、振動位相、振動振幅、振動パワー、相対振動パワー、振動帯域幅、共振、振動線質係数など、監視用途に有用な振動特性を推定できる。

本開示の一態様では、患者の脳の状態を示す少なくとも一つの振動成分を推定する方法を提供する。かかる方法には、脳波記録（ＥＥＧ）データを受け取ることと、少なくとも一つの振動成分からなる状態空間モデルをＥＥＧデータに当て嵌めることと、状態空間モデルに基づきＥＥＧデータの交差周波数結合の少なくとも一つの値を推定することと、交差周波数結合の少なくとも一つの値に基づいて報告を生成することと、表示部とメモリの少なくとも一方に報告を出力することと、が含まれる。

上記方法において、交差周波数結合の少なくとも一つの値は、状態空間モデルに含まれる第１波成分と第２の波成分の間の位相振幅変調の推定値を含み得る。

上記方法において、第１波は徐波成分であり得、第２波はアルファ波成分であり得る。

上記方法において、位相振幅変調の推定値は、徐波成分の位相とアルファ波成分の振幅に基づいて計算できる。

上記方法において、交差周波数結合の少なくとも一つの値は、状態空間モデルに含まれる振動成分とＥＥＧデータの周波数の帯域との相関値を含み得る。

上記方法において、状態空間モデルの当て嵌めには、少なくとも一つの振動成分に含まれる各振動成分の調和振動数の当て嵌めを含み得る。

上記方法において、少なくとも一つの振動成分はアルファ成分を含み得る。そして、アルファ成分は中央周波数に対する事前分布に関連する。事前分布はフォンミーゼスの事前分布であり得る。

上記方法において、状態空間モデルの減衰係数は事前分布で制約し得る。事前分布はベータ分布であり得る。

本開示の別の態様によると、患者の脳の少なくとも一つの振動成分を推定するためのシステムを提供できる。このシステムは、脳波記録（ＥＥＧ）データを受け取るように構成された複数のセンサと、ＥＥＧデータを受け取るよう構成されたプロセッサと、表示部と、を備え、プロセッサは、少なくとも一つの振動成分を含む状態空間モデルをＥＥＧデータに当て嵌めるステップと、状態空間モデルに基づいてＥＥＧデータの交差周波数結合の少なくとも一つの値を推定するステップと、交差周波数結合の少なくとも一つの値に基づいて報告を生成するステップと、を実行するように構成され、表示部は報告を表示するよう構成される。

上記システムにおいて、交差周波数結合の少なくとも一つの値は、状態空間モデルに含まれる第１波成分と第２波成分の間の位相振幅変調の推定値を含み得る。

上記システムにおいて、第１波は徐波成分であり得、第２波はアルファ波成分であり得る。

上記システムにおいて、位相振幅変調の推定値は、徐波成分の位相及びアルファ波成分の振幅に基づいて計算できる。

上記システムにおいて、交差周波数結合の少なくとも一つの値は、状態空間モデルに含まれる振動成分とＥＥＧデータの周波数の帯域との相関値を含み得る。

上記システムにおいて、状態空間モデルの当て嵌めには、少なくとも一つの振動成分に含まれる各振動成分の調和振動数の当て嵌めを含み得る。

上記システムにおいて、少なくとも一つの振動成分は、アルファ成分を含みえる。そして、アルファ成分は中央周波数に対する事前分布に関連する。事前分布はフォンミーゼスの事前分布であり得る。

上記システムにおいて、状態空間モデルの減衰係数は事前分布で制約し得る。事前分布はベータ分布であり得る。

本発明の前述及びその他の利点は以下の説明から明らかになるであろう。以下、本発明の好ましい実施形態を例示するために本願の一部として添付される図面を参照しながら本発明をより詳細に説明する。但し、かかる実施形態は、必ずしも本発明の範囲を全て表すものではなく、本願請求項の範囲及び本明細書に本発明の範囲の解釈を目的として言及されているに過ぎない。

生理学的監視システム１０の一実施形態を示す。 ＥＥＧセンサの代表例を示す。 患者監視用システムの代表例を示す。 データを取得して報告を作成するステップの代表例を示す。 データを取得して報告を作成するステップの別の代表例を示す。 プロポフォール投与中の患者における自己回帰モデルのＥＥＧスペクトルを示す。 休息状態における自己回帰モデルのＥＥＧスペクトルを示す。 複合振動時系列を示す。 図６Ａの徐波振動時系列を示す。 図６Ａのアルファ振動時系列を示す。 バンドパスの徐波成分とアルファ成分に関するマルチテーパースペクトルを示す。 バンドパス法と状態空間法における総信号を示す。 バンドパス法と状態空間法における徐波振動成分を示す。 バンドパス法と状態空間法におけるアルファ振動成分を示す。 バンドパス法と状態空間法における残差を示す。 ＥＥＧデータのＡＲ成分のスペクトルを示す。 ＥＥＧデータの推定振動スペクトルを示す。 フォンミーゼスの事前分布を伴わない徐波・アルファモデルにおけるＡＲ成分のスペクトルを示す。 フォンミーゼスの事前分布を伴わない徐波・アルファモデルにおける推定振動スペクトルを示す。 アルファ無し徐波モデルにおけるＡＲ成分のスペクトルを示す。 アルファ無し徐波モデルにおける推定振動スペクトルを示す。 状態空間モデルによりＥＥＧデータを当て嵌めるステップの代表例１１００を示す。 （Ａ）はＥＥＧの生データを示す。（Ｂ）はＥＥＧデータの６秒の時間窓を示す。（Ｃ）は６秒の時間窓に当て嵌まる徐波振動を示す。（Ｄ）は６秒の時間窓に当て嵌まる速波振動を示す。 （Ｅ）は状態空間モデルの代表形式を示す。（Ｆ）は確信区間での徐波振動位相を示す。（Ｇ）は確信区間での速波振動振幅を示す。（Ｈ）は推定されるＫ^ｍｏｄ及びそれに関連する分布を示す。（Ｉ）は推定されるφ^ｍｏｄ及びそれに関連する分布を示す。（Ｊ）は回帰推定されたアルファ波振幅を示す。 所与のデータに関するパラメトリックパワースペクトル密度データの第１のパス当て嵌めを示す。 該パラメトリックパワースペクトル密度データの閾値を示す。 パラメトリックパワースペクトル密度データの第２のパス当て嵌めを示す。 回帰推定されたパワースペクトル密度データを示す。 パラメトリックパワースペクトル密度データの当て嵌め線を示す。 患者のＥＥＧデータに供給されるプロポフォール濃度を示す。 該患者のＥＥＧデータに対応する反応確率を示す。 該患者のＥＥＧデータの変調回帰値を示す。 該患者のＥＥＧデータのパラメトリックスペクトルを示す。 該患者のＥＥＧデータのモジュログラム（変調グラフ）ＰＡＭを示す。 該患者のＥＥＧデータの変調パラメータを示す。 反応確率を示すプロット図である。 はプロポフォール濃度を示すプロット図である。 ６秒の時間窓を使った標準的な方式の場合のモジュログラムを示す。 図１５Ｃに関連する変調指数を示すプロット図である。 １２０秒の時間窓を使った標準的な方式の場合のモジュログラムを示す。 図１５Ｅに関連する変調指数を示すプロット図である。 ＳＳＣＦＡ法の場合のモジュログラムを示す。 図１５Ｇに関連する変調指数を示すプロット図である。 ｄＳＳＣＦＡ法の場合のモジュログラムを示す。 図１５Ｉに関連する変調指数を示すプロット図である。 反応確率を示すプロット図である。 プロポフォール濃度を示すプロット図である。 ６秒の時間窓を使った標準的な方式の場合のモジュログラムを示す。 図１６Ｃに関連する変調指数を示すプロット図である。 １２０秒の時間窓を使った標準的な方式の場合のモジュログラムを示す。 図１６Ｅに関連する変調指数を示すプロット図である。 ＳＳＣＦＡ法の場合のモジュログラムを示す。 図１６Ｇに関連する変調指数を示すプロット図である。 ｄＳＳＣＦＡ法の場合のモジュログラムを示す。 図１６Ｉに関連する変調指数を示すプロット図である。 （Ａ）は従来の方法による振動成分の分解を示す。（Ｂ）は（Ａ）の従来の方法を用いて求められたパワースペクトル密度を示す。（Ｃ）は（Ａ）の従来の方法による位相振幅相互相関の推定を示す。（Ｄ）はＳＳＣＦＡ方法による振動成分の分解を示す。（Ｅ）は（Ｄ）のＳＳＣＦＡ法を用いて求められたパワースペクトル密度を示す。（Ｆ）は（Ｆ）のＳＳＣＦＡ法による位相振幅相互相関の推定を示す。 振幅変調の動的特性を示す。 位相変調の動的特性を示す。 正面電極に関する相関関係ベースの交差周波数結合を示す。 モジュログラムを組み合せてまとめた図を示す。 交差周波数結合の、周波数関数としての主要周波数モードを示す。 該主要周波数モードの総エネルギーのパーセンテージを示す。 第１主要周波数モードについて推定される交差周波数結合のパターンを示す。 第１主要周波数モードについて推定される異なる状態における脳の異なる個所での低周波活動による広帯域周波数結合を示すグラフである。 当て嵌め処理がなされた状態空間モデルを用いて交差周波数結合分析を行うステップの代表例を示す。

以下、患者監視装置及びその使用方法の代表例について説明する。該システムは、非限定的な一例として、後述の如く、患者の生理学的監視（少なくとも一種の麻酔性薬物投与の影響下にある患者を監視するなど）を行うために利用できるセンサを含む。

例えば、後述の如く、該システムと方法の非限定的な一用途としては全身麻酔、鎮静状態又はその他の薬物誘発性状態の監視が挙げられる。したがって、本開示の一態様では、少なくとも一種の麻酔性薬物投与の影響下にある患者を監視するシステムを提供する。該システムは、少なくとも一種の麻酔性薬物投与中に患者から生理学的データを収集するように構成された複数のセンサと、該複数のセンサから生理学的データを受け取り、該患者が全身麻酔下又は鎮静状態にある時に脳の状態を判定するために使用できる信号マーカを該生理学的データに基づいて決定するように構成されたプロセッサと、を含み得る。

また、該システムと方法は、全身麻酔状態、鎮静状態又はその他薬物誘発性の状態に要される脳状態を作り出し規制するための薬物投与の閉ループ制御を提供し得る。したがって、本開示の一態様では、該少なくとも一種の麻酔性薬物投与の影響下にある患者の閉ループ制御を行うためのシステムを提供する。該システムは、少なくとも一種の麻酔性薬物投与中に患者から生理学的データを収集するように構成された複数のセンサと、該複数のセンサから生理学的データを受け取り、該患者が全身麻酔下又は鎮静状態にある時に脳の状態を判定するために使用できる信号マーカを該生理学的データに基づいて推測するように構成されたプロセッサと、を含み得る。また、該システムは、監視装置により提供される脳状態情報に基づいて薬物を投与するためのドラッグデリバリーシステム（薬物伝達システム）を含み得る。

図面を参照すると、図１Ａは、生理学的監視システム１０の一例を示す。生理学的監視システム１０では、一つ又は複数のセンサ１３を使って医療患者１２を監視している。該センサ１３はそれぞれ、ケーブル１５又はその他の通信リンクや媒体を通じて生理学的監視装置１７に信号を送信する。生理学的監視装置１７は、プロセッサ１９と、任意選択的に表示部１１を含む。一つ又は複数のセンサ１３は、例えば電気的なＥＥＧセンサなどの検出素子を含む。各センサ１３は患者１２の生理的パラメータを測定することによりめいめいの信号を生成できる。その後、生成された信号は、一つ又は複数のプロセッサ１９で処理される。この際、表示部１１が備えられている場合、一つ又は複数のプロセッサ１９は、処理した信号を表示部１１に伝達する。一実施形態において、表示部１１が生理学的監視装置１７に組み込まれているが、別の実施形態においては、表示部１１は生理学的監視装置１７と別々に備えられている。一例として、監視システム１０を携帯可能監視システムとして構成できる。また、別の例として、監視システム１０を、表示部無しのポッド状にして、別途設置される表示部に生理的パラメータデータを提供するように構成することもできる。

説明を明確にする目的で、図１Ａにおいては一つ又は複数のセンサ１３が一つのブロックとして示されているが、当然のことながら、図中のセンサ１３は一つ又は複数のセンサを示すものとする。一実施形態において、該一つ又は複数のセンサ１３は、後述されるいずれかの種類のセンサを一つ含む。他の実施形態において一つ又は複数のセンサ１３は、少なくとも２つのＥＥＧセンサを含む。更に別の実施形態において、該一つ又は複数のセンサ１３は少なくとも２つのＥＥＧセンサと一つ又は複数の脳酸素供給センサなどを含む。前述の各実施形態においては、任意選択的に異なる種類のセンサを追加することもできる。また、生理学的監視システム１０で使用するセンサの数と種類の組み合わせは、上記以外の組み合わせでもよい。

図１Ａに示すシステムの構成例では、信号をセンサから受け取って処理するために使用される全てのハードウェアを同じ筺体に収納している。他の実施形態では、信号をセンサから受け取って処理するために使用されるハードウェアの一部を別の筺体に収納している。更に、特定の実施形態において、生理学的監視装置１７は、センサ１３により送信される信号を受け取って処理するために使用されるハードウェア、ソフトウェア、又はハードウェアとソフトウェアの両方を、一つの筺体又は複数の筺体に収納した状態で含む。

図１Ｂに示されるように、ＥＥＧセンサ１３はケーブル２５を含み得る。ケーブル２５は、電気シールド内に３つの伝導体を含み得る。一つの伝導体２６は、生理学的監視装置１７に電力を供給できる。一つの伝導体２８は、生理学的監視装置１７に接地信号を供給できる。また、一つの伝導体２８は、センサ１３から生理学的監視装置１７に信号を伝送できる。センサが複数の場合は、ケーブル１５を更に一本又は複数本追加できる。

一部の実施形態では、接地信号をアースに接地させるが、他の実施形態では、接地信号を患者に接地させるようになっている。また、実施形態の中には、ケーブル２５が電気シールド層内に２つの伝導体を担持しており、該シールド層が接地線として機能するようになっているものもある。ケーブル２５内の電気的インターフェース２３が、生理学的監視装置１７のコネクタ２０内の電気的インターフェース２１に電気的に接続されることにより、ケーブル２５と生理学的監視装置１７が電気的に接続される。他の実施形態においては、センサ１３と生理学的監視装置１７がワイヤレス方式で通信する。

一部の構成例において、図１Ａと図１Ｂに示すシステムは、基準情報又はその他のデータを含むために、プロセッサ１９によりアクセス可能なメモリやデータベース又はその他のデータ記憶場所（図示せず）を更に含み得る。特に、かかる基準情報には、様々な年齢カテゴリなどの患者カテゴリの一覧を、関連する信号や信号マーカ又はシグネチャと共に含み得る。例えば、信号マーカ又はシグネチャには、様々な信号振幅、位相、頻度、パワースペクトル、又はそれの範囲や組み合わせ、その他信号マーカ又はシグネチャを含み得る。一部の態様において、プロセッサ１９は、かかる基準情報を用いて特定の患者特性や、患者の外見上の年齢、又は推定年齢又はその他患者の状態やカテゴリを決定できる。その他の態様においては、選択且つ又は決定された患者の特徴に基づいてデータ取得を規制又は変調し得る。

特に、ここで図２を参照すると、本開示の態様に則るシステムの代表例２００が示されているが、このシステム２００は、独立型の脳監視装置としても携帯可能装置としても構成することができ、若しくは既存の脳監視装置の中心的要素として組み込むこともできる。後述の説明からも明らかなように、システム２００は、様々な医療処置の実施を伴う手術室又は集中治療の場面（例えば、麻酔薬の投与時、ならびに術前・術後の評価時など）において貴重な役割を果たし得る。

システム２００は、図２に脳波記録（ＥＥＧ）電極アレイとして示される患者監視装置２０２（生理学的監視装置など）を含む。但し、患者監視装置２０２は数々の異なるセンサを含むものとして企図されている。特に、患者監視装置２０２は、例えば、外部刺激による覚醒を測定するために電気皮膚反応（ＧＳＲ）を監視する機構や、その他心電図や血圧監視装置をはじめとする心臓血管系監視装置ならびに固視微動監視装置などの監視システムも含み得る。この設計の一具現化例では、ＧＳＲ且つ又は固視微動を測定するための追加の電極を備えた前頭ラプラシアンＥＥＧ電極レイアウトを利用することができる。この設計の別の具現化例では、前述のＥＥＧシグネチャを最適に検出できることが判明している電極と、ここでも別個のＧＳＲ電極との任意の組み合わせを得るために後処理において組み合わせられ得る、前頭配列の電極を組み込むことができる。この設計の更に別の具現化例は、これもまた別個のＧＳＲ電極を備えた、音源定位のための６４〜２５６個の間のセンサを用いて頭皮表面全体をサンプリングする高密度のレイアウトを利用することができる。

患者監視装置２０２は、ケーブル２０４を介して監視システム２０６と接続されており、かくして監視システム２０６との通信を行う。また、構成部材間のケーブル２０４やそれに類似する接続を、無線接続に置き換えることも可能である。監視システム２０６は、患者監視装置２０２からの生信号（ＥＥＧ電極アレイから得た信号など）を受け取って、該信号を組み立てて処理し、更には様々な形式（時系列波形、スペクトログラムなど）で表示するように構成し得る。一部の作動形態において、監視システム２０６は、患者監視装置２０２のセンサからスカウトデータを生理学データとして又はその他のデータ形式で取得し、それに含まれる信号マーカやシグネチャを特定するように設計し得る。例えば、信号振幅や位相、周波数、パワースペクトル、その他の信号マーカやシグネチャを、様々適した方法を使って、スカウトデータやその他取得したデータにおいて特定し得る。加えて、マルチテーパー分析を行って、数桁におよぶ信号のダイナミックレンジを特定・考慮することもできる。その後、監視システム２０６は、かかる信号マーカ又はシグネチャを使用して、患者の外見上の年齢又は推定年齢などをはじめとする様々な患者特性を決定し得る。

一実施形態において、監視システム２０６による生理学的データの取得は、スカウトデータから決定された患者特性に基づいて調節又は規制し得る。特に、特定の決定された患者特性に応じてスケールを決定し、それ以降のデータ取得を、該決定されたスケール且つ又はユーザにより提供され如何なる指示に基づいて調節するように、監視システム２０６を構成し得る。例えば、一つ又は複数の増幅部利得を、その他のデータ取得パラメータとともに調節することにより、データの取得を規制し得る。なおまた、一部の態様において、取得した様々な生理学的データを、該スケールで表示できる特定の形式にフォーマットするように監視システム２０６を更に構成し得る。このように、年齢に合わせたスケールを、患者の外見上の年齢又は推定年齢に基づいて決定しえるので、それ以降に該選択された年齢適合スケールを使ってデータを収集でき、よって、年齢をもとに補正されたデータを生成・例示することができる。

図にも示される如く、監視システム２０６を更に専用の分析システム２０８に接続し得る。但し、監視システム２０６と分析システム２０８を共通システムに統合する又は組み合わせることも可能である。分析システム２０８は、監視システム２０６からＥＥＧ波形を受け取り、後述の如くＥＥＧ波形とシグネチャを分析する。但し、監視システム２０６と分析システム２０８の如何なる分析又は処理機能を、必要又は所望に応じて共有させることも個別に分散させることも企図される。

一部の態様においては、特定の医療処置を受ける患者に関して決定された特性に関連する情報を、システム２００の操作者又は臨床医に提供し得る。例えば、以前から、高齢の患者は手術室でバーストサプレッションを呈する可能性が高いことが判明している。特に、バーストサプレッションは深刻な脳不活性状態のことであり、電気的活動の群発が一定期間で抑制されて等電期間と交互に出現するものをいう。アルファ波（８〜１０Ｈｚ）と徐波（０．１〜４Ｈｚ）の信号振動を特徴とする麻酔誘発性の無意識な脳状態は、バーストサプレッションを呈するのに要される投与量よりも少ない投与量の麻酔剤で達成できる。これは、麻酔剤の投与量を、高齢患者に対して現在推奨されている投与量よりもかなり少ない投与量に減らすべきであることを意味し得る。一般的に、現在推奨されている投与量では、高齢患者がバーストサプレッションを呈することから、リアルタイムでのＥＥＧ監視に基づいて麻酔量を滴定することにより、十分な全身麻酔状態を得るとともに麻酔ガス被曝を削減することが可能になり得る。よって、システム２００は、適切な麻酔量を選択する際に使われる情報を、定められた患者特性に基づいて提供し得る。このように、例えば、全身麻酔下の高齢患者が手術後に認知障害を呈する確率を低減し得る。

他の例では、監視システム２０６且つ又は分析システム２０８は、若年患者・中年患者・高齢患者ならびに薬物中毒の患者などの特定の患者に応じた術前又は術後評価を行い、患者に関する特定の状態（麻酔薬感受性や、認知障害などの術後に発生し得る如何なる合併症）を特定且つ又は予測するために利用できる先験的情報を見極める機能も有する。なおまた、麻酔ケアや麻酔後ケア、又は集中治療のための特有のレジメンなどの特定の指示を与え得る。

また、システム２００は、ドラッグデリバリーシステム２１０を含み得る。ドラッグデリバリーシステム２１０を分析システム２０８と監視システム２０８に連結して、該システム２００が閉ループ監視・制御システムをなすようにし得る。このような本発明に係る閉ループ監視・制御システムは、多岐にわたる運用ができるが、ユーザーインタフェース２１２を含むことにより、ユーザが閉ループ監視・制御システムのコンフィギュレーションを設定したり、閉ループ監視・制御システムからフィードバックを受けたり、また、必要に応じて、閉ループ監視・制御システムのコンフィギュレーションを変更したり且つ又はオーバーライドしたりできるようになる。一部の構成例において、ドラッグデリバリーシステム２１０は、患者の意識状態を麻酔により低下させる目的（全身麻酔又は鎮静の目的）でなされる麻酔剤の投与を制御する機能を有するだけではなく、患者の意識を増減させるためのシステムと方法を実施且つ反映することもできる。

例えば、本発明の一態様において、ドーパミンやノルエピネフリン再摂取の輸送体の抑制剤としてメチルフェニデート（ＭＰＨ）を使用して、イソフルランによる全身麻酔からの意識回復を積極的に誘発することもできる。メチルフェニデートは、意識を回復させたり、覚醒に似た脳波変動を誘発したり、呼吸ドライブを高めたりする際に使用できる。メチルフェニデート誘発性の行動的効果と呼吸効果はドロペリドールにより抑制できる。これは、メチルフェニデートがドーパミン作用性覚醒経路の活性化により覚醒を誘発するという所見を支持するものである。プレチスモグラフィー（容積脈波記録法）と血液ガス実験では、メチルフェニデートは毎分換気量を増加させ、よって、脳から麻酔が排除される速度を速める、ということが立証されている。また、上記のような制御システムを利用して、エチルフェニデート又はその他の薬剤によって覚醒を促進することにより、イソフルランやプロポフォール又はその他麻酔剤による全身麻酔からの意識回復覚醒を誘発することもできる。

したがって、図２を参照して説明したようなシステムは、ドラッグデリバリーシステム２１０を２つの特定のサブシステムとともに含むことにより、麻酔からの意識回復を誘発するために備えることができる。このように、ドラッグデリバリーシステム２１０は、一種又は複数の麻酔剤を被験者に配送するべく設計された麻酔剤投与システム２２４と、全身麻酔効果を覆したり又は被験者の麻酔からの自然な意識回復を促進したりする一種又は複数の薬剤を被験者に配送するべく設計された覚醒剤投与システム２２６とを含み得る。

例えば、メチルフェニデート及びそれの類似物と派生物は、覚醒と呼吸ドライブを高めることにより、患者の麻酔誘発性無意識状態からの意識回復を誘発する。したがって、手術の終わりに全身麻酔による無意識状態と呼吸抑制作用をくつがえすためにメチルフェニデート、アンフェタミン、モダフィニル、アマンタジン又はカフェインを配送する際に覚醒剤投与システム３２６を使用できる。メチルフェニデートとしては、デクストメチルフェニデート（Ｄ−ＭＰＨ）、ラセミ体メチルフェニデート又はレバメチルフェニデート（Ｌ−ＭＰＨ）、若しくはこれらが同じ又は異なる割合（約５０％：５０％、又は約６０％：４０％、又は約７０％：３０％、又は８０％：２０％、９０％：１０％、９５％：５％など）で含まれる組成物を使用し得る。メチルフェニデートの投与量を注意力欠如障害（ＡＤＤ）又は注意欠陥多動性障害（ＡＤＨＤ）の治療に用いられる投与量よりも多くして、その他の薬剤を投与し得る。かかるメチルフェニデートの投与量としては、約１０ｍｇ／ｋｇ〜約５ｍｇ／ｋｇであり得、また、約５ｍｇ／ｋｇ〜１０ｍｇ／ｋｇの如何なる整数値であり得る。場合によっては、かかる投与量は、約７ｍｇ／ｋｇ〜約０．１ｍｇ／ｋｇ、又は約５ｍｇ／ｋｇ〜約０．５ｍｇ／ｋｇであり得る。その他の薬剤としては、吸入投与可能なものも含み得る。

図３に移ると、本開示の態様に係るステップ３００が示されている。まず、ステップブロック３０２では、如何なる量の生理学的データを取得するが、ここにおける生理学的データは、ＥＥＧ信号などの生理学的信号を表し、例えば、患者監視装置２０２を利用して患者から得られるデータである。一部の態様において、該生理学的データには、様々な患者特性の決定などの目的でスカウトデータを含み得る。その後、ステップブロック３０４では、取得した生理学的データを用いて信号マーカやシグネチャを特定又は決定する。例えば、信号振幅や位相、周波数、パワースペクトル、その他の信号マーカやシグネチャを、様々適した方法を使って、スカウトデータやその他取得したデータにおいて特定し得る。

一部の好ましい実施形態において、かかる信号マーカ又はシグネチャは、患者の外見上の年齢又は推定年齢などをはじめとする様々な患者特性を決定するために使用し得る。加えて、ステップブロック３０４には、決定された患者特性に応じてスケールを決定するステップも含み得る。一態様においては、数桁におよぶ広範な信号のダイナミックレンジを考慮できるスペクトル推定法（マルチテーパー法など）を採用できる。別の態様においては、信号振幅の自動見積りを行って、可視化スケールならびに取得増幅部利得に適した年齢コホートと主治医設定を推測し得る。

次のステップブロック３０６では、スカウトデータに基づいて決定された信号マーカ又はシグネチャを用いて、それ以降に取得される信号データに関して、データ取得ステップの調節又は規制を行い得る。例えば、その他のデータ取得パラメータとともに一つ又は複数の増幅部利得を調節することによりデータ取得を規制し得る。一部の態様において、データ取得の規制は、決定された患者特性に応じて決定されたスケールを使うことと、該決定されたスケール且つ又はユーザにより提供される如何なる指示に基づいて以後のデータ取得を調節することを含み得る。一例として、ステップブロック３０４にて患者の外見上の年齢又は推定年齢に基づいて決定された年齢に合わせたスケールを使うことにより、それ以降、該選択された年齢適合スケールを使って如何なるデータを取得できるので、年齢に基づいて補正されたデータが生成されることになる。

ステップブロック３０８では、上記のように取得したデータを使って、患者の現在又は将来的な脳状態を見極めることが可能になる。例えば、該年齢をもとに補正されたデータを用いて組み立てた分析又は処理済のＥＥＧ波形を利用して現在且つ又は将来的な麻酔又は鎮静作用の深さを評価し得る。また、かかる脳状態の決定には、臨床医又はユーザにより提供される如何なる情報（医療処置に関連する情報など）を含み得る。

その後、ステップブロック３１０では、例えば、印刷された報告書又は好ましくは即時表示などの形で報告が生成される。かかる報告には、生の又は処理済データ、シグネチャ情報、現在且つ又は将来的な脳状態の表示、ならびに、患者の推定年齢且つ又は外見上の年齢などの患者特有の特徴に関連した情報を含み得る。シグネチャ情報又は見極めた状態は、波形やスペクト対数ラム、コヒ対数ラム、確率曲線などの形で表示し得る。一部の態様において、該報告には、スケールに対して表示される形式にフォーマットされた生理学的データを含み得る。その他の態様において、該報告には、麻酔薬感受性や、認知障害などの術後に発生し得る如何なる合併症、更には麻酔ケア、麻酔後ケア、又は集中治療のための特有のレジメンなども示し得る。

図４に移ると、本開示の態様に係る別のステップ４００における各ステップが示されている。具体的に説明すると、ステップ４００はステップブロック４０２から開始し、ステップブロック４０２では、サンプル又はスカウトデータが、たとえば前述のような患者監視システムを用いて取得される。次いで、ステップブロック４０４において、サンプルデータは、様々な調節又は対照分類を用いて分析されて、取得したサンプルデータを表す患者カテゴリを特定する。具体的に説明すると、このステップは、サンプルデータ内の信号マーカ又はシグネチャを特定することと、対照分類に関連付けられた信号マーカ又はシグネチャとの比較を実行することとを含む。例えば、信号振幅、位相、周波数、パワースペクトル、及び他の信号マーカ又はシグネチャを、様々な適切な方法を用いてサンプルデータ内で検出し得る。

ステップブロック４０４において実行された分析は、患者の外見上の年齢且つ又は推定年齢など特有の患者特性を示し得る。一部の態様では、特有の患者特性を示す特定された、又は明白なカテゴリを、ステップブロック４０６において任意選択的に表示し得る。更に、ステップブロック４０８では、ユーザ入力も受け取り得る。

その後、ステップブロック４１０において、様々な通信パラメータに関して決定がなされる。これは、決定された又は推測された患者特性又はカテゴリ、及び任意選択的にユーザ入力を考慮に入れることを含む。例えば、ステップブロック４１０においては、取得したデータの年齢適合スケールが、決定された患者特性且つ又は信号、取得したデータ内に存在する信号マーカ又はシグネチャに基づいて決定し得る。次いで、ステップブロック４１２において、その後のデータ取得を、決定された通信パラメータを使用して、年齢に応じたデータを取得できるように調節し得る。前述の如く、データ取得を調節することは、通信パラメータを用いて様々な増幅部利得を適切に調節又は変更することを含みえる。一部の態様では、決定された通信パラメータを、取得したサンプルデータに直接的に適用し得る。例えば、年齢に応じたスケールが、サンプルデータに適用されて年齢に応じたデータを作り出すことができる。

次いで、ステップブロック４１４において、前述の方法で取得又は処理したデータは、患者の現在の又は将来的な脳状態を決定するために使用し得る。例えば、年齢補償されたデータを用いて組み立てられた分析又は処理済みのＥＥＧ波形は、麻酔又は鎮静の現在且つ又は将来的な深度を査定するために使用し得る。更に、かかる脳状態を決定することには、医療処置に関係する情報などの、臨床医又はユーザによって提供された任意の情報を含み得る。

次いで、ステップブロック４１６においては、任意の適した形状又は形で報告が生成される。一部の態様では、該報告が、スケーリング（縮拡）されたデータ又は該データを表すデータカテゴリを表示しえる。他の態様では、該報告に、麻酔薬感受性や、術中又は術後に合併症が発生する確率、患者の外見上の年齢又は推定年齢、及び本開示の態様に関係する他の情報を示し得る。

患者の監視に有用な指標又はその他の情報を導出するべく、神経細胞信号を評価するために利用できる脳状態の指標を提供する目的で本システムと方法を使用できる数々の様式が存在する。かかる非限定的な一例としては次のものが挙げられる。

徐波振動（０．１〜１Ｈｚ）、デルタ振動（１〜４Ｈｚ）、シータ振動（４〜８Ｈｚ）、アルファ振動（８〜１２Ｈｚ）、ベータ振動（１２〜２５Ｈｚ）、ガンマ振動（２５〜７０Ｈｚ）、及び高ガンマ振動（７０Ｈｚ以上）などの振幅特性の決定且つ又は分析

徐波振動（０．１〜１Ｈｚ）、デルタ振動（１〜４Ｈｚ）、シータ振動（４〜８Ｈｚ）、アルファ振動（８〜１２Ｈｚ）、ベータ振動（１２〜２５Ｈｚ）、ガンマ振動（２５〜７０Ｈｚ）、及び高ガンマ振動（７０Ｈｚ以上）などにおける分析能特性の決定且つ又は分析

徐波振動（０．１〜１Ｈｚ）、デルタ振動（１〜４Ｈｚ）、シータ振動（４〜８Ｈｚ）、アルファ振動（８〜１２Ｈｚ）、ベータ振動（１２〜２５Ｈｚ）、ガンマ振動（２５〜７０Ｈｚ）、及び高ガンマ振動（７０Ｈｚ以上）をはじめとする振動におけるピーク周波数の特性の決定且つ又は分析

徐波振動（０．１〜１Ｈｚ）、デルタ振動（１〜４Ｈｚ）、シータ振動（４〜８Ｈｚ）、アルファ振動（８〜１２Ｈｚ）、ベータ振動（１２〜２５Ｈｚ）、ガンマ振動（２５〜７０Ｈｚ）、及び高ガンマ振動（７０Ｈｚ以上）などに関するセンサ間の又は複数のセンサ間のコヒーレンス特性の決定且つ又は分析

徐波振動（０．１〜１Ｈｚ）、デルタ振動（１〜４Ｈｚ）、シータ振動（４〜８Ｈｚ）、アルファ振動（８〜１２Ｈｚ）、ベータ振動（１２〜２５Ｈｚ）、ガンマ振動（２５〜７０Ｈｚ）、及び高ガンマ振動（７０Ｈｚ以上）などの位相特性の決定且つ又は分析

徐波振動（０．１〜１Ｈｚ）、デルタ振動（１〜４Ｈｚ）、シータ振動（４〜８Ｈｚ）、アルファ振動（８〜１２Ｈｚ）、ベータ振動（１２〜２５Ｈｚ）、ガンマ振動（２５〜７０Ｈｚ）、及び高ガンマ振動（７０Ｈｚ以上）に関するセンサ間の又はセンサ郡間の位相差特性の決定且つ又は分析

徐波振動（０．１〜１Ｈｚ）、デルタ振動（１〜４Ｈｚ）、シータ振動（４〜８Ｈｚ）、アルファ振動（８〜１２Ｈｚ）、ベータ振動（１２〜２５ Ｈｚ）、ガンマ振動（２５〜７０ Ｈｚ）、及び高ガンマ振動（７０ Ｈｚ以上）などを含む複数の振動同士の交差周波数結合関係についての特徴の決定且つ又は分析（例：ある一つのオシレータの位相と別のオシレータの振幅との間の関係）

一例では、全身麻酔又は鎮静状態を監視・追跡するために、シータ帯域（４〜８Ｈｚ）とアルファ帯域（８〜１２Ｈｚ）及びベータ帯域（１２〜２５Ｈｚ）に及ぶ振動のピーク周波数を提示すことができる。また、これに伴い、徐波振動（０．１〜１Ｈｚ）及びデルタ（１〜４Ｈｚ）帯域におけるパワーを提示すことができる。例えば、麻酔専門医やその他の臨床医は、プロポフォールやセボフルランなどの麻酔薬を患者に投与時、本明細書に記載のシステムと方法により得られるこれら特定の指標を使って、徐波デルタ及びアルファ振動に対応する全身麻酔状態（ベータ及びガンマ振動を特徴とする鎮静状態とは異なる）を見定めることができる。

別の例では、全身麻酔又は鎮静状態を監視・追跡するために、徐波振動（０．１〜１Ｈｚ）からデルタ（１〜４Ｈｚ）振動の位相とアルファ振動（８〜１２Ｈｚ）からベータ（１２〜２５Ｈｚ）振動の振幅との間の位相振幅結合を提示し得る。当該関係の位相は、鎮静（例：谷部で最大になる「谷マックス」変調）又は深い無意識（例：ピーク部で最大になる「山マックス」変調）の状態を区別するために使用できる。結合関係の強さによって、深い無意識状態を維持するのに配送すべき麻酔剤の投与量を、閉ループ制御の下に規制することが可能である。

別の例では、ＥＥＧ及びＥＥＧにより導出された数量など、対象とする生理学的信号とは異なる二次的な生理学的信号を提示するために本システムと方法を使用することができる。例えば、低周波の一次自己回帰ＡＲ（１）成分をモデルの状態動的特性成分に追加することにより低周波のドリフトを考慮することが可能になる。同様に、運動誘発性の二次生理学的波形の状態空間モデルも含むことができる。同様に、ＡＲ又は移動平均成分をモデルの状態動的特性成分に組み込んで、筋肉誘発性二次生理学的信号を表現することもできる。このようにして、該システムと方法は、二次生理学的信号を、対象とする生理学的信号と区別するために使用できる。これらの二次生理学的信号は本質的にそれ自体が、異なる状態（意識状態又は動きなど）の指標として有用なものである。

別の作動形態において、該システムは、その他の装置（その他の生理学的監視装置など）からの情報又は医療記録システムからの情報を使用して、ＥＥＧ又は生理学的センサにより記録された脳又は生理学的状態の解釈を通知することができる。かかる情報によって、特定のオシレータモデル並びに特定の二次生理学的モデルの使用を報告し、所望の監視変数を提供することができる。

他の作動形態において、該システムは、臨床現場で発生し得る異なる仮想事態を表現することを意図したモデルのパラメータ化とモデルのデータベースを採用できる。該システムは、医療記録からの情報（患者の人口学的情報又は医療履歴など）や、投与された麻酔薬に関する情報、ならびに手術を行った症例のタイプに関する情報を使って、所定の仮想事態に適用され得る一式の候補モデルを選択できる。該システムは、監視対象の特定患者について記録されたＥＥＧ又はその他の生理学的データを使って、所定の仮想事態に最も適したモデル（複数可）を選択できる。モデルのデータベースとしては、異なる入力変数の関数やルックアップテーブル又はそれに類似した構想として表される異なる形式（異なるパラメータの事前確率分布を含む）を取り得る。

一般的に、ＥＥＧの振動は、ノイズに認められる独特な成分として表現できる。例えば、プロポフォールによる麻酔下における徐波・アルファ振動を、以下の式によりモデル化できる。

以下、更に一般的な事例を説明する。この時点で、状態空間形式で書き込める２次自己回帰（ＡＲ（２））ステップを使用して各振動成分をモデル化できる。なお、各振動成分のモデルは、後述される観測の存在により、厳密に解釈すると自己回帰移動平均（ＡＲＭＡ）モデルであるが、本明細書では、該モデルを言及する際には、モデルの振動構造を強調するために各モデルの根本的な状態動的特性に応じて説明する。差し当たり、以下のことが想定できる。

かかる２次元（２Ｄ）潜伏状態はそれぞれ独立していると想定でき、２次元（２Ｄ）潜伏状態の固定刻み幅に亘る発展は、ｊ＝ｓ，ｆの場合、スケーリングされたノイズの多い回転として次のようにモデル化できる。

式中、ａ_ｊはスケーリングパラメータであり、Ｒ（ω_ｊ）は、角度ω_ｊの２次元の回転（円振動数）であり、σ_ｊ ^２はプロセス雑音変動である。このような状態空間振動分解の一例を図１２−１の（Ａ）〜（Ｄ）に示す。この方式によると、モデルによって徐波成分と速波成分を直接推定できるので、従来の帯域通過フィルタリングが不要になる。恐らく更に重要なことは、振動の各成分を位相ベクトル又は解析信号の実質又は虚数単位と見なすことができる、という点であろう。結果として、観測された実信号から虚数信号成分を合成する従来の方法で使用されるヒルベルト変換の必要がなくなる。したがって、潜在ベクターの極座標は、徐波瞬時位相φ_ｔ ^ｓと速波振動振幅Ａ_ｔ ^ｆ（図１２−２の（Ｆ）及び（Ｇ））の直接的な表現を提供する。ちなみに、ｘ_ｔ＝［ｘ_ｔ ^ｓＴ ｘ_ｔ ^ｆＴ］^Ｔであり、以下の標準的な状態空間表現を得る。

期待値最大化法（ＥＭ）又は対数尤度又は対数事後確率密度を最大にする数値処理を用いて推定値（Φ，Ｑ，Ｒ）を得ることができる。（非線形系に遭遇する場合もあり得るので）この独立した振動成分を更に追加して又は振動の高調波を考慮してモデルを拡張するために、更に一般的な公式に基づいたものを下記のようにして導出できる。

以下、方程式（３）と方程式（４）により具現化される状態空間モデルのより具体的な派生モデルを説明する。

ｊ＝ｓ，ｆ及びｔ＝２．．Ｎの場合、各成分を次のプロセス方程式にしたがって表せる。

前述の如く、各振動の位相φ_ｔ ^ｉ及び振幅Ａ_ｔ ^ｊは、潜在ベクトル極座標を使用して求められる。

各振動は、周波数ｆｊにピークとする広帯域パワースペクトル密度（ＰＳＤ）を有する。更に、このＰＳＤのパラメータ表示を更に以下のように導出できる。

Ｍ＝［１ ０ １ ０］且つｘ_ｔ＝［ｘ_ｔ ^ｓＴ ｘ_ｔ ^ｆＴ］^Ｔと設定し、Ｑ及びΦを、めいめいがブロックＱ_ｉ及びａ_ｊＲ（ω_ｊ）のいずれかからなるブロック対角行列とした場合、方程式（３）及び方程式（４）で表される標準的な状態空間となる。

観測された信号ｙ_ｔを考慮して、隠れ振動ｘ_ｔとそれが生成するパラメータ（Φ，Ｑ，Ｒ）の両方を推定することを目標にする。そのために、期待値最大化（ＥＭ）アルゴリズムを使ってＥステップとＭステップによりかかる推定を行う。一般的に、次のように誘導される。カルマンフィルタ及び固定区間平滑化を使って、期待値最大化（ＥＭ）アルゴリズムのＥステップで隠れ振動ｘ_ｔを推定する一方、Ｍステップの各繰り返しにおいて生成パラメータを推定する。

なお、ＡＲ（２）形式では、異なる振幅、及び異なる共振又は帯域幅のレベルを表現できる。加えて、ＡＲ（２）スペクトルは、発振周波数以上の周波数にて「１／ｆ」のプロファイルを有する。上記の場合以外にも、これらの振動のいずれか又は両方が存在しない状況が発生するのも一般的である。例えば、意識がしっかりしている時の休息状態において、前頭ＥＥＧチャンネルにアルファ振動の存在、又は徐波振動の存在は期待できない。このような状況において、広帯域の「１／ｆ」動的特性が存在する場合もあり、ＡＲ（１）ステップを使用してモデル化し得る。あるいは、ＡＲ（１）とＡＲ（２）の動的特性の組み合わせを状態空間形式で表現することもできる。

以下、方程式（３）及び方程式（４）で表されるモデルにはＡＲ（２）ステップが二つあるため、該モデルを「ＡＲ（２＋２）」モデルという略称で言及する。ＡＲ（２＋２）モデルをＥＥＧデータに当て嵌めることにより、アルファ振動と徐波振動を特定できる。以下、ＡＲ（１）のステップとＡＲ（２）のステップの組み合わせを「ＡＲ（１＋２）」モデルと呼ぶ場合もある。また、ＡＲ（１）又はＡＲ（２）のいずれかの状態動的特性により表される一つの動的成分を伴うモデルについても考えてみる。上記のように、これらのモデルは厳密に解釈すると、方程式（４）における観測ノイズの存在が故に、自己回帰移動平均（ＡＲＭＡ）モデルであると解釈されるが、本明細書では、該モデルを言及する際には、モデルの振動構造を強調するために各モデルの根本的な状態動的特性に応じて説明する。シャムウェイ（Ｓｈｕｍｗａｙ）やストッファー（Ｓｔｏｆｆｅｒ）の期待値最大化アルゴリズムなど、状態空間モデル向けの期待値最大化アルゴリズムを使って、ＡＲ（２＋２）などのモデルを一つ又は複数、ＥＥＧデータに当て嵌めることができる。初期のＡＲパラメータは、ユーザにより設定可能である。例えば、初期徐波振動ピーク周波数を１Ｈｚに設定でき、初期アルファ振動ピーク周波数を１０Ｈｚに設定でき、そして、両振動の極を初期値として、複素平面の発端からの半径を０．９９に設定できる。期待値最大化アルゴリズムでは、カルマンフィルタや固定区間平滑化及び共分散平滑化アルゴリズムを利用して推定パラメータを更に絞り込むことができる。ＥＭアルゴリズムの期待ステップと最大化ステップは、ＡＲパラメータやノイズ共分散などの推定パラメータが出力される前に所定回数（例えば２００回）繰り返すことができる。その際、特定の収束基準（推定対象パラメータの収束や、又は対数尤度又は対数事後確率密度の値における収束を含むが、必ずしもこれらに限定されない）を満たすまで連続して繰り返すことができる。ＡＲパラメータとノイズ共分散の推定後、異なるモデルについて赤池情報量基準（ＡＩＣ）を計算して、各モデルが所与のＥＥＧにどれほどうまく当て嵌まるかを比較する。ＡＩＣは、ＡＩＣ＝−２ｌｏｇ（Ｌ）＋２ｐとして定義され、式中、Ｌは該モデルの尤度値であり、ｐはモデル中のパラメータ数である。

ＡＲ（２＋２）モデルが如何に正確に徐波振動とアルファ振動を決定できるかに関するＡＲ（２＋２）モデルの実行能力を判定するために特定の試験を行ったが、その試験において、１８〜３６歳の健常者にプロポフォールを投与して無意識状態を誘発し、そのＥＥＧを記録した。ＥＥＧ信号は、５０００Ｈｚのサンプリング周波数で、２００Ｈｚまでダウンサンプリングして記録した。前頭のＥＥＧチャンネルの一つを選択して、ベースライン状態（意識がしっかりしているが目を閉じている状態）を、プロポフォール誘発性の無意識状態と比較することにより分析した。プロポフォール誘発性の無意識状態では、前頭ＥＥＧに大きな徐波（０．１〜１Ｈｚ）及びアルファ（８〜１２Ｈｚ）振動が含まれる。ベースライン状態では、かかる大きな振動は認められなかったが、徐波のドリフトが認められる場合もある。ＡＲ（１）モデル、ＡＲ（２）モデル、ＡＲ（１＋２）モデル及びＡＲ（２＋２）モデルの中から最も優れたモデルを選択するためにＡＩＣを使用した。各モデルの状態空間定式化を考慮に入れるために、ＡＩＣにおけるパラメータとして、ＡＲパラメータと状態ノイズ変動と観測ノイズ変動（全て期待値最大化（ＥＭ）アルゴリズムにより推定されたもの）を含んだ。また、ＡＩＣの代わりに、その他ベイズの情報量基準（ＢＩＣ））などの情報を、最も優れたモデルを選択するために使用できる。

帯域通過フィルタを使用して、信号を周波数依存性の成分に分解できる。したがって、帯域通過フィルタ処理した信号を状態空間モデルの性能を比較した。そこで、遮断周波数１Ｈｚの、次数が１００ハミングウィンドウ適用バンドパスＦＩＲフィルタを使って、徐波を表した。通過帯域８〜１２Ｈｚの、次数が１００であるハミングウィンドウ適用低域ＦＩＲフィルタを使って、アルファ波を表した。帯域通過した徐波・アルファ成分の加算により再構築されたバンドパス信号を作成した。

以下の表Ｉは、プロポフォール誘発性の無意識状態のためのＡＲ（１）、ＡＲ（２）、ＡＲ（１＋２）ＡＲ（２＋２）モデルに関するＡＩＣスコアを示す。表からも明らかなように、これらのモデルの中でも、ＡＲ（２＋２）は、プロポフォール誘発性の無意識状態においてＡＩＣレベルが最も低く、他のモデルよりも優れた結果を示す。ベースラインの場合はＡＲ（１）モデルのＡＩＣが最も低いが、プロポフォールの場合についてはＡＲ（２＋２）モデルのＡＩＣが最も低く、上記各モデルをそれぞれの条件に合わせて選択することを可能になった。当て嵌められたモデルのスペクトルは、観測データ（スペクトル分解能３Ｈｚ、データ１０秒分につき２６個のテーパ）のノンパラメトリックマルチテーパースペクトル推定とともに図５に示す通りである。図５Ａは、プロポフォールに誘発された状態におけるＡＲ（２＋２）モデルのスペクトルを示し、図５Ｂは、休息により誘発された状態におけるＡＲ（１）モデルのスペクトルを示す。

休息状態においては、識別可能なアルファ振動はなく、ＡＲ（１）モデルの選択を支持している。休息状態に関するモデル化された観測データと収集データ間の残差は、４．９３４８ｕＶ以下の程度であった。このＡＲ（１）の場合、記録されたデータとほとんど同一の隠れ状態が一つだけ存在する。このベースライン状態のためのＡＲ（１）モデルは次の通りである。

プロポフォール誘発性の無意識状態を表すＡＲ（２＋２）モデルのスペクトルを、同じく観測データのマルチテーパースペクトル推定値とともに図５Ａに示す。ＡＲ（２＋２）モデルの徐波・アルファ振動成分は、マルチテーパースペクトルで観測されたピークと密接に対応している。ＡＲ（２＋２）モデルは徐波・アルファ振動を表す隠れ状態を備えているので、時間領域の振動を区別することが可能になる。図６はＡＲ（２＋２）モデルを使って推定した徐波振動とアルファ振動及び複合振動の時系列を示す。特に、図６Ａは複合振動時系列を示し、図６Ｂは図６Ａの徐波振動時系列を示し、図６Ｃは図６Ａのアルファ振動時系列を示す。残差は０．２１９８ｕＶ以下の程度である。

方程式（５）、方程式（６）及び方程式（７）で表される状態空間モデルのＡＲ成分の振動形式に加えて、従来のＡＲモデル定式化も使用し得る。例えば、当該プロポフォール誘発性状態における各成分を示すＡＲモデルは、次のように表せる。

この形式のモデルでは、振動の特性（ピーク周波数など）を直接モデルパラメータから計算し得る。各ＡＲ（２）成分のラジアンでのピーク周波数は次の式により表される。

式中、ａ_１は初回ラグパラメータ、ａ_２は２番目のラグパラメータである。プロポフォール下における徐波のピーク周波数は１．０１１８Ｈｚであり、アルファ振動のピーク周波数は１１．６９２８Ｈｚである。ＡＲ（１）ステップはゼロ周波数にて実数値の極を一つ有するので、ＡＲ（１）ステップの定義によると、休息状態における徐波成分のピーク周波数は０Ｈｚである。ＡＲモデル成分の極によって、パワースペクトルの形状とピーク周波数が決定される。極の半径によって、ピーク周波数の隆起が決定され、よって実質的に、振動過程の減衰が決定される。このプロポフォール下の徐波振動のためのＡＲ（２）モデルにおいて、極の半径は０．９９であり、アルファ振動用のモデルでは、半径が０．９５６であり、予想通り、徐波の方が若干大きい。

状態空間モデルの性能を、より標準的な従来の帯域通過フィルタ方式と比較した。図７は、徐波バンドパスとアルファ成分、再構築された信号、及び残差のマルチテーパースペクトルを、状態空間モデルの残差とともに示す。帯域通過フィルタは、状態空間モデルを使用して取得した振動に類似した振動を特定するが、バンドパス法による残差は著しい低周波振動構造を示す。それに比べ、状態空間モデルからの残差は約６０ｄＢほど少なく、またそれほど構造化されていない。帯域通過モデルにおいて、残差におけるパワーと構造は任意のバンドパス遮断以外における漏れの結果であるが、本発明に係る方法はＡＲモデルを使って成分の振動を表すことによりそれを避ける。図８は、時間領域におけるバンドパス法による成分の振動、再構築された信号及び残差を状態空間との比較とともに示す。図８Ａは、バンドパス法と状態空間法における総信号を示す。図８Ｂは、バンドパス法と状態空間法における徐波振動成分を示す。図８Ｃは、バンドパス法と状態空間法におけるアルファ振動成分を示す。図８Ｄは、バンドパス法と状態空間法における残差を示す。周波数領域の分析と同様に、帯域通過フィルタにより、かなりの時間領域残差構造（状態空間モデルでは存在しない）が認められる。したがって、状態空間モデルを使用することにより、特定の欠点（バンドパスに基づく方法の場合の比較的大きな残差など）無く、振動を特定できる。

本発明者らが得た結果によると、神経の時系列に振動する構造を特定するために第１次及び第２次システムの組み合わせから成る特定の形式で構成された状態空間モデルが如何にして神経細胞の時系列における振動構造を特定するために使用できるかを示す。神経科学分析において採用される典型的な方式として、スペクトルバンド内のパワー、又は帯域通過フィルタ処理された信号の変動を計算することである。この方式では、対象となる周波数帯域全体にわたる広帯域のノイズを振動と区別できない、という大きな制約が存在する。本発明の方式は、潜在的な振動成分を明示的にモデル化してから、異なる振動構成（即ち、ＡＲ（１）、ＡＲ（２）、ＡＲ（１＋２）とＡＲ（２＋２））を持つモデル一式の中から特定のモデルを選択して信号の組成を決定することにより、かかる制約の解決を図る。この方式によると、特定の振動成分に起因するパワーを、その他の広帯域成分とは独立して計算でき、成分の振動時系列を区別することが可能である。

この状態空間方式は、ドリフト（実数極を伴うＡＲ（１）又はＡＲ（２））又は振動単位（複合極を伴うＡＲ（２））のいずれかを表現するためのＡＲパラメータ構造を指定するので、従来のＡＲ又はＡＲＭＡモデル化と異なる。より一般的なＡＲ又はＡＲＭＡモデルはかかる成分を確実に特定できるが、本発明者らの経験では、かかるモデルのスペクトル表現は、モデル数が増加するにつれて非常に変わり易く、ばらつきが生じ易くなり得る。推定上のドリフト単位又は振動について特定の形式を仮定して置くことにより、かかる機能特性を表現するのに必要となるモデルパラメータの数を減らすことができ、モデル推定のロバスト性（ロバスト性）を向上できる。なおまた、本明細書において提案される構造は、より一般的なＡＲ又はＡＲＭＡ構造よりも簡単に、特定の振動成分について解釈できるものである。

状態空間法は、振動成分を時間領域で区別できるという点で、特定の方法よりも有利である。本明細書に記載の線形システム方式は、ＡＲ成分により示唆される時間的構造又はインパルス応答がオシレータに特有であるため、振動の検出に更なる特異性を提供し得る。更に、状態空間法は、未処理の（即ち、帯域通過フィルタ処理されていない）ＥＥＧデータに使用することができので、帯域通過フィルタのいくつかの欠点を回避できる。

モデルパラメータの事前分布：
特定のＥＥＧデータをモデルに当て嵌める際、観測された振動のパワーが低く、より広域の周波数スペクトルのランドスケープと区別できない状態であると、該観測された振動の一部がモデルによって正確に当て嵌められない場合もある。上記の各方法では、ユーザがモデルにおいて制御できるのは、初期化パラメータの選択のみである。低パワーの成分を表すための成分は、特に該成分が大きく広帯域である場合、当て嵌め手順中に他の振動成分が示す周波数帯域に移行する可能性がある。赤池情報量基準（ＡＩＣ）などの従来のモデル選択手順を使用して評価すると、これら過剰適合モデルは更に低いＡＩＣ値を有し得る。これは、ある程度、ＡＩＣが、より多くの成分を有する状態空間モデルに有利に働くからである。前記状態空間モデルは、初期化に敏感であり得、既に研究されている神経細胞振動の予備的知識に基づいた更なる制約を必要とし得る。

第一の代表例として、ＡＲ（２＋２）状態空間モデルを使用して、強力なアルファ振動を伴うＥＥＧデータを当て嵌める。該ＥＥＧデータのＡＲ成分のスペクトルを図９Ａに示し、推定振動のスペクトルを図９Ｂに示す。状態空間モデルは、強力なアルファ振動を伴うＥＥＧデータうまく当て嵌る傾向がある。

これに対し、アルファ振動が存在しないと、モデルをうまく当て嵌められない場合がある。アルファ振動が存在しない場合、アルファ成分が不正確に徐波デルタ周波数範囲に当て嵌められてしまう。意図した周波数帯域にモデルのアルファ（８〜１２Ｈｚ）成分を停留するための振動力がないと、モデル当嵌処理中に、当該成分の周波数がより大きな徐波成分の方向に移ってしまう。このモデルのＡＩＣは、単に徐波成分を伴うモデルのＡＩＣより低く、よって、「アルファ＋徐波」モデル（ＡＲ（２＋２））を選択すべきであることを示唆する。但し、これは、（ノンパラメトリック）スペクトル推定に明らかなアルファ振動がないので我々の直観に反する。しかるべく、アルファ成分はもはやアルファ周波数レーティングに存在しない。その代わりにアルファ成分は、１．０８７２Ｈｚをピークとする徐波周波数帯域に移動している。後述される中央周波数の事前分布のなき徐波・アルファモデルのＡＲ成分のスペクトルを図１０Ａに示す一方で、事前分布を伴う徐波・アルファモデルについて推定された振動のスペクトルを図１０Ｂに示す。また、アルファ無し徐波モデルのＡＲ成分のスペクトルを図１０Ｃに示す一方で、アルファ無し徐波モデルについて推定された振動のスペクトルを図１０Ｄに示す。

この種のモデルを更に識別し易く、且つ不適なモデルの当て嵌めに対して更にロバストにするために、事前分布を振動周波数と半径又は減衰係数のそれぞれに使用できる。

単一成分の中央周波数の事前確率分布の場合、モデルは、発生する可能性のある様々な周波数にアクセスしつつ、事前に定義されている周波数帯域に近似的に留まる（つまり、特定の周波数帯域に成分を停留させる）必要があり得える。

本発明者らは、事前確率分布としてフォンミーゼス分布を選んだ。フォンミーゼス分布は、０Ｈｚからサンプリング周波数に及ぶ領域での対称確率分布であり、次の形を取る。

式中、ωは、事前分布により制約したい振動の中央周波数を表し、Ｉ_０は、次数０の変型ベッセル関数であり、κ及びμは、分布の形状を示すハイパーパラメータである。パラメータμは、事前分布のための平均又は中央周波数を設定する。集中パラメータκは、分布の幅又は帯域幅を決定し、意図した周波数帯域における停留の強さを決定する。振動ＡＲ成分を特定の周波数範囲に停留させる目的で、その他類似した特性を有する確率分布も考慮に入れることができる。フォンミーゼス事前分布では、限定される周波数帯域に周波数成分が束縛されており、故にモデル当て嵌め又はパラメータ推定ステップ中にその他の振動成分と重複することがないため、モデルの当て嵌めよりロバストにできる。以下、試験事例を使って、フォンミーゼス事前分布を伴わないモデルが適切に機能できない理由を説明する。フォンミーゼス事前分布が無いと、ＡＩＣが当該の被験者について（実際には徐波帯域に単一振動数しかないにもかかわらず）二成分モデルを選択してしまう。とりわけ、該選択されたモデルは、徐波帯域における２つの振動に当て嵌まる。フォンミーゼス事前分布がない状態では、高周波のアルファ成分が１．０７８２Ｈｚを中心として徐波帯域にドリフトする。

中央周波数と中央への集中を制御するフォンミーゼス事前分布パラメータは、所望の周波数帯域に基づいて選択することができる。例えば、神経科学文献において一般に受け入れられている周波数帯域の幅と中心点に基づいて決定できる。例として、下の表ＩＩに示す徐波・アルファモデルは、アルファ周波数成分にフォンミーゼス事前分布を有し、確率値０．３６７にて中央周波数１０Ｈｚ、幅１．２Ｈｚ程度を使用してアルファ成分がアルファ周波数帯域に停留できるようにするための援助を提供する。

フォンミーゼス事前分布を追加することにより、成分が所望の範囲に留まり、データにアルファパワーが欠けているため、成分には全くと言って良い程パワーが割り当てられず、成分が広帯域幅を有し、ＡＩＣ適用時には単に徐波振動のみ（アルファ無し）のモデルよりも高い値を示す。かくして、ＡＩＣは、フォンミーゼス事前分布を使用するモデルの中から正しいモデルを選択する。これは、フォンミーゼス事前分布の使用が如何にしてモデルのロバスト性を向上させるのか、またそれによってどのようにＡＩＣなどのモデル選択基準に基づく正しいモデルの選択を可能にするのか、を示している。以下に示す如く、事前分布を伴わない徐波・アルファモデルとアルファ無し徐波モデル及び事前分布を伴う徐波・アルファモデルのＡＩＣ、徐波周波数、徐波半径（減衰係数）、アルファ周波数及びアルファ半径（減衰係数）を表ＩＩに示す。

同じような理由で、その他のパラメータも初期化に敏感であり得、かかるパラメータに事前分布を追加することによりモデル当て嵌めとモデル選択におけるロバスト性を高めることが可能になり、それによってより完全なモデル選択の自動化を実現できる。

かかるパラメータの一例としては半径又は減衰係数が挙げられる。半径又は減衰係数は、個々の成分の安定性を制御し、よって、全システムの安定性を制御するものである。半径又は減衰係数が１以上になると、振動時系列が不安定になり、また、時系列の振幅が無限に増大することになる。本発明では、０及び１の間のパラメータを保持するようにベータ分布で減衰係数を制御することにより、パラメータを０から１の間に維持することにより、システムの安定性を維持する。ベータ分布の形状は、二つのパラメータで制御され、０から１の領域を有し、次の方程式により表される。

式中、ａは半径又は減衰係数であり、α及びβは、事前分布の形状を決定するハイパーパラメータである。

以上説明したように、ＡＲプロセスの当て嵌めに一般に使用されるモデルの選択は赤池情報量基準（ＡＩＣ）によるものである。モデルの当て嵌め手順と選択手段は次の通りである。最大化ステップと期待ステップの繰返し処理は、期待値最大化アルゴリズムからなる。カルマンフィルタと固定区間平滑化が、推定振動時系列とも呼ばれる隠れ状態の時系列を推定するのに必要である。

独立型及び調和振動分解の期待値最大化（ＥＭ）アルゴリズム：
ここで、独立型及び調和振動を推定できる期待値最大化ＥＭアルゴリズムの代表例を説明する。モデルの当て嵌め中にＥＭアルゴリズムを使用して、所与のモデルの最善の当て嵌めを可能にするモデルパラメータを反復的に検索できる。該モデルは、主に最大化ステップ（Ｍステップ）と期待ステップ（Ｅステップ）の二段階を含む。期待ステップと最大化ステップは、最適解に達するまで所定回数（例えば２００回など）実行することができる。ＥＭアルゴリズムは、初期化されたＡＲパラメータを伴うモデルを受け入れて、該モデルのための推定ＡＲパラメータを出力できる。

ノイズ単位は全て加算性ガウス雑音と想定されるので、長さＮの時間窓当たりの完全データの対数尤度は次のように表せる。

に関して対数Ｌを最大化したいが、隠れ振動ｘ_ｔにアクセスできないので、ＥＭアルゴリズムを使用して、対数尤度（Ｅステップで）推定し推定値（Ｍステップで）最大化するという過程を交互に繰り返す。繰返し処理ｒではカルマンフィルタを使用することにより、セット

を前提としてｘ_ｔを推定し、それによって以下のアクセスを得ることができるようになる。

カルマンフィルタは、観測結果とモデルパラメータを前提にして隠れ振動を推定するために使用できる。まずは次の時点における状態を予測してから観測結果と比較する。この時点において、最新の観測データに基づく推定の更新情報を生成できる。全観測時系列を前提として、逆方向平滑化（即ち固定区間平滑化）を適用することで、全観測系列を考慮に入れるように更新を絞り込むことができる。

なお、ｔ，ｔ_１，ｔ_２＝１．．Ｎについては、以下の方程式（１８）で表されるようにし、次いで、方程式（１９）で表される順方向平滑化アルゴリズムと逆方向再帰を使ってかかる数量を計算する。

ここで、ＥＭアルゴリズムのＥステップを説明する。上記カルマンフィルタ推定法は、ＥＭアルゴリズムのＥステップに含み得る。ここで、調波成分による状態空間振動分解に関する

を説明する。単一振動は、回転行列Ｒ（ω）、スケーリングパラメータａ及びプロセスノイズ共分散行列Ｑ＝σ^２Ｉ_２×２によって定義される。以下において、ｄは、調波ｈ_１，．．．，ｈ_ｄにそれぞれ関連する独立型振動とする。また、ｊ＝１．．ｄの場合、基本周波数ω_ｊの単一振動が、Ｒ（ｈω_ｊ）、ａ_ｊ，ｈ及びＱ_ｊ，ｈ−σ^２ _ｊ，ｈＩ_２×２によりめいめい定義される総和Ｆである。なお、振動成分の総数は以下のとおりとする。

Φ及びＱは、ａ_ｊ，ｈＲ（ｈ_ｊ，ｈ）及びＱ_ｊ，ｈの対角ブロックからなるブロック対角行列である。

Ｍステップでは、Ｒ及び（Φ，Ｑ）について独自にＧ_ｒを最大化でき、次の式を書き込める。

また、Ｑはブロック対角線であるので、次の式が成り立つ。

また、Ａは左右対称であり、Φは、２×２回転行列を要素として成るブロック対角行列であるので、方程式（２３）が成り立ち、次のものが得られる。

プロセスノイズ共分散σ^２ _ｊ，ｈと、スケーリングパラメータａ_ｊ，ｈと、基本周波数ω_ｊを微分することにより次のものが得られる。

方程式（２８）の一番上の式を下の（３番目の）式内に置換すると次の式が得られる。

三角関数の公式を使用すると、方程式（２８）を最終的に次のように書き直すことができる。

総じて、ｈ_ｊ＞１のときにｊ＝１．．ｄであるため、方程式（３０）を使用して、ω_ｊを数値的に解決して、ｈ＝１．．ｈ_ｊについてａ_ｊ，ｈ及びσ^２ _ｊ，ｈを演繹する。

ｈ_ｊ＝１ならば、直ちに次の式が成り立つ。

図１１を参照すると、ＥＥＧデータを状態空間モデルと当て嵌めるステップ１１００の代表例が示されている。ステップ１１００では、上記のようにアルファ振動を当て嵌めるために事前分布を使用することと、高調波を検出できるＥＭアルゴリズムを使用することとを含み得る。状態空間モデルは、数個の振動成分と、かかる成分はいずれも、ＥＥＧデータ依存性の（つまり、患者が休息状態であるか、覚醒状態であるか、無意識状態であるかにより異なる）調波振動を含み得る。プロポフォールを用いた代表事例のための状態空間モデルは徐波、アルファ振動、及びアルファ振動の調波振動を含み得る。また、状態空間モデルは、ＥＥＧデータによっては更なる非調波振動を含み得る。例えば、患者にプロポフォールを投与した場合、ステップ１１００では、調波振動伴う又は伴わないアルファ振動と徐波振動成分を当て嵌め得る。ステップ１１００では、少なくとも一つのメモリについての指示として実施できる。その後、当該の指示を少なくとも一つのプロセッサにより実行できる。

ステップ１１００のステップ１１０４においてＥＥＧデータを受け取ることができる。ＥＥＧデータは、図２を参照して上述した患者監視装置２０２などの患者監視装置から受信できる。ＥＥＧデータを波形として解釈することもできる。その後、ステップ１１０８に移ることができる。

ステップ１１０８で、出発モデルを選択できる。一部の実施形態において、ステップ１１００において、字数の大きいＡＲモデルを当て嵌めて、最も低いＡＩＣで標的字数を特定してから、該標的字数の最も突き出た周波数ピーク（複数可）を特定してＡＲ（２）成分の出発パラメータを決定することにより、出発モデルを選択できる。プロポフォールを用いた代表事例の開始パラメータは、各ＡＲ（２）成分（つまり、振動）の初期の極（初期の徐波振動ピーク周波数と初期ピークのアルファ振動ピーク周波数を定義するもの）を含み得る。初期の極は、ピーク周波数を決定するために使用できる情報を含む。一部の実施形態において、ステップ１１００では、成分の数が最も少ないモデル（徐波振動を表現するＡＲ（２）成分を一つのみ含むモデル、など）を選択できる。ステップ１１００における後段のステップでは、成分の数が多いかどうかを判断し、多い場合は更なる成分をモデルに追加する。ＡＲ（２）成分はいずれも、モデルパラメータに一つ又は複数の事前分布を有し得る。事前分布によって、より正確にモデルを所与の振動に当て嵌めることが可能になる。以上説明したように、アルファ振動領域内のＡＲ（２）成分は、特定のアルファ振動（比較的弱いアルファ振動など）の当て嵌めを改善するために、振動周波数のフォンミーゼス事前分布を有し得る。フォンミーゼス事前分布は、確率値０．３６７にて、１０Ｈｚの中央周波数及び約１．２〜２Ｈｚの幅を有し得る。以上説明したように、モデルの減衰係数に関する事前分布を表するためにベータ分布を使用できる。ベータ分布は、０〜１の領域を有するため、減衰係数を０から１の間に維持できる。その後、モデルはステップ１１１２に移ることができる。

ステップ１１００のステップ１１１２では、モデルをＥＥＧデータに当て嵌めることができる。ステップ１１００ではアルファ周波数帯域以上で発生し得る雑音周波数を固定することによりラインノイズを隔離できる。ステップ１１００では、ＥＭアルゴリズムを使用してモデルの当て嵌めを行える。当該ＥＭアルゴリズムとしては、上記「独立型及び調和振動分解の期待値最大化（ＥＭ）アルゴリズム」の欄で説明したＥＭアルゴリズムを使用し得る。ＥＭアルゴリズムは、それぞれが一つ又は複数の高調波と関連する任意数の独立型振動を当て嵌めできる。このように、独立型振動（即ち、徐波及びアルファ振動）における高調波の存在を、たとえ高調波の近くに他の独立型振動が存在していても検出できる。以上説明したように、期待ステップとＥＭアルゴリズムを、ＥＭアルゴリズムの期待ステップと最大化ステップ、即ち、ＡＲパラメータ及びノイズ共分散などの推定パラメータが出力される前に例えば２００回といった所定回数繰り返し得る期待ステップと最大化ステップを通じて繰り返し得る。繰返し処理は、特定の収束基準（推定対象パラメータの収束や、又は対数尤度又は対数事後確率密度の値における収束を含むが、必ずしもこれらに限定されない）を満たすまで継続し得る。その後、医療従事者により、又はオシレータを分析する他のステップにおいて、オシレータ周波数などの最終的な推定パラメータを使用できる。ＥＭアルゴリズムが期待ステップと最大化ステップを所定回数繰り返した後、ステップ１１００ではステップ１１１６に移ることができる。

ステップ１１００のステップ１１１６では、更に多くの振動が存在するか否かを判断できる。一部の実施形態において、ステップ１１００では、ＥＥＧデータにホワイトノイズという理由だけでは説明できない残差や新奇性が存在するかどうかも判断できる。ステップ１１００では、予測された推定時系列（即ち上記カルマンフィルタからのｘ_ｔ ^ｔ-１）をＥＥＧデータの時系列から差し引いて、新奇性の時系列を得ることができる。ステップ１１００では、新奇性の周波数スペクトルに著しいピークが存在するかどうか判断できる。通常、ホワイトノイズには一定の周波数スペクトルがある。ステップ１１００では、新奇性の時系列の累積パワー度に関するプロットを計算し、該累積パワー度に関するプロットが、絶えず増加するパワープロットと比べて統計的な有意差を示すかどうか判断できる。新奇性とホワイトノイズの間に統計的な有意差が認められない場合、ステップ１１００では、ＥＥＧデータに更なる振動が存在せず、ステップ１１１６において当て嵌められたモデル（以下、「当て嵌めモデル」）をＥＥＧデータの振動の近似に使用すべきである、と決定する。一つ又は複数の統計的な有意差が存在する場合、該ステップでは、ＥＥＧデータに更なる振動が存在する、と決定する。また、ステップ１１００では、新奇性の累積パワー度が（通常ホワイトノイズに見られるように）絶えず増大しているか判断することにより、新規性がホワイトノイズと同じ様子かどうか否かを判断する。一部の実施形態において、ステップ１１００では、単にモデルのＡＩＣスコアを計算し、当該ＡＩＣスコアを前回（即ち、事前に実行したステップ１１１６にて）データに当て嵌められたモデルのＡＩＣスコアと比較できる。ステップ１１００では、後のステップにおいて、更なる振動が存在する場合にモデルのより良い当て嵌めを達成するために当該モデルに対して更なる振動成分を追加することになる。振動成分の追加によって各モデルのＡＩＣスコアが徐々に改善されていくと、ステップ１１００では、ＡＩＣスコアの悪化が認められるまで、振動成分を追加し続ける。現在のモデルのＡＩＣスコアが前回のモデルよりも優れている場合、ＥＥＧデータに更なる振動（複数可）がある、と該ステップにて判断され得る。現在のモデルのＡＩＣスコアが前回のモデルよりも低い場合、当該ＥＥＧデータには更なる振動（複数可）が存在せず、よって当該ＥＥＧデータにおける振動を推定するための当て嵌めモデルとして前回のモデルを使用すべきである、と該ステップにて判断され得る。その後、ステップはステップ１１２０に移ることができる。

ステップ１１００において、ステップ１１１６で更なる振動が存在しないと判断されると（ステップ１１２０において「いいえ」の場合）、ステップ１１２０からステップ１１２８に移ることができる。ステップ１１１６で更なる振動が存在しないと判断されると（ステップ１１２０において「はい」の場合）、ステップ１１２４に移ることができる。

ステップ１１００ではステップ１１２４において、更なるＡＲ（２）成分及び関連するパラメータを推定できる。字数の大きいＡＲモデルを新奇性の時系列に当て嵌めることができ、ＡＲ成分の最も高いピークに最も寄与する１対の複素極は選択できる。その後、該１対の複素極に基づいて初期パラメータを決定できる。あるいは、該更なるＡＲ（２）成分の中央周波数を残差の時系列から推定できるが、かかる推定は、残差の時系列のパワーを推定し、そして当該推定パワーに基づいてピーク周波数を推定することにより行われる。あるいは、新奇性のパワースペクトルに基づいても該中央周波数を推定できる。また、該更なるＡＲ（２）成分の減衰係数は、前回のデータセットの平均として初期化したり、又は累積パワー度プロットの形状から特定したりできる。推定された中央周波数と減衰係数は、該一対の複素極により提供される情報と同じ情報を提供できる。中央周波数且つ又は減衰係数が導出される前に、時系列を低域フィルタでフィルタ処理できる。あるいは、該ステップでは、ハラー（Ｈａｌｌｅｒ）らの記載によるＦＯＯＯＦアルゴリズムを使用し、残差に基づいて中央周波数と減速パラメータを推定できる。その後、ステップはステップ１１１２に移ることができる。

ステップ１１００では、ステップ１１２８において、推定ＡＲパラメータやノイズパラメータを含む当て嵌めモデルを出力できる。ステップ１１００では、当て嵌めモデルをユーザが見れるように表示部に出力できる、且つ又は別のステップにおける保管且つ又は使用を目的としてメモリに当て嵌めモデルを出力できる。該当て嵌めモデルは、ステップ１１２０で更なる振動が無いと判断される前にステップ１１１６で最後に当て嵌められたモデルであり得る。ステップ１１００においてＡＩＣスコアの悪化が見られるまで更なる成分が繰り返しモデルに追加されるのであれば、当て嵌めモデルは最も優れたＡＩＣスコアを有するモデルであり得る。ステップ１１００では、当て嵌めモデルを人間にとって更に見やすい状態にするべく、当て嵌めモデルに基づいて報告を生成し、該報告を出力し得る。その後、ステップ１１００を終了できる。

状態空間モデル及び上記それに関連する当て嵌め方法は、ＥＥＧデータの交差周波数分析を実行するために使用できる。交差周波数分析には、振動同士の関係、例えば、第１の振動の位相又は振幅と該振動よりも高い周波数を有する振動の位相又は振幅の間の関係の分析の如何なるものを含む。第１振動の位相が第２振動の振幅を変調するところの位相振幅結合（ＰＡＣ）が一般的な例である。ＰＡＣは、様々な周波数の範囲にわたって、又は振動の位相と広帯域の（非振動の）信号の振幅との間で、実行できる。また、本発明では、第１の振動（以上説明したように解析信号として示されるもの）の実際又は仮想の部分と、第２の振動の振幅、様々な周波数の範囲、又は広帯域信号の関係を推定する新規な方法を提示する。

ＰＡＣ分析の標準的な方式では、位相及び振幅の関係を定量するためにビン化されたヒストグラムを使用することにより、統計的非効率性の主な原因となる位相と振幅間の関係を定量化する。生信号ｙ_tは、徐波振動と速波振動を分離するためにフィルタ処理される最初のバンドパスである。その後、徐波振動の即時位相φ_ｔ ^ｓと速波振動の振幅Ａ_ｔ ^ｆを推定するためにヒルベルト変換が適用される。時間ｔにて、アルファ振幅Ａ_ｔ ^ｆが、徐波振動位相φ_ｔ ^ｓの即時値に基づいて、等間隔に（通常１８個に）隔離されている長さδψの位相のビンの一つに割り当てられる。ヒストグラムは、特定の観測時間窓Ｔ（例えば、２秒程度のエポック）にかけて構築され、それによって以下の位相振幅モジュログラム（ＰＡＭ）を得ることができる。

所与の時間窓Ｔについては、ＰＡＭ（Ｔ，．）が、速波振動振幅が徐波振動位相に対してどのように分布しているかを評価する確率分布関数である。通常は、その後に、カルバック・ライブラー情報量を利用して一様分布にて変調の強さが計測される。それによって以下の変調指数（Ｍ１）を得ることができる。

最後に、この標準的な方式では、ランダム置換なの代替データを使って、観測されたＭ１の統計的有意性を評価する。問題の動的特性に依存する間隔での分布からランダムな時間的推移Δｔを導出し、また、移行した速波振幅Ａ^ｆ _ｔ-Δｔ及び元の徐波位相φ_ｔ ^ｓを使用して位相振幅結合を推定する。その後、当該置換られた時系列のＭ１を計算し、このステップを繰り返してＭ１の帰無分布を構築する。元のＭＩが、置換られた値の９５％を超える場合、当該元のＭ１は有意性を持つと見なされる。総じて、この方法によると、モジュログラムを合理的に適切と見なされる程度に推定し、十分な同程度のデータセグメントを有意性評価のために置換えることができるようにするためには、根本的なステップが、十分な時間にわたって定常状態に維持されることが要求される。

その代わりとして、本発明では、少ないパラメータ数で交差周波数結合関係を定量することができる交差周波数分析のパラメトリック表明を導入する。そのために、以下の制約付き直線回帰問題を考慮する。

なお、かかる制約付き直線回帰問題は、最終的には以下のように書き直すことができる。

φ^ｍｏｄは、速波振動ｘ_ｔ ^ｆの振幅が図１２−２の（Ｈ）〜（Ｊ）に示される如く最大となる好ましい位相である一方、Ｋ^ｍｏｄは変調の強を制御する。図１２−１の（Ａ）は生のＥＥＧデータを示す。図１２−１の（Ｂ）は、ＥＥＧデータの６秒の時間窓を示す。図１２−１の（Ｃ）は、６秒の時間窓に当嵌まる徐波振動を示す。図１２−１の（Ｄ）は、６秒の時間窓に当嵌まる速波振動を示す。図１２−２の（Ｅ）は、状態空間モデルの代表的な形式を示す。図１２−２の（Ｆ）は徐波振動位相を示す。図１２−２の（Ｇ）は速波振動振幅を示す。図１２−２の（Ｈ）は推定されたＫ^ｍｏｄと分布を示す。図１２−２の（Ｉ）は推定されたφ^ｍｏｄと分布を示す。図１２−２の（Ｊ）は回帰推定されたアルファ波振幅を示す。例えば、Ｋ^ｍｏｄ＝１且つφ^ｍｏｄ＝０の場合、速波振動は、徐波振動のピークで最も強くなる。一方、φ^ｍｏｄ＝πの場合、速波振動は徐波振動の谷又は底で最も強くなる。

最後に挙げるべき点として、代替データに依存して統計的有意性を決定すると更に効率が下がるので、その代わりに、本発明のモデル定式化では、図１２−２の（Ｆ）〜（Ｊ）に示す如く、変調パラメータの事後分布ｐ（Ｋ^ｍｏｄ，φ^ｍｏｄ｜｛ｙ_ｔ｝_ｔ）を推定して、関連する確信区間（ＣＩ）を推定することを可能にする。

本発明者らは、本発明の方式を状態空間交差周波数分析（ＳＳＣＦＡ：Ｓｔａｔｅ−Ｓｐａｃｅ Ｃｒｏｓｓ−Ｆｒｅｑｕｅｎｃｙ Ａｎａｌｙｓｉｓ）法と呼ぶ。生理学系は時間とともに変化する。そのため、複数の重複しない時間窓に生理学系を適用する。本発明に係る方法の変形においては、第２の状態空間モデルを使用することにより様々な窓にわたって変調パラメータに対して時間的連続性制約を課すこともでき、それによって、本発明者らが二重状態空間交差周波数分析（ｄＳＳＣＦＡ：Ｄｏｕｂｌｅ Ｓｔａｔｅ−Ｓｐａｃｅ Ｃｒｏｓｓ−Ｆｒｅｑｕｅｎｃｙ Ａｎａｌｙｓｉｓ）と呼ぶものを得ることができる。

本発明に係る方法の性能を実証するべく、まずは、図１４に示される如くプロポフォールを投与して鎮静及び無意識状態を誘発した被験者ボランティアにおけるＥＥＧデータを分析した。図１４Ａは患者のＥＥＧデータに供給したプロポフォール濃度を示す。図１４Ｂは、患者のＥＥＧデータの反応確率を示す。図１４Ｃは、患者のＥＥＧデータの変調回帰値を示す。図１４Ｄは、患者のＥＥＧデータのパラメトリックスペクトルを示す。図１４Ｅは、患者のＥＥＧデータのモジュログラムＰＡＭを示す。図１４Ｆは、患者のＥＥＧデータの変調パラメータを示す。予想通りに、プロポフォール濃度が増加するにつれて、被験者の聴覚的な刺激に対する反応の確率が減少した。この間、図１４Ｄに示される如く、反応確率が減少し始めるにつれてパラメトリックパワースペクトル密度（方程式（５９）を参照）が変化して、ベータ（１２．５〜２５Ｈｚ）振動が生じ、その後、反応確率が０（ゼロ）になった時に徐波（０．１〜１Ｈｚ）とアルファ（８〜１２Ｈｚ）振動が生じる。時間窓Ｔについては、図１４Ｆに示される如く、ｄＳＳＣＦＡを使って変調強度Ｋ^ｍｏｄ位相φ_Ｔ ^ｍｏｄ（及びＣＩ）を推定し、これらの推定を位相振幅モジュログラム（図１４ＥのＰＡＭ（Ｔ，．））に集める。所与の時間窓Ｔにおいて、ＰＡＭ（Ｔ，）は、関数の台［−π，π］を有する確率密度関数（ｐｄｆ）であり、それによって、速波振動の振幅が徐波振動の位相に対してどのように分布しているか評価する。反応確率が０（ゼロ）の時、速波振動振幅が徐波振動のピークで最大になる強い「山マックス」Ｋ^ｍｏｄ≒０．８，φ_Ｔ ^ｍｏｄ≒０）パターンが認められる。無反応状態への又は無反応状態からの遷移中には、速波振動振幅が徐波振動のピークで最大になる「谷マックス」（山マックスより弱い）パターン（Ｋ^ｍｏｄ≒０．２５，φ_Ｔ ^ｍｏｄ＝±π）が認められる。本発明に係るモデルは、結合が強い場合にＡ_ｔ ^ｆをより正確に予測できるので、変調関係の決定係数Ｒ^２が、結合強さＫ^ｍｏｄを正確に映し出すことに留意されたい。

継続的に変化せず安定した長期時間窓を通じて平均すると、従来の方法はＰＡＣの定性的評価に優れている。しかし、多くの場合、実験条件や臨床的症状が急激に変化するような状況においては、短い時間窓における分析を要し得る。過去の研究では、比較的長いδｔ＝１２０秒の時間窓で従来の方法を使用してＰＡＣを分析した。この場合、１４分程度の間隔をおいて一定の速度でプロポフォールを投与したのでかかる時間窓が適切と見なされる。ＳＳＣＦＡ法とｄＳＳＣＦＡ法では統計的効率が向上されるので、δｔ＝６秒といったはるかに短い時間窓を対象とした分析が可能になる。患者２名から得たかかる分析の結果を図に示す。ちなみに、被験者２名のうち１名においては（図１５に示すように）カプリング強度が強く、他の１名においては（図１６に示すように）カプリング強度が弱かった。図１５Ａは、反応確率を示すプロット図である。図１５Ｂは、プロポフォール濃度を示すプロット図である。図１５Ｃは、６秒の時間窓を使った標準的な方式の場合のモジュログラムを示す。図１５Ｄは、図１５Ｃに関連する変調指数を示すプロット図である。図１５Ｅは、１２０秒の時間窓を使った標準的な方式の場合のモジュログラムを示す。図１５Ｆは、図１５Ｅに関連する変調指数を示すプロット図である。図１５Ｇは、ＳＳＣＦＡ法の場合のモジュログラムを示す。図１５Ｈは、図１５Ｇに関連する変調指数を示すプロット図である。図１５Ｉは、ｄＳＳＣＦＡ法の場合のモジュログラムを示す。図１５Ｊは、図１５Ｉに関連する変調指数を示すプロット図である。図１６Ａは、反応確率を示すプロット図である。図１６Ｂは、プロポフォール濃度を示すプロット図である。図１６Ｃは、６秒の時間窓を使った標準的な方式の場合のモジュログラムを示す。図１６Ｄは、図１６Ｃに関連する変調指数を示すプロット図である。図１６Ｅは、１２０秒の時間窓を使った標準的な方式の場合のモジュログラムを示す。図１６Ｆは、図１６Ｅに関連する変調指数を示すプロット図である。図１６Ｇは、ＳＳＣＦＡ法の場合のモジュログラムを示す。図１６Ｈは、図１６Ｇに関連する変調指数を示すプロット図である。図１６Ｉは、ｄＳＳＣＦＡ法の場合のモジュログラムを示す。図１６Ｊは、図１６Ｉに関連する変調指数を示すプロット図である。

そのために、δｔ＝１２０秒又はδｔ＝６秒で用いられたＳＳＣＦＡとｄＳＳＣＦＡ及び標準方法をモジュログラムと変調指数（Ｍ１）の推定に基づいて比較した。後者では、いずれかの時間窓Ｔについて、ＰＡＭ（Ｔ，．）がどれほど一様分布と異なるかを測定することにより変調の強さを評価した。かかる評価には、通常カルバック・ライブラー情報量を使用した。したがって、推定されたＰＡＭにおける不規則な変動の如何なるものによってＭ１が増加し、それによってバイアスが発生する。本発明に係るモデルパラメータ化をＰＡＭとＭ１及び関連するＣＩを導出するために使用するが、標準的なノンパラメトリック分析については一般的にビン化ヒストグラムに依存する。その結果として、代替データセットを構築してｐ値を報告することにより統計的有意性を推定する。

前記被験者２名は両者とも、無意識状態に又は無意識状態に移る時に、前記の典型的な位相振幅交差相関プロファイルを示した。ＳＳＣＦＡ法は、短期の時間窓においても、より効率的に位相と振幅を推定して、平滑化されたＰＡＭ推定を作成するので、ＳＳＣＦＡ法により導出されたＭＩ推定は標準的な方式よりもバイアスが少ない。同様の理由で、φ^ｍｏｄの推定も、標準的な方式の場合よりも変動が少ない。ｄＳＳＣＦＡアルゴリズムは、ＰＡＭに対して時間的連続性制約を加えるので、更にＰＡＭ推定の変動を抑えつつ、ＰＡＣの時間依存性変化を追跡することが可能になる。最後に挙げるべき点として、本発明に係るパラメトリック戦略は、Ｋ^ｍｏｄ、φ^ｍｏｄ及びＭＩの事後分布を提供するので、変数ごとにＣＩを推定することが可能になるとともに、代替データによる方法に頼ることなく有意性を評価できる。

動物モデルにおける侵襲的記録の代表例として異なる仮想事態を背景に、本発明に係る方式の性能を説明するために、シータ（６〜１０Ｈｚ）振動及び低ガンマ（２５〜６０Ｈｚ）振動が関与することを仮定したうえで、所与の課題を学習中のラットから得たデータを分析した。ｄＳＳＣＦＡを２秒の時間窓に適用し、ラットが弁別課題を学習するにつれて海馬のＣＡ３領域におけるシータ／ガンマカプリングが増加したことを確認した。ｄＳＳＣＦＡ法を使った分析を、対象となる周波数をあらかじめ選定したり、帯域通過フィルタ処理遮断周波数を指定したりせずに行った。むしろ、該データによる位相と振幅について、根本的な特定の振動周波数をＥＭアルゴリズムにより（初期始点がシータ及びガンマの範囲にあることを前提として）推定できた。したがって、本発明に係る方法をＬＦＰデータの分析に効果的に適用できることと、本発明に係る方法では、一定の周波数や周波数帯域を指定することなく根本的な振動周波数を特定できることと、を次に説明する。

複数の模擬データセットを分析して、異なる変調機能と信号対雑音比レベルの関数としてより系統的な方法で本発明に係るアルゴリズムの試験を行った。設計上、かかる模擬データは、状態空間オシレータモデルとは異なる生成ステップ又はモデルを使って構築された（つまり、当該の模擬データ生成処理は、本発明に係る方法で使用される「モデル等級」に含まれない）。ここで、徐波成分とアルファ成分に焦点を置いて本発明に係る主な実験データの事例を再現するが、その目的は、当が信号対雑音比、信号形状などに強く依存するので、本発明に係る意図は、本発明に係るアルゴリズムの精度や正確さなどの特徴を完全に定義することではなく、むしろ、短期且つノイズを多く含む時間変動性のデータセットの場合において本発明に係るアルゴリズムの性能が如何にして且つどのような理由で標準的な分析法の性能を上回るかを例示的に説明することを標的としている。

まず、複数の異なる時間スケールで変調パラメータを使った広帯域信号におけるｄＳＳＣＦＡの分解能とロバスト性を従来の手法と比較する。その結果を、異なる生成パラメータ（後述を参照のこと。ｆ_ｓ ^ｇｅｎΔ＝ｆ_ｆ ^ｇｅｎ、σ_ｓ及びσ_ｆ並びに関連する信号）について報告する。標準的な技法は、結合パラメータが定常的で長期にわたる時間窓で平均した場合にはロバストであるが、変調パラメータが時間窓を隔てて急速に変化すると効果を失う。時間窓が長いと変調を解決できない。しかし、窓のサイズを減らして補正すれば、推定のバラツキを著しく削減出来る。妥協して両立できる方法を実験に基づいて見出す必要がある。その一方、ＳＳＣＦＡもｄＳＳＣＦＡも、６秒の時間窓で適用した場合、信号対雑音比が低くても振幅変調の迅速な変化を追跡できることが判明している。また、ｄＳＳＣＦＡアルゴリズムは、変調パラメータの事後分布の推定も提供するので、ＣＩの構築と統計的推定を単刀直入に行えるようになる。それに比べて、代替データ方式では、シャッフルできるデータセグメントがどんどん少なくなるので実施できない。

ある研究の報告では、非線形ＰＡＣ定式化が、変調された信号が徐波振動の多項式関数であるところの被駆動型自己回帰（ＤＡＲ）過程として説明されている。後者は、ドライバとも呼ばれ、フィルタ処理により事前設定周波数周辺の観測から取り除かれて、ＤＡＲ係数を推定するために使用される。次いで、信号のパラメトリックスペクトル密度を、徐波振動の関数として導出する。その後、速波振動パワーが優先的に分布される位相によって変調を表現する。ドライバについてグリッド検索を行い、様々な高速周波数にわたって各徐波中央周波数のモジュログラムを得る。その後、最も可能性の高い且つ又は最も強い結合関係に関連した周波数が、最終的な結合推定として選択される。

このようなパラメータ表現は、特に信号時間窓が短い場合に効率を向上させるが、初期フィルタ処理ステップに依存するので、従来の技法と同じ限界をいくつか有する。後述の如く、不要波のＣＦＣは、突然の変動を示す信号又は非線形性から現れる可能性がある。加えて、初期フィルタ処理ステップにおいて、短いデータセグメントからのＰＡＣ推定に広帯域の徐波振動が混入し得る。

本発明に係る方法を標準的な技法及びＤＡＲ法と比較するために、異なる対象周波数（ｆ_ｓ及びｆ_ｆ）とスペクトル幅（Δｆ_ｓ ^ｇｅｎ）及び信号対雑音比（ＳＮＲ）を使って変調信号（方程式（８８）、λ＝３及びφ^ｍｏｄ＝−π／３）を生成した。方程式（８７）を使用して、当該生成パラメータの典型的な信号トレースを生成できる。その後、これらの方法について、どれほどうまく徐波振動及び速波振動又は変調位相を回収できるかを比較した。ここに記載されるその他の方法とは反対に、ＳＳＣＦＡ法は、対象とすべき周波数を選択するために、コモジュログラム（共変調グラフ）を全て計算せず、むしろ状態空間オシレータモデルの当て嵌めによって特定する。結合は、第２ステップにおいてのみ推定される。前記において、アルゴリズムを初期化するために具体的な予備知識を使ったが、公正な比較を提供するために初期化手続きを採用する。各条件について、６秒の時間窓を４００個生成した。必要に応じて、Δｆ_ｓ ^ｆｉｌｔ＝Δｆ_ｓ ^ｇｅｎを使ってドライバを抽出した。

本発明に係るアルゴリズムによって、特に徐波振動の帯域がさらに広い場合に、高速周波数をより好ましく回収できることが判明した。また、本発明に係るアルゴリズムは、変調位相の推定において他の方法を上回る。（本発明に係るアルゴリズムは、広帯域（Δｆ_ｓ ^ｇｅｎ＝３Ｈｚ）の又は弱い（σ_ｓ，σ_ｓ）＝（０．７，０．３））徐波振動の場合に安定しており、考え得るほぼ全ての状況において異常値が極めて少なく標準偏差も小さい状態で正確にφ^ｍｏｄを推定できる。

ＣＦＣが神経系の調整において中心的な役割を果たす可能性が高いにもかかわらず、標準的なＣＦＣ分析法は、未解決のままとなっている懸念をもたらす多くの条件により影響される。基本となる信号に著しい遷移又は非線形性がある場合に不要波の結合が生じ得る。一方、帯域通過フィルタ処理の周波数帯域が適切に選択されていないと、真の基本的な結合を見逃す場合もある。ここで、ＳＳＣＦＡ法が如何にして、かかる全制約に対してロバスト性を発揮するかを説明する。また、本発明に係る方法が、線型モデルの使用により、直観に反するように信号の非線形特性を如何にして抽出し得るのかも説明する。

振動神経波形は、急激に変化したり非対称になったりして、狭帯域に閉じ込められないことを特徴とする。そのような場合、スペクトルのコンテンツを標準の帯域通過フィルタで切り捨てると、信号の形状がゆがんで、誤って結合と解釈された可能性のあるアーティファクト的成分が導入される可能性がある。

状態空間オシレータモデルは、帯域通過フィルタ処理に取って代わる、非正弦波に対応できるものを提供する。本項では、明示的に徐波振動信号の高調波を表現するようにモデルを拡張することにより、該モデルが、急激な遷移がある振動又は非線形系により生成された可能性のある振動をより好ましく表現できるようにする。そのために、同じ基本周波数ｆ_ｓに関してｈ振動を最適化する。これについては、さらに詳しく後述する。さらに、（ｉ）徐波高調波ｈの数、及び（ｉｉ）速波振動の有無を決定するために、このモデルを情報基準（赤池情報量基準（ＡＩＣ）又はベイズ情報量基準（ＢＩＣ））と組み合わせる。ΔＩＣ＝ＩＣ−ｍｉｎ（ＩＣ）を最小化することにより最も適したモデルを選択する。ＡＩＣもＢＩＣも同様の性能を発揮するが、ここではＡＩＣに基づくＰＡＣ推定のみを報告する。複数の徐波高調波が好ましい場合には、基本振動の位相を使用してＰＡＣを推定する。

まず、観測ノイズ（ν_ｔ〜Ｎ（０，Ｒ）、√Ｒ＝０．１５）を追加したファン・デル・ポール振動子（方程式（８９）、ε＝５，ω＝５_S ^-1）を使って非対称的で突然の変動を示す信号をシミュレーションした。その後、変調された高速正弦波を伴う場合（図２７Ａ、Ａ_ｔ ^ｆ＝Ａ_０（１＋ｃｏｓφ_ｔ ^ｓ）、Ａ_０＝２√Ｒ及びｆ_ｆ＝１０Ｈｚ）と伴わない場合（図２７Ｂ）の２種類の仮想事態を検討した。本発明に係るモデルは著しい遷移に対応できるので、ＡＩＣとＢＩＣのいずれによっても正しい数の独立成分を特定できる。結果として、明らかな速波振動を検出されないと、ＰＡＣの計算は行われない。一方、速波振動が存在しない場合、標準的な技法によって、突然の変動を示す信号徐波振動に由来する速波成分を抽出し、結果として、不要波の結合が検出されることになる。

信号伝導高調波から生じる非線形入力は、ＣＦＣ分析においても同様に乗り越えるべきハードルである。徐波振動ｘ_ｔ ^ｓ＝ｃｏｓ（ω_ｓｔ）が、ｙ_ｔ＝ｇ（ｘ_ｔ ^ｓ）として非直線的に変換されたとすると、次のような第二次近似が成立する。

ω_ｓ／２π＝１Ｈｚの場合、ｆ_ｆ＝２Ｈｚをピークとする振動を抽出するためにｙ_ｔを０．９〜３．１Ｈｚのあたりで帯域通過フィルタ処理することにより次のものが得られる。

かかる場合において、振動分解によってＣＦＣを使わずに高調波成分を正確に特定できる一方、標準のＣＦＣ分析では、かなりの結合が伴うことが察しられる。

狭帯域成分よりも（拡張された）状態空間オシレータの方が、生理学的信号のモデル化に適していることが認められた。加えて、信号のコンテンツ（例えば、プロポフォール誘発性麻酔の徐波・アルファ又は齧歯動物の海馬に見られるシータガンマ振動）に関する予備的知識と組み合わせたモデル選択パラダイム（理論枠組み）によって、原則に基づいてより理にかなった方法で交差周波数分析を研究することが可能になる。

過度な狭帯域幅の帯域通過フィルタを変調信号に適用すると、変調構造が壊滅される可能性がある。本発明に係るＳＳＣＦＡアルゴリズムは、変調サイドローブを生じさせる構造的関係を明示的にモデル化しない状態空間オシレータ分解を使用する。試験中、当て嵌められた時系列とスペクトルに認められるように、変調がうまく抽出されることが判明した。図１７（Ａ）は、標準的な方法による振動成分の成分の分解を示し、図１７（Ｄ）は、所与の信号に関するＳＳＣＦＡ方法による振動成分の成分の分解を示す。図１７（Ｂ）は、標準的な方法を用いて求められたパワースペクトル密度を示し、図１７（Ｅ）は、ＳＳＣＦＡ法を用いて求められたパワースペクトル密度を示す。図１７（Ｃ）は、標準的な方法による位相振幅変調推定を示し、図１７（Ｆ）は、所与のＥＥＧ信号に関するＳＳＣＦＡ法による位相振幅の推定を示す。ＳＳＦＣＡ法で使用される振動分解方法は、ｄＳＳＦＣＡ法で使用される振動分解方法と同じである。これは、モデルにおいて、速波成分の周波数反応を、サイドローブを包含するのに十分に広げることにより実施できる。これは、該アルゴリズムが、ノイズ共分散Ｒとσ_ｆ ^２及びσ_ｆ ^２を急激且つ大幅に増加させる一方で、ａ_ｆを急激かつ大幅に収縮させることによって実現される。ＡＲＭＡ（４，２）（２つの振動の積を表し、極がａ_ｓａ_ｆｅ^{±ｉ（ωｆ±ωｓ）}にある）のような高次モデルを使用することにより、又は非線形観測により直接結合をモデル化することが可能になり得る。但し、いずれの場合においても、かかるモデルを、データに当て嵌めづらく、ノイズを多く含み非定常的な非正弦波の生理学的信号に適用された時に十分に律則できなくなる。それよりむしろ、本発明に係るよりシンプルなモデルの方が、変調された高周波数成分をロバストに引き抜けることが分かった。

ＥＭアルゴリズムでは収束を保証できるが、最大化すべき対数尤度は必ずしも凹形とは限らない。字数が

である最も適した自己回帰（ＡＲ）ステップのパラメータにより、ｄ振動からなる信号成分を初期化できるそれでもなお、かかる処理においては、電気生理学的信号の非周期成分が理由で、初期化にバイアスが生じ得る。実際、非周期成分は通常、１／ｆ^ｘのべき乗関数（ＡＲステップにより回帰推定し得るもの）で表される。そのような場合、初期化において、根本的な実際の振動を考慮できない場合がある。

このような潜在的な問題を和らげるために、ハラーらのアルゴリズムを状態空間振動フレームワークに採用する。本発明に係る初期化アルゴリズムは、結果として得られるスペクトルと、更に詳しく後述される振動（方程式（ＺＺＨＨ））のパラメトリックパワースペクトル密度（ＰＳＤ）とを当て嵌める前に、振動成分を非周期的成分から解き放つことを目的とする。

観測データ信号ｙ_ｔのＰＳＤは、マルチテーパー法を使用して推定できる。その後、一般的に周波数分解能ｒ_ｆを１Ｈｚに設定する。この場合、

の時間帯域幅積が得られる。その後、Ｋ＜＜［２ＴＷ］＜−１の条件を満たすようにテーパ数Ｋを選択する。

観測ノイズＲ_０（Ｒの初期化に使用されたもの）は、次の方程式により推定できる。

その後、観測ノイズＲ_０をＰＳＤから取り除くことができる。

非振動成分は、ＰＳＤから回帰推定できる。その後、周波数ｆのｄＢにおける非周期的信号ＰＳＤを次の方程式によりモデル化する。

は、非周期的信号の勾配を、ｇ_０はオフセットを、そしてｆ_０は「ニー（ｋｎｅｅ）」周波数を制御する。非振動成分に相当する周波数を特定するために第１のパス当て嵌めを適用する。図１３Ａに示すようにχとｇ_０をそれぞれχ＝２とｇ_０＝ＰＳＤ（ｆ＝０）に設定する間、ｆ_０のみを当て嵌める。図１３Ｂに示すように、非周期的信号に関連した周波数を特定するために閾値を固定する（通常、残差の０．８分位）。その後、図１３Ｃに示すように、ｆ_０とχとｇ_０を演繹した周波数にのみ第２のパス当て嵌めを適用する。図１３Ｄに示すように、ｄＢにおける生ＰＳＤからｇ（ｆ）を取り除いて、アルゴリズムの第２ステップのために矯正済ＰＳＤを生成する。その後、当て嵌められたパラメータを使用して、ＥＭアルゴリズムを初期化できる。

矯正済ＰＳＤから、方程式（ＺＺＨＨ）により提供される理論的ＰＳＤを使用して、独立型振動の所与の数値ｄ_０（例えば、ｄ_０＝４）を当て嵌める。そうするために、十分な幅（ｒ_ｆ／２以上の幅）のＰＳＤピークを、当該ピーク周辺の幅２ｒ_ｆの近隣における振動の理論的スペクトル）を当て嵌める前に特定する。振動ｊについては、

を演繹する。

は、非周期的成分を削除した後の所与の振動のオフセットを表すので、それを調節してσ_ｊ ^２を推定する。

その後、結果として得られるスペクトルＰＳＤ_ｊを差し引いて、図１３ＥのＤＳＳ初期化当て嵌め線１３００によって示されるように全ての振動が推定されるまでステップが繰り返される。そして最後に、（ｆ_ｊ）_０の近隣における振動ｊのパワーＰ_ｊを推定し、次の方程式により総パワーＰ_０に対するそれの寄付を推定する。

振動が仕分けられ、結果として得られるパラメータを、

の第１振動によりＥＭアルゴリズムを初期化するために使用する。

統計的効率を上げるために、ＰＡＣのパラメータ表現を導入する。所与の時間窓について、次の（制約付き）直線回帰問題を考慮する。

方程式（４１）は次のものと同等であることが分かる。

と設定することにより、モデルの一貫性を保証できる、つまり、変調包絡線が搬送波信号の振幅を超えることが決してないようにできる。しかし、これを達成するには、計算コストの高い制約を課す必要があり得る。データの信号対雑音比が高く、Ｋ^ｍｏｄが偶然にも１を超える可能性が低い場合には、無制約問題

を解決することにより達成することも可能である。束縛解の下では、βの事後分布が切断多変量ｔ分布である。

尤度、共役事前・事後パラメータβ、Ｖ、ｂ、ｖ及び正規化定数Ｚは、次の方程式に基づいて更に正当化され導出される。ここでは、かかる推定を状態空間相互周波数分析（ＳＳＣＦＡ）とし、次のように表す。

標準的な方式は、統計的有意性を決定する際に代替データに依存するが、それによって効率が更に低下する。本発明は、その代わりに、変調パラメータＫ^ｍｏｄとφ^ｍｏｄの確信区間（ＣＩ）を得る元となる事後分布ｐ（β｜｛ｙ_ｔ｝_ｔ）を推定する。そして、事後分布を推定するために、前記のような（ｉ）状態空間オシレータモデルと（ｉｉ）パラメトリックＰＡＣにより提供される事後分布から標本をサンプリングする。

（ｉ）上記ＥＭアルゴリズムのｒ番目のＥステップにおいて使用したカルマンフィルタによって、ｔ，ｔ’＝１．．Ｎの場合のモーメントが次のように求められる。

（ｉｉ）、各χについて、方程式（８）を使用することにより、再サンプリングした徐波振動の位相φと速波振動の振幅Ａを計算する。その後、方程式（４４）を使って、ｐ（β｜Ａ，φ）からｌ_２個のサンプルを導出する。その結果、ｌ_１×ｌ_２個のサンプルが作成されて、次の推定が行われる。

最後に、Ｌ_２ノルムを使用してβ^{ＳＳＣＦＡ}周辺にＣＩを構築し、それをもって、Ｋ^ｍｏｄとφ^ｍｏｄのＣＩを導出する。

前記分析を適用する対象となる複数の重複しない時間窓（長さ＝Ｎ）に時系列を分ける。故に、変調を考慮に入れた

における一式のベクトル{β^{ＳＳＣＦＡ}}_Ｔ（式中、Ｔは長さＮの時間窓を示す）を求める。

第２の状態空間モデルを使って変調の動的特性を表現できる。ここで、観測ノイズを有する次数ｐの自己回帰（ＡＲ）モデルを、時間窓にわたって変調ベクトルβ^{ＳＳＣＦＡ}に当て嵌める。それによって、二重状態空間相互周波数分析（ｄＳＳＣＦＡ）による推定を得ることができる。

ここでは、ユール−ウォーカーのような方程式を数値的に解決且つ最適化するにより次に進む。これについては更に詳しく後述する。そして、ベイズの情報量基準で次数ｐを選択する。最後に、当て嵌め済のパラメータを使用して、ｌ_１×ｌ_２個の再サンプリングされたパラメータをフィルタ処理でき、それによって必要に応じ｛β^{ＳＳＣＦＡ}｝_ＴのＣＩを構築できる。

標準的な技法をＳＳＣＦＡとより良く比較するために、時間窓Ｔについて、本発明に係るパラメトリックモデルによるＰＡＭの近似式を導出する。

プロポフォールを麻酔剤として投与した被験者の意識の消失及び回復中に当該患者のヒトＥＥＧデータを分析した。短期間ではあるが、健常者ボランティア１０名（１８〜３６歳）に、プロポフォールを、６段階の標的効果部位濃度範囲（０μｇ／Ｌ、１μｇ／Ｌ、２μｇ／Ｌ、３μｇ／Ｌ、４μｇ／Ｌ及び５μｇ／Ｌ）で徐々に投与量を増やしながら注入投与した。注入投与はコンピューター制御のもとに行い、各濃度をそれぞれ１４分間維持した。意識の消失及び回復を行動的に監視するために、被験者に対して４秒ごとに聴覚刺激（クリック音又は口頭命令）を与え、それに対して、被験者はボタンを押すことによって反応を示すように指示を受けた。モンテカルロ法によりかかるデータに状態空間モデルを当て嵌めることにより、反応の確立と関連する９５％の確信区間を推定した。最後に、ＥＥＧデータを、アンチエリアジングフィルタにより前処理して、２５０Ｈｚまでダウンサンプリングした。

振動分解に関連するパラメータ表示を使用してＥＥＧのスペクトログラムを計算した。０．１〜１Ｈｚ及び９〜１２Ｈｚのあたりで帯域通過フィルタの使用によりアルファ及び徐波成分特定且つ抽出されると想定されるところの６秒と１２０秒の時間窓で標準的なＰＡＣ分析法を適用した。標準的なＰＡＣ法の有意性を、ランダムな置換（パーミュテーション）を２００回行うことにより評価した。

以下、同様のスペクトルのコンテンツで自己回帰移動平均過程（ＡＲＭＡ）を構築することにより、振動ｘ^ｊ _１，ｔのパワースペクトル密度（ＰＳＤ）のパラメータ表示を導出する。便宜のため、以下の説明においては添え字のｊを省略する。ちなみに、振動は漸近的に二次定常性であるということをまず述べておく。では、それの自己共分散列を計算してみる。Ｒは回転行列であるので、Ｒ（ω）^ｋ＝Ｒ（ｋω）及びＲＲ^Ｔ＝Ｉとなる。したがって、方程式（２）から、次のものが得られる。

故に、

について、書き込むことができる。結果として、二次定常ステップにより振動を近似することができる。また、それの理論的パワースペクトル密度は、ウィーナー・ヒンチンの定理によると次のようになる。

ここで、以下のＡＲＭＡ（２，１）を考慮する：

上記のものから、階数２の二つの負ルートΨ^± _１を演繹し、最終的に同じ自己共分散系列につながる。総じて、本発明では次のものを選択する。

離散的フーリエ変換を方程式（５４）収穫率に適用すると、次の方程式が得られる。

最終的に、ｆ_０を中心とした振動のＰＳＤを次のように表せる。

本発明ではτ_β＝σ_β ^-２を使用し、方程式（４１）に定義されるモデルの尤度を次のように想定する。

式中、｛Ａ_ｔ ^ｆ｝_ｔ＝Ａ_ｆ及び｛φ_ｔ ^ｓ｝＝φ^ｓとする。共役事前分布は次のとおりである。

ここで、所与のＡＲ次数ｐについて、方程式（４９）により定義されたパラメータＲ_β、Ｑ_β及び｛ｈ_ｋ｝^ｐ _ｋ＝１を推定する。方程式（４５）で推定した変調ベクトルの自己共分散列を

とする。結果として、次の方程式が成り立つ。

観測ノイズ候補Ｒ_β ^＊については、方程式（７３）を反転できれば、直ちにＱ_β ^＊及び｛ｈ_ｋ ^＊｝^ｐ _ｋ＝１にアクセスする。よって、カルマンフィルタを使用して、候補モデルの尤度を演繹する。

したがって、テプリッツ行列Ｃ＝（Ｃ_{｜ｉ−ｊ｜}）_{ｉ，ｊ＝０..ｐ}最小の固有値Ｒ_β ^ｍを入力し、（Ｃ−ＩＲ_β ^＊）が最大階数であると判明している（０，Ｒ_β ^ｍ）におけるＲ_β ^＊に対するモデル尤度を数値的に最適化する。

Ｒ_βからＱ_βと｛ｈ_ｋ｝^ｐ _ｋ＝１を得てから、再びカルマンフィルタを使用してβ_Ｔ ^{ｄＳＳＣＦＡ}を推定する。

ここで、τを中心とした長さδｔ＝Ｎ／Ｆ_Ｓの時間窓について、近似されたパラメータのモジュログラムを導出する。その際、次の条件を使用する。

明確にするために述べておくが、ｔ又はκを無差別に使い、次の方程式を念頭に置いて説明する。

シミュレーション：
本発明は、異なるモデルにより生成された模擬データセットを使って、本発明に係るＳＳＣＦＡ法とｄＳＳＣＦＡ法の試験を行った。単位分散ガウス雑音、（別途記載なき限り）ｆ_Ｓ＝１Ｈｚを中心とする徐波振動、及び（別途記載なき限り）ｆ_ｆ＝１０Ｈｚを中心とする変調徐波振動を組み合わせることにより、各シミュレーション信号を構築した。但し、ここで注記すべき重要な点は、これらのシミュレーション信号を生成する際に、本発明においてデータを分析するために使用する状態空間オシレータモデルとは異なる方法又は「モデル等級」を使用することにした、ということである。標準処理については、ランダムな置換（パーミュテーション）を２００回行うことにより有意性を評価し、ｆ_Ｓとｆ_ｆが分かっているという前提で、通過帯域が余波成分の場合は０．１〜１Ｈｚに、速波成分の場合は８〜１２Ｈｚに設定された帯域通過フィルタを使って成分を抽出した。

徐波振動のシミュレーション：

神経細胞振動は完璧な正弦波でなく、むしろ、広帯域特性を有する。本発明は、互いに独立で同一の分布に従うガウス雑音を、次のインパルス応答関数で畳み込む（フィルタ処理する）ことにより広帯域徐波振動をシミュレーションした。

最後に、結果として得られる系列が正規化され、それによって標準偏差がσ_Ｓに設定される。

本発明に係る方法及び標準方法の時間的分解能を評価するために、異なる時間的変動性変調速度で模擬データセットを生成した。まず、変調された速波振動を構築するために、上記の如くω_ｆ＝２πｆ_ｆを中心とする速波振動を構築してσ_ｆに正規化してから、それを次の方程式により変調した。

式中、Ｋ_ｔ ^ｍｏｄ及び位相φ_ｔ ^ｍｏｄは時間的に変化し、図１８Ａと図１８Ｂに示す動的特性にめいめい従い得る。以下の代替変調機能を使用しても模擬データを生成できる。

式中、ｕ^Ｔ _φｍｏｄ＝［ｃｏｓ（φ^ｍｏｄ）−ｓｉｎ（φ^ｍｏｄ）］である。

急激又は著しい遷移を伴う信号は、アーティファクト的な位相振幅相互相関又は位相振幅変調につながる可能性がある。かかる条件の下における本発明に係る状態空間ＰＡＣ法のロバスト性を評価するために、ファン・デル・ポール振動子を使って急激な変化を伴う信号を生成した。ここで、振動χは微分方程式により規制される。

一定の時間ステップでオイラー法を採用することにより方程式（８９）の解を得た。

要するに、本発明に係るアルゴリズムは、変調に由来する非線形を第２段階で回帰モデルと当て嵌める前に、その第１段階で当該非線形を抽出できる。当方式の結果として生じる深刻な事態の主な例としては、速波成分の推定にけるバラツキが急激に増大することが挙げられる。速波振動の推定における幅広のＣＩの例については、図１２―２の（Ｇ）を参照されたい。次に、実際の場合よりも広範な分布から速波振動振幅を再サンプリングする。これはφ^ｍｏｄの推定に影響を及ぼさないが、Ｋ^ｍｏｄを再サンプリングする場合には控えめな推定がなされる、即ち、Ｋ^ｍｏｄを再サンプリングしない場合と比べて確信区間が広くなる。たとえそうであっても、依然として本発明に係る方式の方が従来の方法よりも優れた性能を発揮する、という結論に達した。

本発明は、振動の状態空間モデルを位相振幅結合（ＰＡＣ）のパラメトリック定式化と統合する新規な方法を提供する。本状態空間モデルでは、各振動を解析信号として表すことにより、位相又は振幅を直接推定する。その後、パラメトリックモデルを使用して、制約付き直線回帰によりＰＡＣ関係を特徴づける。回帰係数（結合関係をほんの少数のパラメータによって効率的に表現する回帰係数）を第２の状態空間モデルに組み込んで、ＰＡＣにおける時間依存性の変化を追跡できる。以上、この方法の効率の高さを、数々の用途において収集した神経系の時系列データを分析することより実証するとともに、異なる生成モデルに基づいたシミュレーション研究によると、統計的有効性が標準的な技法と比べて向上することを例示的に説明した。最後に、本発明に係る方法が如何にして、標準的な交差周波数結合分析法にまつわる多くの制限に対してもロバスト性を発揮するかを説明した。

本発明に係る方法の効率の高さは数々の要因に起因する。まず、状態空間解析信号モデルは、分析の対象となっている振動の位相と振幅への直接のアクセスを提供する。また、この線型モデルは、データから推定される「ソフトな」周波数帯域制限を課すことにより非線形特性（変調）を抽出できる卓越した能力を備えている。よって、振動分解は、厳格な帯域制限に束縛されない生理学的信号を分析するのに良く適している。また、本発明では、非線形振動を表現できる高調波拡張を提案し、それによって、帯域通過フィルタ処理によるアーティファクトの結果生じる不要波のＰＡＣと、真のＰＡＣとをより良く区別できるようにした。結合関係のパラメータ表現によって、異なる変調形状に対応できるようになり、モデル効率が更に向上する。

上記説明を通じて、本発明は、ＰＡＣ分析など、特定の種類の交差周波数分析に現在用いられている方法にまつわる大きな制限の大多数を取り上げた。神経系の時系列をより効率的に処理し、対象とする周波数帯域を自動的に選択、抽出、モデル化する。標準方法とは反対に、本発明は、時間を隔ててＰＡＣに関連する数量を平均する必要がないので、要求される連続時系列データの量を削減できる。なおまた、本発明で提案されているモデルでは、対象とする信号の事後分布を明確に定義できる。かかる事後分布から標本をサンプリングすることにより、（データ内の非定常構造を識別しづらくするとともに偽陽性率を過小評価する）代替データを構築する必要性を回避することができる。逆に、ＳＳＣＦＡ法では変調パラメータの事後分布を推定するので、ＣＩを報告し、本発明による結果の統計的有意性についての情報並びに変調の強さと方向を提供する。よって、本発明に係るＰＡＣの動的推定によって、徐波（０．１〜１Ｈｚ）信号の場合６秒という遥かに短い時間窓に基づいて推測することが可能になる。被駆動自己回帰モデル（ＤＡＲ）や一般化線形モデル（ＧＬＭ）など、ＰＡＣを表現するその他新規の方法も提案されている。先に述べたように、ＳＳＣＦＡ法は、特に信号対雑音比が低い場合に、ＤＡＲや標準の方式よりも優れた性能を発揮する。ＧＬＭは、本発明同様に変調関係をパラメータで表現するが、その際、統計的な面であまり効率的とは言えず、ブートストラップを使って信頼区間を提供するさらに一般的な形態で該パラメータ表現を行う。ＤＡＲとＧＬＭはいずれも、依然として信号の抽出を帯域通過フィルタに依存して行うため、やはり、かかるフィルタによりもたらされる重大問題に弱い。本発明に係る方法は、状態空間モデルを変調のパラメトリックモデルと組み合わせて使用する初の方法である。

かかる改良によって、将来ＣＦＣをめぐる研究を行う際に採用される分析を大幅に改善することが可能になり、ＣＦＣをほぼリアルタイムで追跡することが要求される医療用途も実施可能となり得る。かかる用途の一例としては、ＥＥＧに基づく麻酔誘発性無意識状態の監視が挙げられる。プロポフォール誘発性の麻酔状態において、周波数帯域を極めて明確に定義できるだけではなく、ＰＡＣシグネチャにより、深い無反応状態（山マックス）と遷移状態（ｔｈｒｏｕｇｈ−ｍａｘ）とをはっきりと区別することに加え、視床の極性レベルにおいて根本的な変化を反映できる可能性が最もあり得る。したがって、ＰＡＣにより、高感度且つ生理学的に妥当な麻酔誘発性脳状態マーカを提供し得るので、スペクトル特性だけでは提供できない更なる情報を提供できるようになる。しかるべく、ある研究によると、全身麻酔下にある患者で無意識状態をスペクトル特性だけでは完璧に予測できなかった事例が報告されている。かかる同様のデータにおいて、ＣＦＣの測定（山マックス）を使えば、患者が覚醒不能能な完全無意識状態をより正確に示すこと可能になる。手術室では、ボーラス注入法によって薬物を急速に投与する場合があるので、薬剤注入速度も急激に変化することがあり、患者が外科的刺激により覚醒し、それに応じて患者の脳状態が、時間とともに（秒レベルの時間尺度で）変化してしまう事態につながり得る。脳状態におけるかかる急激な遷移は、従来の方法を使用して推定された変調パターンではぼやけてしまう場合がある。このため、より高速且つより信頼性の高い変調分析は、全身麻酔管理におけるすさまじい改善につながり得る。従来の方法は、一つの推定値を生成するために数分間分のデータを必要とするので、非実用的である一方、本発明に係る方法は、かかる用途にふさわしい秒レベルの時間スケールでＣＦＣを推定できる。

当初ＣＦＣ分析法が神経科学に導入されて以来、健常者と疾病者のいずれにおいてもＣＦＣが脳の協調の基本的機構であることを示唆するデータが豊富に存在する。本発明に係る方法は、統計的効率と時間分解能を著しく向上するとともに、既存の技法にて遭遇する難関な問題の多くに対処する。このような性能の改善により、分析法が不十分であるために限界に達していた重要且つ新たな発見への道が切り開かれるとともに、ＰＡＣ分析の信頼性と効率の向上を通じて医療用途での活用が実現され得る。

以下、振動波と様々な周波数の範囲との間の交差周波数分析を実行する方法を説明する。一般的に、以下の方法では、フィルタリング技術を使用して徐波周波数と広帯域の周波数を分離するが、当然のことながら、上記の状態空間モデリング法は、フィルタリング技術の少なくとも一部の代わりとして使用することにより、振動をより好ましく当て嵌めるとともに、一可能性としてより、余波振動などの振動と５Ｈｚ〜５０Ｈｚ程度の広帯域の周波数範囲との間の関係について正確な推定を提供することもできる。以下、上記の状態空間モデル及び関連する当て嵌め技法を使用して、余波振動と広帯域波との間の関係を推測する方法についても説明する。徐波及び広帯域の波の関係を推定するためにどのように使用できるかについて説明するだろう。

過去数年にわたり、前頭皮質と後部皮質の意識状態媒介における役割についての議論がなされている。「後部皮質ホットゾーン」の仮説では、後部感覚皮質と感覚連合皮質が、（前頭皮質とは別に）意識状態を媒介する主な役割を果たしている、と提案されている。一方、前頭・後部相互作用が、意識的知覚表象を引き起こす重要な役割を果たしている、と主張する研究報告もある。近年、かかる議論が更に拡大し、無意識状態に関する様々な研究や、前頭部と後部における機能障害が如何にして意識の消失に関係しているのかについても話し合われている。

異なる無意識状態（麻酔状態、ＮＲＥＭ睡眠状態及び昏睡状態）は独特な電気生理学的シグネチャをしている。但し、これら全てに共通する一つの特徴としては、脳波記録（ＥＥＧ）において約１Ｈｚのあたりで交互の大きな偏向として認められる徐波活動が挙げられる。神経細胞は所定時間持続する脱分極期（活動期）と過分極期（静止期）の間を繰り返すが、上記徐波は、かかる根本的な皮質活動期と静止期の大縮尺の指標である、と考えられる。神経細胞は、脱分極された活動期には活性し、過分極化された静止期には不活性となる。

無意識状態の徐波活動は空間的に広範囲に分布するが、近年の研究において、徐波動的特性の空間的分布が無意識状態と如何にして関係するのかについて検討されている。ＮＲＥＭ睡眠中、後部の徐波パワーが 睡眠を夢見ることと比較して、夢見ない睡眠中に後部野に対する徐波パワーが、無夢睡眠時の方が夢見睡眠時よりも強くなる。対照的に、前頭部は、麻酔誘発性の無意識状態と密接に関係している。麻酔下にある動物の試験を行ったところ、前頭前野皮質に薬理学的刺激を与えると、徐波振動パワーの明らかな減少が認められるとともに、行動覚醒が回復される。麻酔誘発性の徐波振動は、麻酔深度に応じ、異なる位相（谷で最大となる「谷マックス」動的特性とピークで最大となる「山マックス」動的特性）にて前頭部アルファ振動を変調する。つまり、前頭部のアルファパワーが徐波のピークにて最大となる時（山マックス）、痛刺激によっても無意識状態からの覚醒が不可能な深い麻酔状態を反映すると考えられる。これらのデータは、徐波よりも周波数が高い律動に対する徐波の変調効果に基づいて、類似した徐波パワーを有する異なる行動状態をお互いに分離できることを示唆する。徐波パワーから徐波変調に視点を移すことによって、無意識状態の媒介における後部皮質と前頭皮質の役割を和解できるのではないだろうか。

前回の無意識時交差周波数変調分析はアルファ律動に対する徐波の影響に焦点を置いているが、活動期と静止期は律動性神経活動と非律動性神経活動の両方に影響するので、皮質活動期と静止期に基づいて徐波を解釈すると、如何なる周波数においてもＥＥＧパワーを余波にカプリングすべきである、ということを示唆する。この考え方に沿って、ＥＥＧで測定された幅広い周波数の範囲にわたって徐波の変調効果を追跡するために新規分析法を導入する。本発明では、プロポフォール誘発性無意識状態の異なる段階にわたって異なる皮質領域における余波変調を分析するためにこの方法を使用する。まず、前頭における山マックスの動的特性は広帯域の現象であることが分かり、山マックスの動的特性は根底にある皮質の活動状態と静止状態を反映することが示唆される。なおまた、異なる無意識状態（つまり、無意識状態が浅い時の後部野と無意識状態が深い時の前頭野）において、後部及び前頭野ともが徐波による広帯域変調を示すことも分かる。この結果は、麻酔誘発性の無意識状態が一種の状態に限られず、むしろ、めいめい異なる度合いの後部皮質と前頭皮質の妨害を示すとともにそれに対応して異なる意識処理及び行動覚醒能を示す複数の状態からなる、という考え方を支持する。

図１９Ａに示される如く、帯域通過フィルタ処理を施した余波電圧（０．１〜４Ｈｚ）と帯域通過フィルタ処理を施した高周波信号の即時振幅との間の相関関係を使用してクロス周波数結合を定量化した。交差周波数結合を一次元的に要約して示す。つまり、正の相関関係にある場合は、徐波電圧が正電圧（ピーク最大結合山マックス）の時に高周波振幅が高くなり、負の相関関係にある場合は、徐波電圧が負電圧（谷最大結合谷マックス）の時に高周波振幅が高くなることを意味する。その後、本発明は、異なる周波数が徐波とはどのように関係するか見るために高周波振幅に対して使用された周波数帯域を変えることができる。この方式を採用すると、ピーク最大結合山マックスと谷最大結合谷マックスの動的特性が、無意識状態への遷移中に電極や周波数によって如何に変化するかを見ることができる。

図１９Ｂは、代表的な被験者において山マックスと谷マックスの動的特性がどのような挙動を示すのかを示す。この際、該被験者が、ボタンを押すことによって反応をしめす聴覚課題に取り組んでいる最中に、プロポフォールの投与量を１４分ごとに増量しながら投与した。プロポフォール投与量が十分増量された時点で、刺激に対する被験者の応答が途絶えた（意識消失状態：ＬＯＣ）。投与量を５段階に分けて増量した後、均一に減量を開始してから間もなく、被験者が刺激に対して再び応答を示した（意識回復状態：ＲＯＣ）。

前頭電極での相関関係に基づく交差周波数結合によって、先行する研究により得られたアルファ律動に関する基本的な結果を再確認できた。つまり、１０Ｈｚにて、徐波への結合が谷マックス動的特性として始まり、投与量が高くなった時点で山マックスに切り替わり、そして意識回復中に再び谷マックスに戻る。しかし、１０Ｈｚ以上を見てみると、約２０Ｈｚにおよぶ高い周波数（プロポフォールで鎮静状態を誘発した際に被験者が意識を保っている状態に対応する高いパワーを示す周波数）にてＬＯＣよりもはるか前に谷マックス動的特性が現れる。これは、谷マックス動的特性が示される間は意識がありえる、という洞察とも一貫している。

恐らく、交差周波数結合における最も顕著な結果は、山マックスパターンが全周波数（５〜５０Ｈｚ）に広がるということである。この結果は、アルファ波と余波以外の信号に狭帯域律動がないので、非律動性の広帯域活動が余波に結合していることを示唆する。これは、山マックス動的特性が根本的な皮質の活動期と静止期の挙動（皮質活動が、律動性又は非律動性にかかわらず、徐波の谷部で中断されること）を反映する、という解釈と一貫している。

前頭電極と後部電極との間の交差周波数結合を比較すると、意識消失後まもなく、前頭の谷マックスと同時に、後部チャンネルにおいて広帯域の山マックスパターンが生じることが分かる。即ち、前頭のチャンネルが山マックス動的特性に遷移する前に、後部のチャンネルは既に徐波への広帯域結合を示す。全ての被験者において、後部山マックス動的特性の発現が、ＬＯＣ後、前頭の山マックス動的特性の発現前に認められる。かかる結果は更に、８〜１６Ｈｚにおける交差周波数結合を濃度別に示すトポグラフィックマップにも示される。つまり、意識消失後まもなく、前頭電極群が依然として谷マックス結合を示す一方で、後部野の電極は山マックス結合を示す。プロポフォール濃度の上昇につれて、前頭電極は山マックス結合に関与し始める。このことは、皮質活動が余波により広帯域で静止する現象は、プロポフォールの高投与量で前頭に広がる前に、まず後部野から始まることを示唆する。

以降、濃度を次の４段階に分けて分析を行った。
（１） ベースライン：プロポフォール投与前
（２） 鎮静 ：被験者の反応が途絶える前の濃度
（３） 低無意識用量：被験者の反応が途絶えた後の第１濃度
（４） 高無意識用量：所望の無意識状態を得るために要される以上の量であるが、バーストサプレッションや昏睡のような状態をもたらさずに投与できるプロポフォールの最大投与量
これによって、ＬＯＣ及びＲＯＣの時間的経過が被験者ごとに異なると仮定した場合に、当該各被験者の情報を組み合わせることが可能になる。

周波数にまたがる交差周波数結合パターンが電極やプロポフォール濃度及び被験者によって異なるので、この活動を最も良く要約できる結合パターンの特定を試みた。そのために、非心性の主成分分析（ＰＣＡ）を使用して、図２０Ａに示すように交差周波数結合パターンを基本周波数モードに分解した。該主な基本モードは、周波数にまたがる直交パターンであり、交差周波数結合を表す行列において総エネルギーの徐々に少なくなる割合を、周波数、電極、プロポフォール濃度と被験者別に捉えるものである。

第１の基本モードは、全周波数に対して正であり、したがって、交差周波数結合における広帯域の影響を表すとともに、総エネルギーの７８％を捉える。第２と第３の基本モードは、めいめいアルファ及びベータ周波数帯域に狭帯域ピークを有し、第１の基本モードよりもかなり少ないエネルギーを捉える。よって、広帯域皮質活動を徐波に同調することが、プロポフォール誘発性麻酔下における交差周波数結合に最も大きく寄与する。各モードの総エネルギー割合（％）を図２０Ｂに示す。

この第１主要周波数モードは、個々の被験者の要約情報に認められる広帯域の山マックスパターンを捉えることができるので特に関心を引く。特に、該第１主要周波数モードは、プロポフォールの投与量を次第に増やしながら投与する際に、皮質のＥＥＧ発生源（ジェネレータ）における広帯域山マックス現象の空間的分布を説明するために使用することが好ましい。そのために、音源定位されているＥＥＧ信号について交差周波数結合分析を行い、皮質個所の全体にわたって結果として得られる交差周波数結合パターンを、上述のセンサ空間分析から得た第１の基本モードに射影した。図２１は、結果として得られた投影（フリーサーファー：ＦｒｅｅＳｕｒｆｅｒ）平均表面にモーフィングして被験者の全体にわたって平均したもの）を示す。かかる表面において、正の投影（赤色）は、徐波に対する広帯域の山マックス結合を表す。

個々の被験者についてのセンサ空間結果と一貫して、低無意識用量条件下においては、徐波に対する広帯域山マックス結合が後部皮質領域に存在し、高無意識用量条件下において強化されて、前頭皮質領域に広まる。異なる状態における脳の異なる個所での低周波活動による広帯域周波数変調を示す図２１のグラフから、低無意識用量条件下では頭頂葉、側頭葉と後頭葉がいずれも広帯域山マックス結合を有するが、前頭葉は徐波に対する広帯域結合を実質的に持たないことが分かる。対照的に、高無意識用量条件下では、前頭葉が広帯域山マックス結合においてその他の葉と結合される。

これらの結果はともに、図２０からも明らかなように、徐波に対する広帯域山マックス結合がプロポフォール誘発性麻酔下において交差周波数結合動的特性に大きく貢献することを示唆する。更に、図２１と図２２からも明らかなように、意識消失後、広帯域山マックス結合がまずは後部皮質領域で認められ、その次にプロポフォールに投与量が高い時に前頭野で認められる。

ここでは、徐波に対する広帯域結合を、ＥＥＧ信号に見られる局所的な皮質の活動期と静止期のマクロスケールの指標として説明する。集団スパイキング活動は、神経細胞パワースペクトルに対して広帯域作用を有するので、皮質の活動期と静止期に集団スパイキング活動を余波に同調することにより、広帯域パワーが徐波の特定の位相において最も強くなるはずである。この考えは、活動期が集団スパイキング活動と広帯域パワー（＜５０Ｈｚ）が同時に増加する現象に関連することを示す脳波測定及び局所電場電位データにおけるマイクロスケールの証拠により立証されている。ＥＥＧにおける余波パワー増大が常に皮質興奮の活動期と静止期に対応する、と一般的に想定されているが、本発明者らが収集したデータによると、徐波は、広帯域変調を生じるのに十分な集団スパイキング活動の同調前に、臨界振幅（場所によって異なり得る）に達しなければならないことが示唆される。後部野では、前頭野よりも低いプロポフォール投与量で閾値を越える。この余波に関する臨界閾値の考え方は、徐波パワーの飽和に伴って視床と感覚皮質との間の切断が生じるという徐波活動飽和の仮説と一貫している。本発明者らは、皮質の活動期と静止期自体がかかる感覚隔離の仕組みを提供すると仮定している。この解釈によると、広帯域山マックスの動的特性は、局所的な皮質の活動が活動期と静止期によりＥＥＧで検知し得る十分に大きな規模で妨害されていることを表す。

全身麻酔下にある患者の無意識状態を監視できれば、明らかな実益をもたらすことができるが、それを実現できる理に適った方法を特定するといった課題は未解決のままである。本発明では、外部刺激によって患者の意識を回復できる（覚醒可能な）無意識状態と、患者が外部刺激によっても目を覚まさない（覚醒不能の）無意識状態とを区別できる。最近のデータの中には、徐波の谷部で前頭のアルファ律動が最も強い（谷マックス）時には患者を覚醒可能であるが、徐波のピークで前頭のアルファ律動が最も強い（山マックス）時には覚醒不能のであることを間接的に示すものもある。ここで、前頭のアルファ律動が徐波のピークに結合している（山マックス）時、前頭及び後部の皮質の広帯域スペクトルパワーも徐波のピークに結合されることを説明する。よって、前頭及び後部の皮質の余波に対する広帯域結合に関連して覚醒不能状態になることが示唆される。麻酔量が低く、患者は無意識であるが恐らく覚醒可能な状態において、後部皮質が広帯域山マックス動的特性に関与する一方で、前頭皮質が（アルファ帯域の）谷マックス動的特性に関与するのが認められる。したがって、これら交差周波数結合の異なる空間的パターンが、異なる無意識状態を示すと考えられる。麻酔専門医が、これらの状態に関する知識を活用することにより、覚醒可能又は覚醒不能の異なる無意識状態を、各臨床状況に応じて確立することが可能になる。

「後部皮質ホットゾーン」の仮設及び前頭−後部の仮説によると、後部回路と前頭回路の意識状態の媒介における役割が明確に異なる、とされている。麻酔はこの問題に新たな次元をもたらす。何故なら、麻酔剤の投与量によっては、無意識の患者を外部刺激により覚醒可能な状態又は覚醒不能の状態にし得るからである。この次元に沿って考えると、前頭前野皮質が、前頭前野皮質への刺激による意識活性化における中心的な役割を果たすと考えられる。本発明らが得た結果は、後部皮質と対比して前頭皮質の妨害が異なる無意識状態と同時に起こるという点で上記の見解を調和させるように思われる。第一に、後部皮質ホットゾーンの仮設と前頭−後部の仮説はいずれも、後部野の妨害によって意識の内容が劣化するはずであると予測する。本発明らが得た結果によると、プロポフォールの投与量が低い（被験者が無意識であるが覚醒可能であり得る）時にまずは後部野が最初に広帯域山マックスを示す。この時の後部野は、非夢見睡眠時の徐波パワーが高い時の状態に非常によく似ている。このことは、後部皮質領域が、無意識状態中に活動期と静止期を交互に切り替えることにより妨害を受けた意識の内容を媒介する、という解釈と一貫している。第二に、前頭野では、プロポフォールの投与量が多い（覚醒不能状態に関連する条件）時に、広帯域山マックス結合が生じる。このことは、前頭前野皮質が意識の回復に伴う行動覚醒の不可欠であり、前頭前野皮質における活動期と静止期を交互切り替えによってかかる覚醒を防止できる、という解釈と一貫している。

本研究の結果で重要な点は、脳全体の交差周波数結合分析によって、無意識状態への遷移中及びプロポフォール誘発性麻酔下の覚醒不能状態において、空間に依存して変化する脳の動的特性を確実に捉えることができるということである。本発明では、皮質機能を妨害する機能を果たす、マクロスケールの皮質活動における活動期と静止期のマクロスケールの指標として徐波ピークに対する広帯域結合を解釈する。この機能的妨害は、異なる無意識状態において前頭及び後部皮質領域の両方を網羅する。本発明らが得た結果によると、麻酔下の無意識状態は、特異な現象ではなく、むしろ脳状態におけるいくつかの特有な変化に関与しており、後部皮質による意識の「内容」の処理（前頭葉前部の覚醒とは異なる）と及びそれらの内容に関する認識度とを妨害することが示唆される。

前記データを分析するために研究を実施した。１８歳〜３６歳の健常者ボランティア１０名が本研究に参加した。研究前に、各被験者の構造ＭＲＩ（シーメンス・トリオ３テスラ、Ｔ１強調磁化準備形高速グラディエントエコー、スライス厚１．３−ｍｍ、面内分解能１．３ １ｍｍ、ＴＲ／ＴＥ＝２５３０／３．３ミリ秒、フリップ角＝７）を得た。フリーサーファー画像解析スイートを使って皮質の復元を行った。

少なくとも１４分間目を閉じた状態で休息した後の被験者に、プロポフォール麻酔剤を投与し、その際、プロポフォールの投与量を１μｇ／ｍＬ、２μｇ／ｍＬ、３μｇ／ｍＬ、４μｇ／ｍＬ、又は５μｇ／ｍＬを標的効果部位濃度として１４分ごとに増加した。意識回復時には、プロポフォールの投与量を、１４分ごとに効果部位濃度０．５μｇ／ｍＬ、１．０μｇ／ｍＬ、１．５μｇ／ｍＬと３段階に分けて、聴覚合図に対する被験者の反応が途絶えた投与量未満に減少した。高密度（６４チャンネル）ブライアンビジョンＭＲＩプラスシステム（ブライアン・プロダクツ（Ｂｒａｉｎ Ｐｒｏｄｕｃｔｓ）製）を使って、毎秒標本５，０００個のサンプリングレートで、脳波（ＥＥＧ）記録データを収集し、電極位置をディジタル化した（ポルヘムス・ファーストラック３Ｄ）。

研究を通じて、被験者は、クリック音や言葉による指示などの聴覚的合図に対してボタンを指で押すことによって応答する課題に携わった。状態空間モデルで正しい応答の確率が５％未満と推定され、少なくとも５分間その状態を維持した時に意識消失（ＬＯＣ）とした。同様に、状態空間モデルで意識回復中に正しい応答の確率が初めて５％に達し、少なくとも５分間その状態を維持した時に意識回復（ＲＯＣ）とした。

ＥＥＧ前処理：
重度のアーティファクトが持続的に認められるチャンネルを目視検査で特定して分析対象から取り外いた。また、ｅブリッジ・ソフトウェアパッケージによって電気的にブリッジされていると思われる電極も取り除いた。

後述される交差周波数結合分析には、徐波帯域（０．１〜４Ｈｚ）で帯域通過フィルタ処理されたバージョンの信号と、１Ｈｚの遷移帯域を有する帯域幅における一式の高帯域（４〜５０Ｈｚ、それぞれ２Ｈｚ）を使用する。零位相ＦＩＲフィルタを使用してこれらの帯域に全セッション分のデータをフィルタ処理し、毎秒標本２００個の速度でダウンサンプリングした。ラインノイズを避けるとともに、最大ＥＥＧパワーを有する周波数を捉えるために高帯域の上限を５０Ｈｚに設定した。

センサレベル分析については、本発明は、フィルタ処理後最初の最近傍に基づいてラプラスの基準を適用した。

プロポフォールを段階的に投与したため、１回分の量を一定の予測効果部位濃度で被験者に投与する時間は約１２分であった。ここでは、かかる期間を「レベル」と呼び、関心の対象となるベースラインレベルと、鎮静レベルと、低無意識用量レベルと、高無意識用量レベルの４つのレベルを定義した。ベースラインを、プロポフォール投与開始前の期間として定義し、鎮静を、ＬＯＣ前のレベルとして定義し、低無意識用量を、ＬＯＣ後の初めて総量が投与されるレベルとして定義し、そして高無意識用量を、プロポフォールの可能最高量が投与されるレベル又は（バーストサプレッションを呈した被験者場合）バーストサプレッション前の段階として定義した。これらのレベルについて行なわれた分析には、アーティファクトを伴わない長さ３０秒のエポックを、各被験者について各レベルから１０個選択した。被験者のうちの２名は、プロポフォールの第１レベル中に反応が途絶えた。したがって、当該２名については鎮静期間がなかった。被験者間の結果をまとめるデータによると、ベースライン、低用量及び高用量レベルの被験者が１０名、及び鎮静レベルの被験者が８名という結果になった。

交差周波数結合分析：
プロポフォール投与下において、高周波振幅が徐波のピーク又は谷に結合する可能性の方が、上昇相又は下降相に結合する可能性よりも高い。故に、正（山マックス：高周波振幅が徐波のピークで更に高くなる）又は負（谷マックス：高周波振幅が徐波の谷で更に高くなる）のメトリクスを使用して交差周波数結合をまとめる。特に、帯域通過した高周波信号の即時振幅と帯域通過した余波の電圧（０．１〜４Ｈｚ）との間のピアソンの相関係数を使用する。

式中、Ｖは、徐波電圧の時系列を表す列ベクトルであり、Ａは、高周波帯の即時振幅の対応する時系列を表すベクトルである。即時振幅は、帯域通過フィルタ処理された信号の解析信号の大きさに基づいて計算した。相関関係を計算する前に、３０秒エポックごとの即時振幅を、平均を差し引くことにより中心化した。Ｖの期待値が偽（ゼロ）であるので、Ｖを明確に中心化できなかった。

空間（電極、音源箇所）且つ又は時間（複数の３０秒エポック）にわたって平均相関関係を見つけるために、平均すべき次元を方程式１のＶ及びＡの列ベクトルにおいて縦に積み重ねる。該相関関係が、各信号（分母に含まれる、√（Ｖ^ＴＶ）及び√（Ａ^ＴＡ））の標準偏差により正規化された信号（分子に含まれる、Ｖ^ＴＡ）同士の共分散を表すので、信号を縦に積み重ねることにより、正常化を行なう前に共分散及び標準偏差を別々に推定できるという効果を奏する。

ちなみに、この分析には、位相振幅結合分析に通常含まれる徐波の位相を推定する過程がないことに留意されたい。ここでは、以下３つの理由で位相の推定を行わないことにした。第１の理由は、０周辺（徐波ピーク、正の信号の場合）とπ周辺（徐波谷部、負の信号の場合）のプロポフォール・クラスタにおける位相振幅結合である。第２に、徐波は、正弦波により明確に記述されているか否かを問わず活動期と静止期を反映すると考えられている。多くの場合、プロポフォールに誘発された無意識状態の初期において、徐波は、波形形状と周波数の面でかなり不規則である。ここでは、波形の非正弦波特性をさらに多く捉えるために、徐波（０．１〜４Ｈｚ）に帯域幅の広域を使用することにした。プロポフォール誘発性麻酔下では、狭帯域デルタ周波数の律動的活動（１〜４Ｈｚ）は存在しないので、この周波数帯域の帯域通過信号に対する寄与が非律動性であり、波形の形状に影響を与える。したがって、相関関係を使って高周波振幅の非正弦波低周波波形に対する関係を記述することにより、広帯域信号に明らか認められる活動期と静止期の変動と高周波活動との間の関係をより良く捉えることが可能になる（図１９Ａ）。最後に、問題を一次元化すること（正と負の結合）によって単一の個所と周波数の分析を簡易化できるので、分析をその他の周波数と個所に展開して行えるようになる。

図１９Ａは、２つの時間窓における余波の電圧（低周波の活動、ＬＦＡ）とアルファ帯域（α、９〜１１Ｈｚ）との相関関係を計算することにより行う前頭電極の計算を示す。左側の時間窓では、余波の電圧が負の時にアルファ振幅が高くなり、負の相関（谷マックス）になる。右側の時間窓では、余波の電圧が正の時にアルファ振幅が高くなり、正の相関（山マックス）になる。このような分析は、セッションに含まれる全ての３０秒エポック及び様々な振幅周波数（４〜５０Ｈｚの間の２Ｈｚ帯域）について行なうことができる。これにより、選択された電極について、当該セッションにおける時間を隔てた各周波数と徐波との間の交差周波数結合を示すモジュログラムを得ることができる。同様に、プロポフォールレベル内の各電極について相関関係を計算し、プロポフォールレベルの関数としてアルファ帯域と徐波の間の交差周波数結合の空間的分布を示すことができる。

振幅帯域、時間、電極及び被験者にわたって交差周波数結合分析を行なった。ここでは、時間を２つの方法で特徴づける。第１の方法は、セッション・ベースの分析であり、セッションに含まれる３０秒エポック毎に結合を別々に評価し、結果として、前頭・後部要約情報を示す図１９Ｂのモジュログラムを得ることができる。これらのモジュログラムは、右側の頭皮マップに示されるようないくつかの電極間の平均相関関係を示す。時間を特徴づける第２の方法は、プロポフォールレベルに基づいており、各プロポフォールレベルにおけるアーティファクトを伴わない３０秒エポック１０個について結合を評価する方法である。故に、それぞれの結合値は３００秒分のデータについての相関関係を表し、各電極とレベル及び被験者についてかかる値が存在する。このようなレベルベースの分析は、図１９Ｂのトポプロットに用いられており、後述の主成分分析にも採用される。

音源定位：
ＭＮＥパイソン（ＭＮＥ−ｐｙｔｈｏｎ）ツールボックス（２７．ｍａｒｔｉｎｏｓ．ｏｒｇ／ｍｎｅ／）を使って最小ノルムを推定することにより音源定位を行った。ノイズ共分散を対角と想定し、ＥＥＧをもとに６５Ｈｚ以上のスペクトルパワーを使って推定した。音源空間は、フリーサーファーによる皮質表面再構築によるとｉｃｏ−３デシメーション（間引き）であり、皮質表面に垂直になるように音源を制約した。３層境界要素モデルを使用して、頭脳内面から５ｍｍ以内に存在する音源をすべて取り除くようにフォワードモデルを計算した。

音源空間において交差周波数結合分析を行なうために、上述の方法を使用して、徐波帯域と４〜５０Ｈｚ間の２Ｈｚ帯域のそれぞれにおける帯域通過信号の音源定位を行った。皮質表面を示す図（図２１）においては、上述の交差周波数結合分析を、上記４段階のレベルについて各被験者における各音源箇所に適用して行った。各葉における平均交差周波数結合（図２２）にも同様の技法を使用したが、各葉内の全ての音源の信号を組み合わせた。つまり、該各葉における全ての音源について、（各レベルの）１０個のエポックを方程式１のＶ及びＡの列ベクトルにおいて縦に積み重ねた。

非心性主成分分析及び音源空間への射影：
観測された特徴の大部分を説明する周波数にわたる結合のパターンを特定するために、センサ空間結合パターンについてプロポフォールレベル間及び被験者間の非心性主成分分析を行なった。４つの対象レベル（ベースライン、鎮静、低無意識用量と高無意識用量）に関するレベル別結合結果を集計行列Ａ（ｆ，ｉ）に組み合わせた。該集計行列Ａ（ｆ，ｉ）において、ｆは振幅帯域の指標となり、ｉはプロポフォールレベルと電極と被験者（即ち、積み重ねられた全てのレベルと電極と被験者を含む行）の指標となる。その後、全被験者に関し、上記４つの対象レベルで２個の電極について行われた交差周波数結合の結果を生成した。当該分析は、これらのパターンを、周波数にわたり当該パターンの重要な特徴を捉える成分に分解することを目的とする。主成分分析には特異値分解が関与している。

上記分解において、Ｕは、列を基本モードとする行列である。１番目の列は第１の基本モード、つまり、Ａ行列におけるほとんどのエネルギーを捉える周波数全体に及ぶパターンである。２番目の列は第２の基本モードであり、第１モードが取り除かれた後にデータに含まれるほとんどのエネルギーを捉える。基本モードはいずれも直交式であり、単位長さを有する。正値の山マックスへのマッピングと負値の谷マックスへのマッピングを維持するために、必要に応じて基本モードを反転して（−１を乗じて）、最大の要素が正になるようにした。Ｓ行列は対角であり、ｊ番目のモードにより示される総エネルギーの割合（％）を推定するために使用することができる。

ちなみに、ここでは、分散を分解する中心つきＰＣＡではなく、総エネルギーを分解する非心ＰＣＡを使用することにする。特に、原点に特殊な有意性があり、平均値を差し引きことによりデータが中心化されると該有意性が失われてしまうような状況において、非心ＰＣＡが有用である。交差周波数結合のパターンについては、原点が、ＥＥＧ信号が如何なる頻度においても余波への結合を含んでいないような状況を示す。対照的に、平均値は、（振幅周波数の関数として）被験者間、プロポフォールレベル間及び電極間の余波に対する平均結合を表す。原点は平均値よりも有意であり、よって、より一般的な中心つきＰＣＡよりも非心ＰＣＡの方が解釈し易い、と考える。

基本モードが音源空間において如何にして表現されているかを定量化するために、各被験者に関する音源空間パターンを第１のセンサ空間基本モード上に射影した。まず、音源空間分析から得た結合結果を集計行列Ａ^{ｓｏｕｒｃｅ}（ｆ，ｉ）に組み合わせる。ここで、ｉは、プロポフォールレベル、音源箇所（又は葉）と被験者の指標となる。その後、Ａ^{ｓｏｕｒｃｅ}のＵに対する射影は次のとおりになる。

結果として得られるＰ行列の最初の行は、（各プロポフォールレベル、音源箇所と被験者の）音源空間結合パターンの第１の基本モードへの投影を含む。第１の基本モードは周波数間で比較的一定であるので（図２０を参照）、このモードに対する投影は、広帯域である交差周波数結合を示す。また、第１の基本モードが全周波数に対して正であるので、正の投影は広帯域山マックス結合（徐波のピークに結合された全ての周波数）を反映し、負の投影は広帯域谷マックス結合（徐波の谷に結合された全ての周波数）を反映する。

（被験者×レベル×葉）結合パターンの第１基本モードへの投影を被験者間で平均し、当該葉における交差周波数結合に対する第１モードの平均的な寄与を表す図２２に示す推定を求めた。被験者間でブートストラップを使用して９５％の信頼区間を得た。

（被験者×レベル×音源箇所）結合パターンの第１基本モードへの投影を、ＭＮＥパイソンを使ってフリーサーファー平均表面にモーフィングし、結果として得られるマップを被験者間で平均した。結果として得られるマップは、各皮質個所における、選択されたモードの交差周波数結合への寄与を推定する（図２１を参照）。

上記のように、余波成分と様々な広帯域周波数との間の交差周波数関係を検出する際に状態空間モデルを使用することができる。状態空間モデルに基づく交差周波数結合分析法において状態空間モデルを使用するために、方程式（３）と（４）で次のようなパラメータ置換を行い、方程式（３）と（４）を以下のように置き換えることにより、余波成分だけを得ることができる。以下、参照のために方程式（３）と（４）を再び記載する。

式中、Ｍ＝［１ ０］、ｘ_ｔ＝ｘ_ｔ ^Ｓ、Ｑ＝Ｑ_Ｓ、Φ＝ａ_ＳＲ（ω_Ｓ）である。

その後、上記の如くＥＭアルゴリズムを使用してモデルを当て嵌めることができる。次いで、ベクトル時系列ｘ_ｔ ^Ｓの２つの要素（めいめいの実部と虚部）を、上述の交差周波数結合推定方法において、帯域通過フィルタ処理済み余波電圧（即ち余波信号）として使用することができる。高周波帯内の振幅も依然として、帯域通過フィルタ処理を使用し、解析信号の大きさを捉えることにより（即ち、ヒルベルト変換法の適用により）計算することができる。

その後、方程式（９０）を使用して、交差周波数結合を定量することができる。以下、参照のために方程式（９０）を再び記載しておく。

式中、Ｖは、ｘ_ｔ ^Ｓの第１又は第２の要素のいずれかを含む列ベクトルであり、Ａは、それに対応する、高周波帯の即時振幅の時系列を表すベクトルである。相関関係を計算する前に、即時振幅を、平均値を差し引くことにより中心化すべきである。Ｖの期待値は偽（ゼロ）のであるので、Ｖを必ずしも明確に中心化する必要はない。

その後、ｘ_ｔ ^Ｓの各要素について、図１９Ｂに似たグラフを生成できる。ｘ_ｔ ^Ｓの第１要素（即ち実部）に対応する分析によって、徐波のピーク（正、山マックス）と谷（負、谷マックス）に対する高周波活動の結合を捉えることになる。また、パターンは、時間と周波数と電極個所にわたる元の分析に定性的に非常に類似するはずである。ｘ_ｔ ^Ｓの第２要素（即ち虚部）に対応する分析によって、徐波の下降相と上昇相に対する高周波活動の結合を捉えることになる。一部の実施形態において、余波の振動を当て嵌めた後に残る残差を、帯域通過させて、帯域通過された元のＥＥＧ信号の代わりに使用することができる。残差の使用によって、推定される徐波の高周波要素をすべて広帯域信号から取り除き得る。

上述の状態空間モデルベース交差周波数結合分析法について考えられる問題の一つとしては次のことが挙げられる。つまり、余波範囲外の強い振動又は局所的なスペクトルパワーが信号に含まれていると、ＥＭアルゴリズムは所望の性能を発揮できない場合がある。例えば、徐波に関連する範囲外に（約１Ｈｚを下回る範囲に）周波数ｆ_Ｓがズレたり、且つ又は分散σ_Ｓ ^２が増加したりして、高周波パワーと徐波の両方を捉える可能性がある。

上記問題に対する解決法の一つとして、状態空間モデルのための振動成分の数を選択して、結果として得られる隠れ状態を使用することができる。その際、当該の隠れ状態の一つは徐波に対応するべきであり、後の残りは高周波数振動に対応する。その後、高周波成分の即時振幅と徐波成分の実／虚部との間の相関関係を計算できる。次いで、（信号の「非振動」成分反映する）モデルの残差を決定でき、ヒルベルト変換法を実行することにより（上記の如く、一式の周波数帯域にける）高周波残差の即時振幅を得ることができる。

解決法の一つとしては、高周波信号の即時振幅を計算するためにヒルベルト変換を使用する代わりに、高周波数をタイル表示するところの一式の高周波成分の使用する状態空間モデルからも即時振幅を得ることができる。つまり、正準状態空間方程式（３と４）を、徐波成分ｘ_ｔ ^Ｓと一式の高周波隠れ状態ｘ_ｔ ^ｊ（ｉは、例えば５〜５０Ｈｚの周波数の指標となる）で当て嵌めることができる。この場合、一般的に各々全ての高周波成分に狭帯域振動が存在するわけではないので、ＥＭアルゴリズムが思いがけない挙動を示すことがある。したがって、この方法が正常に機能するためには、ユーザが、周波数ｆ_ｉ（周波数帯域のピークの個所）と自己回帰係数ａ_ｊ（ピークの鋭さ）を、高周波帯のために、且つＥＭアルゴリズムによる更新を伴わずに、先験的に選択しなければならない。

図２と図１１ならびに図２３を参照すると、当て嵌められて状態空間モデルを使用する交差周波数分析を行うプロセス２３００の代表例が示されている。ステップ２３００は、少なくとも一つのメモリに関する指示として実施できる。その場合、かかる指示を少なくとも一つのプロセッサにより実行できる。

ステップ２３００のステップ２３０４では、ＥＥＧデータを受け取ることができる。ＥＥＧデータは、図２を参照して上述された患者監視装置２０２のような患者監視装置から受け取れる。ＥＥＧデータは、波形として解釈することができる。その後、ステップ２３００はステップ２３０８に移ることができる。

ステップ２３００のステップ２３０８において、状態空間モデルをＥＥＧデータに当て嵌めることができる。ステップ２３０８には、少なくともステップ１１００のステップ（例えばステップ１１０８〜１１２８）のサブセットが含まれている。その後、ステップ２３００では、それ以降のステップにおいて、当て嵌めモデルの使用により交差周波数解析を実行する。その後、ステップ２３００はステップ２３１２に移ることができる。

ステップ２３００のステップ２３１２において、当て嵌めモデルの徐波成分とアルファ波成分に関する交差周波数解析を実行することができる。ステップ２３００では、２つの波成分間の位相振幅結合のパラメトリック表示を一つ又は複数計算することができる。ステップ２３００では、方程式（３３）を使ってＳＳＣＦＡを計算できる。あるいは又は更に、ステップ２３００においては、方程式（４９）を使ってｄＳＳＣＦＡを計算することができる。ＳＳＣＦＡとｄＳＳＣＦＡはいずれも、アルファ波成分と徐波成分との間の位相振幅結合の正確な推定、より具体的には徐波成分の位相が如何にして、時間のビン化に依存せずにアルファ波成分の振幅を変調するのかについての正確な推定を提供できる。上記の如く、ｄＳＳＣＦＡはＥＥＧデータの時間窓にわたって時間的連続性制約を変調パラメータに課す。つまり、ｄＳＳＣＦＡは、時間的な制約に基づいて計算される。ＳＳＣＦＡ及びｄＳＳＣＦＡから得られるパラメトリック表明（複数可）ＰＡＣにより、徐波成分と速波成分との間の山マックス且つ又は谷マックス結合の指標を得ることができる。その後、ステップ２３００はステップ２３１６に移ることができる。

ステップ２３００のステップ２３１６において、徐波成分とＥＥＧデータの周波数範囲との間の交差周波数分析を実行できる。例えば、徐波成分と周波数５Ｈｚ〜５０Ｈｚの関係を推定し得る。徐波成分とＥＥＧデータの周波数範囲との相関関係は、方程式（９０）を使用して計算できる。方程式（９０）において、Ｖは、当て嵌めモデルから得たｘ_ｔ ^Ｓの第１又は第２の要素のいずれかを含むベクトルである。Ａは、それに対応する高周波帯（即ち５Ｈｚ〜５０Ｈｚ）の即時振幅の時系列を表すベクトルである。高周波帯を通過するように構成された帯域通過フィルタをＥＥＧデータに適用することにより、Ａを決定できる。余波成分をＥＥＧデータに当て嵌めた後に残る残差を、高周波帯を通過するように構成された帯域通過フィルタを帯域通過させることができる。該残差を使ってＡを決定することによって、推定される徐波の高周波要素をすべて広帯域信号から取り除くことができ、可能性として、より良い信号を提供できる。その後、ステップ２３００では、方程式（９０）を使って徐波と高周波帯の間の相関関係を推定できる。その後、ステップ２３００はステップ２３２０に移ることができる。

ステップ２３２０において、ＳＳＣＦＡ又はｄＳＳＣＦＡ、若しくは徐波と高周波帯モデルの間の相関関係をユーザが見れるように表示部に出力、且つ又は別のステップにおける保管且つ又は使用を目的としてメモリに出力できる。その他のステップで、計算により得られた交差周波数解析値（即ちＳＳＣＦＡ、ｄＳＳＣＦＡ、又は徐波と高周波帯の間の相関関係）をさらに処理し、使用者が各数値を簡単に理解できるような形式（例えば、各電極の個所に基づいて、複数の電極測定値から計算した複数のＳＳＣＦＡ値を脳のモデルにマッピングするなどの方法）で報告を作成しすることも可能である。かかる報告を、ユーザが見れるように表示部に出力、且つ又は別のステップにおける保管且つ又は使用を目的としてメモリに出力できる。その後、ステップ２３００が終了する。

上記において様々な構成について説明したが、かかる構成は単なる一例にすぎず、本開示の範囲を限定することを意図したものではない。本願に記載される構成を、本発明の要旨を逸脱しない範囲において様々な態様で実施できることは、通常の技術を有する当業者にとって自明なことである。特に、上述される一つ又は複数の構成から所与の特徴を選択して、上記において明示的に記載されていない特徴の部分的組み合わせからなる代替構成を生み出し得る。更に、上述される一つ又は複数の構成から所与の特徴を選択して組み合わすことにより、上記において明示的に記載されていない特徴の組み合わせからなる代替構成を生み出し得る。当業者にとっては、全体として本願を精査することによって、かかる組み合わせや部分的組合せに適した特徴を容易に理解できるであろう。本明細書及び特許請求の範囲に記載される本願発明の主題は、該当する如何なる技術の変更も網羅及び包含するものである。

Claims (20)

1. 患者の脳状態における少なくとも一つの振動成分を推定する方法であって、
   脳波記録（ＥＥＧ）データを受け取ることと、
   少なくとも一つの振動成分を含む状態空間モデルをＥＥＧデータに当て嵌めることと、
   前記状態空間モデルに基づいて前記ＥＥＧデータの交差周波数結合の少なくとも一つの値を推定することと、
   前記交差周波数結合の少なくとも一つの値に基づいて報告を生成することと、
   表示部とメモリの少なくとも一方に前記報告を出力することと、を含む方法。
2. 前記交差周波数結合の少なくとも一つの値が、前記状態空間モデルに含まれる第１波成分と第２波成分の間の位相振幅変調の推測値からなることを特徴とする請求項１に記載の方法。
3. 前記第１波が徐波成分であり、前記第２波がアルファ波成分であることを特徴とする請求項８に記載の方法。
4. 前記位相振幅変調の推定値が、前記徐波成分の位相と前記アルファ波成分の振幅に基づいて計算されることを特徴とする請求項３に記載の方法。
5. 前記交差周波数結合の少なくとも一つの値が、前記状態空間モデルに含まれる振動成分と前記ＥＥＧデータの周波数の帯域の間の相関値からなることを特徴とする請求項１に記載の方法。
6. 前記状態空間モデルの当て嵌めが、前記少なくとも一つの振動成分に含まれる各振動成分の調波振動数を当て嵌めることを含むことを特徴とする請求項１に記載の方法。
7. 前記少なくとも一つの振動成分がアルファ成分を含み、前記アルファ成分が中央周波数の事前分布に関連することを特徴とする請求項１に記載の方法。
8. 前記事前分布がフォンミーゼス事前分布であることを特徴とする請求項７に記載の方法。
9. 前記状態空間モデルの減衰係数が事前分布により制約されることを特徴とする請求項１に記載の方法。
10. 前記事前分布がベータ分布であることを特徴とする請求項９に記載の方法。
11. 患者の脳状態における少なくとも一つの振動成分を推定するためのシステムであって、
    脳波記録（ＥＥＧ）データを受け取るように構成された複数のセンサと、
    前記ＥＥＧデータを受け取るよう構成されたプロセッサと、
    表示部と、を備え、
    前記プロセッサは、
    少なくとも一つの振動成分を含む状態空間モデルを前記ＥＥＧデータに当て嵌めるステップと、
    前記状態空間モデルに基づいて前記ＥＥＧデータの交差周波数結合の少なくとも一つの値を推定するステップと、
    前記交差周波数結合の少なくとも一つの値に基づいて報告を生成するステップと、を実行するように構成され、
    前記表示部は前記報告を表示するよう構成された、システム。
12. 前記交差周波数結合の少なくとも一つの値が、前記状態空間モデルに含まれる第１波成分と第２波成分の間の位相振幅変調の推測値からなることを特徴とする請求項１１に記載のシステム。
13. 前記第１波が徐波成分であり、前記第２波がアルファ波成分であることを特徴とする請求項１８に記載のシステム。
14. 前記位相振幅変調の推定値が、前記徐波成分の位相と前記アルファ波成分の振幅に基づいて計算されることを特徴とする請求項１３に記載のシステム。
15. 前記交差周波数結合の少なくとも一つの値が、前記状態空間モデルに含まれる振動成分と前記ＥＥＧデータの周波数の帯域の間の相関値からなることを特徴とする請求項１１に記載のシステム。
16. 前記状態空間モデルの当て嵌めが、前記少なくとも一つの振動成分に含まれる各振動成分の調波振動数を当て嵌めることを含むことを特徴とする請求項１１に記載のシステム。
17. 前記少なくとも一つの振動成分がアルファ成分を含み、前記アルファ成分が中央周波数の事前分布に関連することを特徴とする請求項１１に記載のシステム。
18. 前記事前分布がフォンミーゼス事前分布であることを特徴とする請求項１７に記載のシステム。
19. 前記状態空間モデルの減衰係数が事前分布により制約されることを特徴とする請求項１１に記載のシステム。
20. 前記事前分布がベータ分布であることを特徴とする請求項１９に記載のシステム。

Images