TWI745321B - 決定網路連結之方法及系統 - Google Patents

Abstract

本發明係關於用於決定網路連結之方法及系統，特指但不限於決定稀疏網路中之網路連結之方法及系統，並特別應用於對於腦波資料（EEG data），於一方面，本發明提供一種識別方法，於同時產生訊號之互動節點之一網路中，識別所述節點間之連結以及估計經識別為相互連結之節點間之連結係數之方法，包含以下步驟，設置為零連結係數，其所計算之同調性（coherence）或部分同調性（partial coherence）係低於一第一預定閾值（threshold），隨後設置為零連結係數，其所估計之連結係數低於一第二預定閾值；以及重新估計連結係數尚未被設置為零之節點組合之連結係數。

Description

決定網路連結之方法及系統

本發明係關於決定網絡連結之方法及系統，特指但不限於決定稀疏網路中之網絡連結之方法及系統，並係可特別應用於腦波資料（EEG data）。

互動節點之網路中，各節點具有自身之動態（dynamics），此種網路係一用於描述複雜系統之重要計算工具（Strogatz, 2001）。根據特定用途，個別節點之動態、其耦合結構或其集群行為皆決定該系統之動態，在神經科學中，例如用於偵測訊號間之交互作用，特別是節點間之耦合結構等。了解大腦網路有助於揭露基於自然行為或特定疾病之生物學基礎（e.g. Hesse et al., 2003; Tass et al., 1998; Pitzalis et al., 1998; Keyl et al., 2000; Nollo et al., 2005; Bowers and Murray, 2004）。迄今已有許多技術被提出，藉由觀察訊號推斷複雜系統之網路結構。這些包含但不限於轉置熵（transfer entropy）（Schreiber, 2000; Staniek and Lehnertz, 2008）、狀態空間之循環（Arnhold et al., 1999; Chicharro and Andrzejak, 2009; Romano et al., 2007）、共同資訊（mutual information）（Pompe et al., 1998; Paluš and Stefanovska, 2003; Paluš and Vejmelka, 2007; Vejmelka and Paluš, 2008; Frenzel and Pompe, 2007）、相位動態（phase dynamics）（Rosenblum and Pikovsky, 2001; Rosenblum et al., 2002)、同調性（coherence） (Halliday and Rosenberg, 2000; Dahlhaus, 2000; Nolte et al., 2008）、福克-普朗克方程（the Fokker Planck formalism）（Prusseit and Lehnertz, 2008; Bahraminasab et al., 2009)、壓縮感知（compressed sensing） (Lee et al., 2011）或自我迴歸模型（autoregressive modelling）（Dahlhaus and Eichler, 2003; Eichler, 2000; Korzeniewska et al., 1997; Kamiński et al., 1997; Kamiński and Blinowska, 1991; Arnold et al., 1998）。

近年來可發現資料之可用性大幅增加，許多同時經記錄之通道亦增加。於網路分析方面，其導致估算高維度網路之挑戰。該目標係為僅估計網路中之直接連結。此外，其亦欲得出關於連結方向之論點。研究一影響方向之方法係使用因果關係之概念。許多方法（e.g. Hesse et al., 2003; Geweke, 1982, 1984; Chen and Wasterlain, 2006; Dhamala et al., 2008; Baccalá and Sameshima, 2001; Sameshima and Baccalá, 1999; Eichler, 2006; Kamiński and Blinowska, 1991）係立基於格蘭傑因果關係之定義（Granger, 1969）。簡言之，該定義表示若一過程x1係可用於預測x2之未來，則過程x1係可構成其他過程x2之原因。線性格蘭傑因果關係通常係藉由向量自我迴歸程序所建立，其透過多變量尤勒-沃爾克方程式（multivariate Yule-Walker equations）或相似方法所估計而得（Lütkepohl, 2005）。於大多數大型網路中之相鄰矩陣（adjacency matrix）係為稀疏，此意味於所有可能之連結中僅表現出其中少數者。

本發明人係已實現如何假設一稀疏網路可經使用於改善向量自迴歸程序之參數估計。

現有可被使用於決定一耦合結構之常用方法包含同調性（coherence）或部分同調性（partial coherence），其可如Schad等人所示（2009）受到估算。一進一步之技術係為定向淨相關（Eichler, 2005, 2006）。此等方法皆有其限制，如下更詳細敘述。特別地，現有方法一般係可在該方法開始明顯失去準確性和/或其長度成為不可接受或不符合計算要求前，僅處理潛在網路中10個最大節點。

於該整份文件中，儘管仍有其他方法可使用，但於此係根據魯克波爾（Lütkepohl，2005）估算自我迴歸係數。

於一第一情況下，分析一耦合白噪訊程序之15維度網路。圖1顯示該模擬網路之圖形。使用模擬資料，對於同調性與部分同調性之重建圖形分別顯示於圖2（a）與圖2（b）。除了與原始圖形非常不同外，兩重建圖形皆顯示相同數量之子圖。

於一第二情況下，係考量一階6維度自我迴歸程序，

(1)

(2) 其ε係一多變量高斯白噪訊程序（multivariate Gaussian white noise process）而模擬每一N=200數據點之M=100實現值（realisations）。該模擬網路係概括於圖3之圖形中，將真實自迴歸程序階層（p=1）用於估算。

對於所有36個係數而言，真實與平均估算係數間差異之絕對值係顯示於圖4中。誤差條（Error bars）係指100個實現值之平均值之標準差。該結果顯示即便在模擬中僅使用200個數據點，該所有估計之係數係非常接近其自身之真值（true values）。於此所研究之系統係一稀疏系統，即便其相對低維度，一般參數估計仍可進行處理。

於一第三種情況下係使用一較高維度系統。該係為一耦合白噪訊程序之40維度網路。連結係存在於落差（lag）1或落差2其中之一。圖5中所顯示係該模擬網路之圖形。模擬該系統之N=10,000數據點。基於同調性與部分同調性所估計之網路係分別顯示於圖6（a）與圖（b）中。其皆未獲得一該基礎網路之意義表徵（meaningful representation）。

對於定向淨相關（directed partial correlation）分析（Eichler, 2005, 2006），使用p=2階之真程序。所估計出之網路係顯示於圖7中。

一些附加連結顯示於該分析中（於圖7中以虛線箭頭表示）。該些係為偽陽性（false positive）結論。其發生係因大量需於高維度網路中估計之係數。

腦波圖（Electroencephalography，EEG）提供大腦活動之多通道資料，該大腦活動係由一多個（通常至少為20個）貼附至一個人頭皮之小型感測器所偵測，該感測器係偵測當腦細胞傳遞訊息於彼此時，腦神經元中離子流動所產生之電壓波動（voltage fluctuations）。由於多個感測器，來自腦波資料具有一相對應數個通道。腦波圖目前常用於幫助診斷及監測影響大腦之眾多病症，尤其係指癲癇（epilepsy）。

一段時間以來，腦波圖已於臨床上應用作為一大腦功能之測量，希望能確定及區分大腦特定功能狀態，然而至今仍進展緩慢。

本發明人已實現精確地決定腦波資料中因果關係之網路其可能性，可進一步解釋該資料以用於臨床目的。

腦波資料係多維度且易於以網路方式分析。現今許多腦波資料類型具有大量通道（如20個或更多）。傳統分析方法難以提供有意義的資訊或對於該多變量數據之解釋，因此需要新方法來處理所觀察之多階系統。

特別於阿茲海默症之診斷和治療方面，現已認為在明顯症狀被觀察出前，阿茲海默症係長時間發展出（可能約20年）。 因此，有相當多的關注於能提供對於潛在患者或其表現特別易於罹患或風險因子之可靠的早期識別之技術。該些技術旨在識別完全前臨床階段（通常係發病前10-20年）或前驅階段之早期阿茲海默症發病。

若存在用於穩健型資料分析之合適方法，使用腦波資料用來識別阿茲海默症早期症狀或警訊具有相當吸引力，因其係一標準且廣泛使用和取得之技術。

本發明之目的在於提供關於網路連結和係數準確且可信之預測之方法及系統，特別於稀疏網路中。

本發明進一步目的在於提供預測網路連結和係數之有效方法。

本發明進一步目的在於提供用於處理腦波資料之方法及系統，藉由大腦活動之網路總覽提供有意義資訊和/或對於資料上解釋，以及允許該資料後續之應用。

廣泛而言，本發明提供用於識別一網路中節點間連結之方法及系統，其係識別可能為零連結之係數並將其設置為零以便於後續處理。

第一方面，本發明提供一種於同時產生訊號之互動節點之一網路中，識別所述節點間之連結以及估計經識別為相互連結之節點間之連結係數之方法，該方法包含步驟：於一段預定時間內定期記錄每一節點之訊號以形成一資料集；計算該資料集中每一節點組合間之同調性（coherence）與部分同調性（partial coherence）；檢驗每一節點組合，假如該同調性或部分同調性其中之一係低於一第一預定閾值（threshold），則對於所有後續步驟其對應之連結係數設置為零；一第一估計步驟，估計該資料集中連結係數尚未被設置為零之該節點組合之連結係數；對於經該所述第一估計步驟所估計低於一第二閾值之每一連結係數，於所有後續步驟將所述之係數設置為零；以及一第二估計步驟，重新估計該資料集中連結係數尚未被設置為零之該節點組合之連結係數。

該方法於估計連結係數之最後步驟前包含兩個「歸零（zeroing）」步驟。此等步驟係以移除間接連接（考量部分同調性與同調性），可藉此減少或消除來自該經測定網路之偽陽性（false positives）。

此方面之方法僅使用資料，其不依賴任何關於該基本模型之預測或假設。

該方法方面於假設一稀疏程度於該網路中，此意味該相鄰矩陣（adjacency matrix）中某些係數為零。基於該假設，藉由識別係數以增進估計程序，該係數係作為歸零與進行進一步計算或估計前將該些設置為零之候選。因此，該方法方面比起現有方法更有效於預測一稀疏網路之連結。因而該方法較佳應用於已知或預測為稀疏連結之網路上。藉由稀疏連結，我們意指該網路未表現（如該連結係數為零）之成對節點間具有至少50%之潛在連結，較佳至少為60%，以及於某些實施例中至少為75%。實際上，隨著網路變得稀疏，該方法變得更加有效，因此可應用於未表現之80%或90%之潛在連結之網路上。

該方法能使格蘭傑因果關係（Granger-causality）應用於高維度系統（特別是該些具有10個或更多節點之系統）上。雖然現有之格蘭傑因果關係推論通常於低維度系統上運作良好，但本方法於減少相關係數數量之附加步驟使格蘭傑因果關係能運用於更高維度系統上，特指（但不限於）該些稀疏連結之系統。

於詳細說明中所列之模擬於一模擬研究中證明，根據實施例之方法勝過標準方法且避免關於格蘭傑因果關係之偽陽性結論。與格蘭傑因果關係推論之簡易應用相比較，該些模擬展現該方法之優越性。

因此，該方法可實現格蘭傑因果關係之一可靠估計。該方法可容易地應用各種不同測量方法於格蘭傑因果關係上以及其他基於向量自我迴歸模型之方法。

較佳地，該檢驗步驟包含計算每一節點組合所計算之同調性與部分同調性之乘積，以及決定所述乘積是否低於所述第一預定閾值。假如該同調性或部分同調性其中之一係為零，或接近零，而該乘積結果將為零，或接近零。這意指對於每一係數而言，僅需要與該閾值進行單一比較。

如Schad等人（2009）所定義，該第一預定閾值可以是部分同調性之臨界值（critical value）。或者，該第一預定閾值可以是同調性之臨界值。又或者，該第一預定閾值可以是該部分同調性之臨界值與該同調性之臨界值兩者之乘積。

該估計連結係數之第二閾值可由以下步驟決定： 根據該所估計之係數之歐幾里德距離（Euclidean distance）之平方，將該所估計之係數分成兩群組，一包含該些具有高數值之係數之第一群組，與一包含該些具有低數值之係數之第二群組；以及設置該所述第二閾值作為一大於該所述第二群組中所有係數值之一數值。

以此種方式，該所估計之連結係數可被分成兩群組，以及弱連結與強連結間之分別，因此可於其可能是完全歸因於噪訊所估計之連結與該些代表一真正連結間取得一明確區分。因此，該第二閾值於適當程度上可變地經過選擇，以分離該兩群組。或者，可預先設置該第二閾值。

較佳地，估計該連結係數之該第一與第二估計步驟係估算該資料集之自我迴歸係數。

較佳地，該方法進一步包含用以移除離群值之篩選數據集步驟。隨著該方法對於所測量之資料產生作用，其可能易受資料中離群值所影響。離群值係資料中之假影（artefacts），其通常係因非為欲測量以作為記錄數據之一部分所造成。舉例而言，於腦波資料中，眨眼（eye-blinks）可導致該種類之假影產生，於資料集中移除該種離群值係可因此提高該方法之準確性。

較佳地，該方法進一步包含用以移除噪訊之過濾資料集步驟。同樣地，該方法對於所測量之資料產生作用，於該些測量中其可能對於噪訊敏感。因此，過濾數據集以移除噪訊可提高該方法之準確性。

該篩選或過濾動作可於計算步驟前進行，或可併入該係數之實際估算中。

於特定實施例中，產生訊號之互動節點網路係為一種腦波系統。於腦波資料上之網路結構分析可提供對於大腦活動與肌肉活動之深入了解，且所得之網路係可使用於比較目的，例如，針對特定人群之樣本網路，或作為相同個體之未來研究比較。

腦波圖之時間解析度（temporal resolution）係於毫秒（millisecond）範圍內。已知大腦處理時間約為500毫秒， 因此本發明之方法係更可應用於此種資料。然而，該等技術同樣適用於具有較低解析度之其他資料（如功能性磁振造影（functional magnetic resonance imaging）或具有一約為2秒時間解析度之fRMI）。

目前腦波資料一般係記錄超過20分鐘之週期。這可能於該段時間內觀察處於一恆定狀態（或多狀態）患者時導致實際資料收集之問題，以及增加假影出現之可能性。假如該時間段可進一步減少至或許少數幾個一百秒之單位內，則係可減少與/或避免此等問題。

本方法可由相對少量資料提供一網路之穩健預測（robust prediction），可潛在減少該腦波資料所需數量（因而該時間長度）。

本方法可包含或不包含上述較佳與選擇性特徵之全部或其等之組合。

第二方面，本發明提供一監測一病患大腦功能之方法，該方法包含以下步驟：於一段時間內對於患者執行一腦波記錄；使用根據前述第一方面之一方法識別腦波圖上節點訊號間之網路連結，可包含或不包含該些選擇性或較佳特徵之全部或部分。

上述之方法較佳者係透過根據本發明第三方面之一系統實施，如下所述，但非必須。

本發明進一步包含，運行於電腦系統上之電腦程式，該電腦程式能執行上述之方法，其可包含或不包含該些較佳與選擇性特徵之全部或部分。

第三方面，本發明提供一系統，該系統係用於識別網路連結與估計大腦活動資料記錄中節點間之連結係數，該系統包含：複數感測器，於一預定時間內記錄一個體於不同位置之腦部活動以產生一資料集；以及一處理器，其經設置以：計算該資料集中每一節點組合間之同調性與部分同調性；檢驗每一節點組合，假如該同調性或部分同調性其中之一低於一第一預定閾值，則對於所有後續步驟其對應之連結係數設置為零；估計該資料集中連結係數尚未被設置為零之該節點組合之連結係數；每一所估計低於一第二閾值之連結係數，對於後續步驟設置所述係數為零；以及重新估計該資料集中連結係數尚未被設置為零之該節點組合之連結係數。

該系統處理所記錄之資料並於估計該連結係數之最後步驟前實施兩個「歸零」步驟，其等係移除間接連接（藉由將部分同調性與同調性納入考量）。此可減少或消除來自該所決定之網路之偽陽性。

該系統係僅使用資料，其並不依賴任何關於基本模型之預測或假設。

該系統假設於該網路中有一稀疏程度，此意味該相鄰矩陣中某些係數為零。基於該假設，藉由識別係數以增進該估計程序，該係數係作為歸零與進行進一步計算或估計前將該些設置為零之候選。因此，該方面之方法比起現有方法更有效於預測一稀疏網路之連結。因而該方法較佳應用於已知或預測為稀疏連結之網路上。藉由稀疏連結，我們意指該網路未表現（如該連結係數為零）之成對節點間具有至少50%之潛在連結，較佳係至少為60%，以及於某些實施例中至少係為75%。實際上，隨著網路變得稀疏，該方面之方法系可更加有效，因此可應用於未表現之80%或90%之潛在連結之網路上。

該系統之處理器能應用格蘭傑因果關係於高維度系統（特別是該些具有10個或更多節點之系統）上。雖然現有之格蘭傑因果關係之資料可於低維度系統上運作良好，但減少相關係數數量之程序係使格蘭傑因果關係可運用於較高維度系統上，特指（但不限於）該些稀疏連結之系統。

因而，該系統能可靠地估計格蘭傑因果關係與可容易地應用於多種格蘭傑因果關係之測量方法，以及其他立基於向量自我迴歸模型之方法。

大腦活動資料上之網路結構分析可提供對於大腦活動與肌肉活動一深入了解，以及所得之網路可作為比較目的上使用，如針對特定人群之樣本網路，或作為對於相同個體未來研究之比較器（comparators）。

較佳地，該處理器係設置用以計算每一節點組合所計算之同調性與部分同調性之乘積，以及決定所述乘積是否低於所述第一預定閾值。假如該同調性或部分同調性其中之一係為零，或接近零，而該乘積結果將為零，或接近零。這意指對於每一係數而言，僅需要與該閾值進行單一比較。

如Schad等人（2009）所定義，該第一預定閾值可以是部分同調性之臨界值。或者，該第一預定閾值可以是同調性之臨界值。又或者，該第一預定閾值可以是該部分同調性之臨界值與該同調性之臨界值兩者之乘積。

該處理器可被組設用以決定所述第二閾值，藉由：根據該所估計之係數之歐幾里德距離之平方，將該估計係數分成兩群組，一包含該些具有高數值之係數之第一群組，與一包含該些具有低數值之係數之第二群組；以及設置該所述第二閾值作為一大於該所述第二群組中所有係數值之一數值。

以此種方式，該估計連結係數可被分成兩群組，以及弱連結與較強連結間之分別，因此可於該其可能完全歸因於噪訊之估計連結與該些代表一真正連結間取得一明確分別。因此，該第二閾值於適當程度上可變地經過選擇以分離該兩群組。或者，可預先設置該第二閾值。

較佳地，該處理器係經設置估計該連結係數，藉由估計該資料集之自我迴歸係數。

較佳地，該處理器係經設置用以篩選該資料集以移除離群值。該處理器係處理所測量之腦部活動資料，其可能易受該資料中離群值影響。離群值係該資料中之假影，其通常係因非為欲測量以作為記錄數據之一部分所造成。舉例而言，於腦波圖資料中，眨眼特別可能導致該種類之假影產生。於資料集中移除該離群值可因此提高該系統之準確性。

較佳地，該處理器係經設置用以過濾該資料集以移除噪訊。同樣地，作為該處理器處理所測量之腦部活動資料，其可能易受該些測量方式中噪訊影響。因而，過濾該資料集以移除噪訊可提高該系統之準確性。

該篩選或過濾動作可於計算步驟前進行，或可併入該係數之實際估計中。

較佳地，該系統係應用於腦波資料，而該等複數感測器係一腦波儀。

腦波圖之時間解析度係於毫秒範圍內。已知大腦處理時間約為500毫秒，因此本發明之方法係更佳能應用於該資料上。然而，該些技術同樣適用於具有較低解析度之其他資料（如功能性磁振造影或具有一約為2秒時間解析度之fRMI）。

腦波資料一般係以超過20分鐘之週期方式記錄。這可能導致於該段時間內觀察處於一恆定狀態（或多狀態）之患者時導致實際資料收集之問題，以及增加假影出現之可能性。假如該時間段可進一步減少至少數幾個一百秒單位內，則係可減少與/或避免此等問題。

本發明之系統可由相對少量資料提供一該網路之穩健型預測，該腦波資料所需數量（因而該時間長度）可潛在地減少。

本發明之系統可包含或不包含上述該些較佳與選擇性之特徵。

本發明之系統可藉由執行根據上述本發明之第一或第二方面之方法來操作，但非必須。

於以下討論中，該結果係顯示對於格蘭傑因果關係之一特定測量，即所謂的定向淨相關（DPC）（Eichler，2005，2006）。然而，該結果應用於任何基於向量自我迴歸程序之格蘭傑因果關係測量。

根據本發明一實施例之一方法係圖表式地顯示於圖8之流程圖中。該方法涉及一三步驟方法以估計稀疏自我迴歸程序。該基本原則係於該實際擬合過程（actual fitting procedure）前排除一些係數。

該第一分析步驟（S102），係對於根據Schad等人研究之方法下之X與Y（Schad等人，2009）估計同調性（Priestley，1981）

(3) 與部分同調性（Halliday等人，1995）

(4)

使用一閾值，對於同調性與部分同調性乘積係為與零相容之所有係數（如，落入一預定上限內，例如固定為零（S103））。同調性與部分同調性之乘積將為零，假如兩者其中之一係為零。該些係數係維持於零作為該程序之提示。

使用作為決定一係數是否與零相容之該閾值，係可為該部分同調性之臨界值（critical value）或同調性之臨界值，如Schad等人所定義（2009）。

於該第二步驟（S104）自我迴歸係數係由下式估計（Lütkepohl，2005）

(5) 與

(6) 而

。於前一步驟識別為具有同調性與部分同調性非重要性乘積之係數保持為零。

該所得之係數根據其歐幾里德距離平方分成兩個群集。將該群集中具有較小值之係數設置為零（S105）。該步驟考量到同調性與部分同調性係對稱測量（symmetric measures）之事實，因而不可排除單一方向連結。

於第三步驟，再次估算自我迴歸係數（方程式（5）與（6）－S106）。此次，所有於該第一或第二步驟中被識別為與零相容之所有係數保持為零。

該方法之重要優點在於僅有非零係數於第三步驟受到估計。由於所估計之係數之數量於該程序中大量減少，得以提升該估計之準確性，意即從相同數量之資料點係有較少之係數被估計。

該方法之績效可由下列模擬數據發現。模擬（ Simulations ）

對於先前分析之耦合白噪訊程序之15維度與40維度網路（圖1和圖5），使用上述實施例之方法估計定向淨相關。該結果圖形分別顯示於圖9和圖10中。兩個例子皆正確揭露真實基礎網路結構。

應強調者，該等用於模擬之網路係為隨機生成。本發明人已於超過100種不同稀疏隨機網路中測試上述實施例之方法，並皆正確地重建（結果未顯示於此）。

該稀疏估計技術亦改進於較低維度稀疏系統中對於向量自我迴歸係數之估計。例如，如下所述之一階6維度自我迴歸程序（方程式（1）與（2）），圖11係與圖4相同，並且顯示對於36個係數而言，使用根據上述實施例方法於真實與平均所計算之係數間差異之絕對值。誤差條係指100個實現值之平均標準差。未顯示誤差條係該些固定於零之係數，於模擬中，該些係數皆為真正零。結果顯示該6個非零係數之估計係非常接近其各自之真值（true value）（並且比前述之程序更加接近），同時該30個為零之係數皆經正確識別恰好為零。

圖12係耦合白噪訊程序一更進一步模擬之15維度網路。

圖13顯示僅使用同調性預測圖12之網路其連結與係數。由此可知，除了於該網路中更為複雜部分加入數個附加連結，該方法亦錯誤地表徵節點1、6、8與10間之關係。

圖14顯示使用部分同調性預測圖12之網路之連結與係數。該方法係更加成功，因其移除數個間接連接。然而，其仍錯誤地預測該網路中較為複雜部分之連結，該連結並未出現於基礎網路中（偽陽性），以及如同調性方法，錯誤表徵結點1、6、8與10間之關係。

圖15顯示使用上述實施例所述之方法計算圖12之網路其連結與係數。與圖12相比較後可以看出該網路完整且精確地表示出。

圖16係顯示耦合白噪訊程序一更進一步模擬之40維度網路。

圖17顯示使用同調性所預測之網路，由此可知，隨著網路維度增加，該方法中偽陽性所產生之數量完全覆蓋該真實連結。圖18顯示使用部分同調性之預測網路。如同前述之模擬，該方法去除大量偽陽性連結，但仍無法完全再現該基礎網路。該上述實施例所述之方法準確地再現圖16所示之網路。用途範例（ Example Application ）

於一用途範例中，分析一健康自願者之腦波圖資料。腦波圖之記錄係於闔眼期間所獲得。該訊號係於512Hz進行採樣。根據10-20系統將20個電極（electrodes）放置於頭皮上。資料係往下採樣至200Hz以及於稀疏估計中使用一p=2階之模型。分析一100秒之片段。圖19顯示該交互結構。於一健康個體中，如預期的，於大腦許多部分間存在一寬範圍之互連（interconnections）。

圖20與圖21示意性地描述一健康個體與一具有輕度認知損傷個體所預期之網路類型。該等節點經簡化至額葉（frontal）區、中間（central）區、顳葉（temporary）區與枕葉（occipital）區（4通道腦波儀）。如圖20所示，一健康個體預期於該所有區域與每一個別區域間具有互連（因節點有效地相同，其並未標記於圖20中）。

圖21顯示，於一具有輕度認知損傷個體中，雖然每一區域內存在連結，但每一區域間之互連明顯減少，並且通常僅涉及從額葉至枕葉單一方向之少許流動，或無從枕葉區或顳葉區之回流，以及非相鄰區間無連結。

雖然4通道腦波儀並不產生前述之高維度資料，但分析該腦波監測之資料可能特別有用，因4通道腦波儀係易於取得，其容易使用於初級護理裝置上與透過個人自行使用。進一步發展（ Futher developments ）

本發明之稀疏估計技術係可容易地應用於格蘭傑因果關係之頻域測量（frequency-domain measures），如再標準化之部分定向同調性（partial directed coherence）（Schelter等人，2009）。這將能研究不同頻率下之網路，特別是關於腦波圖分析之頻率。

應強調的是上述實施例之方法係由資料所驅動。此意味著該稀疏係數矩陣係基於所測量節點所得之資料集所估計出。根據估計資料集之同調性與部分同調性來放置零。若基礎網路非稀疏，該演算法將會以基礎估計方式簡易地對於所有係數進行估計。

然而，測量方法係永不準確且觀察噪訊會影響格蘭傑因果關係推論（Sommerlade等人，2015）。為了處理觀察噪訊，於此之稀疏估計技術可與Sommerlade等人（2015）所提出之狀態空間模型（state space model ）方法相結合。

除了觀察噪訊外，受測數據會受到離群值影響。例如，對於腦波資料而言，該離群值包含眨眼之假影。當估計自我迴歸程序時，如眨眼該離群值可藉由使用加權穩健型卡曼濾波器（weighted robust Kalman filtering）（Ting等人，2007）或多個一般離群值穩健型卡曼濾波器（Agamennoni等人，2011）移除。於此所述之方法可與該些方法結合以進一步增進該估計技術。進一步用途（ Further uses ）

根據本發明實施例之預測方法可於多種裝置上使用。如上所述，其特別使用於腦波資料之分析上，根據所預測之網路，其可進行進一步評估或確定所得資料來源之個體其認知功能。

該網路預測方法可以多種方式使用。除了診斷之外（例如藉由將來自健康個體以及具認知損傷個體之腦波圖資料所取得之網路相互比較），該網路之預測可使用於監測一個體之反應以對於認知損傷進行治療。舉例來說，一來自一受治療個體其衍生自腦波資料之網路之一記錄可被維持並隨時間監測其變化。一有效之治療可減緩或停止該個體認知功能之惡化（藉由該網路中連結之數量與/或長度顯示），或使前述下降反轉（藉由該網路中連結之數量與/或長度顯示）。

該網路預測方法亦可使用作為一鑑別器（discriminator）用以識別或篩選用於特別治療試驗之候選者，藉由識別具有神經層面上認知功能之特定特徵之個體。該網路預測方法亦可使用於篩選以確保所有受試者具有相同或相似認知損傷，或選擇具有一廣範圍認知損傷之參與者。

該網路預測方法亦可結合使用於目標反應問題以測試特定反應以及決定於測試過程中是否影響大腦內其反應之模式。一般條件（ General provisions ）

除了所述結構元件與使用者相互作用外，上述實施例之系統及方法可於一電腦系統上全部或部分執行（特指於電腦硬體或電腦軟體）。

所謂「電腦系統」，包括硬體、軟體與儲存裝置係用於實施一系統或執行根據上述之一方法。例如，一電腦系統可包含一中央處理器（central processing unit，CPU），輸入方法（input means）、輸出方法（output means）與資料儲存（data storage）。較佳地，該電腦系統具有一顯示器以提供一視覺輸出顯示。該資料儲存可包括隨機存取記憶體（RAM），磁碟驅動器（disk drives）或其他電腦可讀取媒體（computer readable media）。該電腦系統可包含多個與網路相連接之電腦裝置並能藉由該網路彼此溝通。

上述實施例之方法可提供作為電腦程式（computer program）或電腦程式產品或執行一所安裝電腦程式之電腦可讀取媒體，當運行於一電腦上以執行上述之方法。

所謂「電腦可讀取媒體」，包含但不限於，任何可由一電腦或電腦系統直接讀取與進入之非暫態介質（non-transitory medium）或媒體。該媒體可包含但不限於，磁性儲存媒體，如軟磁碟（floppy discs）、硬碟儲存媒體（hard disc storage media）與磁帶（magnetic tape）；光學儲存美體（optical storage media），如光碟（optical discs）或唯讀式光碟（CD-ROMs）；電子儲存媒體（electrical storage media），如記憶體（memory），包含RAM、ROM與快閃記憶體（flash memory）；以及上述如磁性/光學儲存媒體之混合與組合。

雖然透過上述之例示性實施例之結合描述本發明，許多同等之修改與變化對於所屬技術領域之人係顯而易見。因此，上述所列本發明之例示性實施例係為說明性質而非用以限制。在不脫離本發明之精神與範圍下可對於所述實施例進行各種改變。參考文獻（ REFERENCES ）

離群值-穩健卡曼濾波器（Agamennoni, G., Nieto, J., Nebot, E., May 2011. An outlier-robust kalman filter. In: Robotics and Automation (ICRA), 2011 IEEE International Conference on. pp. 1551-1558.） 一種偵測相互依賴之穩健方法：應用於顱內記錄之腦波圖（Arnhold, J., Grassberger, P., Lehnertz, K., Elger, C. E., 1999. A robust method for detecting interdependencies: Application to intracranially recorded EEGs. Physica D 134, 419-430.） 藉由卡曼濾波器之非平穩時間序列之適應性AR模型（Arnold, M., Miltner, W. H. R., Witte, H., Bauer, R., Braun, C., 1998. Adaptive AR modeling of nonstationary time series by means of Kalman filtering. IEEE T. Bio-Med. Eng. 45, 553-562.） 部分定向同調性：確定神經結構之新概念（Baccalá, L. A., Sameshima, K., 2001. Partial directed coherence: A new concept in neural structure determination. Biol. Cybern. 84, 463-474.） 大腦動態物理學：福克-普朗克分析揭示麻醉期間腦波δ與θ之變化（Bahraminasab, A., Ghasemi, F., Stefanovska, A., McClintock, P. V. E., Friedrich, R., 2009. Physics of brain dynamics: Fokker-Planck analysis reveals changes in EEG delta and theta activity during anaesthesia. New J. Phys. 11, 103051.） 心臟跳動間隔變化與收縮壓變化間之交互作用（Bowers, E. J., Murray, A., 2004. Interaction between cardiac beat-to-beat interval changes and systolic blood pressure changes. Clin. Autonom. Res. 14, 92-98.） 癲癇狀態：成人病理生理與管理（Chen, W. Y., Wasterlain, C. G., 2006. Status epilepticus: Pathophysiology and management in adults. Lancet Neurol. 5, 246-256.） 使用秩統計量定向耦合之可信賴偵測（Chicharro, D., Andrzejak, R. G., 2009. Reliable detection of directional couplings using rank statistics. Phys. Rev. E 80, 026217.） 多變量時間序列之圖形交互模型（Dahlhaus, R., 2000. Graphical interaction models for multivariate time series. Metrika 51, 157-172.） 時間序列之因果關係及圖形模型（Dahlhaus, R., Eichler, M., 2003. Causality and graphical models for time series. In: Green, P., Hjort, N., Richardson, S. (Eds.), Highly Structured Stochastic Systems. Oxford University Press, pp. 115-137.） 從傅立葉與小波轉換之時間序列資料估計格蘭傑因果關係（Dhamala, M., Rangarajan, G., Ding, M., 2008. Estimating Granger causality from Fourier and wavelet transforms of time series data. Phys. Rev. Lett. 100, 018701.） 圖形時間序列模型之馬可夫特性（Eichler, M., 2000. Markov properties for graphical time series models. Preprint University of Heidelberg.） 用於評估神經系統內有效性連結之圖形方法（Eichler, M., May 2005. A graphical approach for evaluating effective connectivity in neural systems. Philos Trans R Soc Lond B Biol Sci 360 (1457), 953-967.） 多變量時間序列中圖形模型之動態關係（Eichler, M., 2006. Graphical modeling of dynamic relationships in multivariate time series. In: Schelter, B., Winterhalder, M., Timmer, J. (Eds.), Handbook of Time Series Analysis. Wiley-VCH, Ch. 14, pp. 335-372.） 用於多變量時間序列耦合分析之部分共同資訊量（Frenzel, S., Pompe, B., 2007. Partial mutual information for coupling analysis of multivariate time series. Phys. Rev. Lett. 99, 204101.） 多重時間序列間線性相關與反饋之測量（Geweke, J., 1982. Measurement of linear dependence and feedback between multiple time series. J. Am. Stat. Assoc. 77, 304-313.） 時間序列間條件線性相關與反饋之測量（Geweke, J., 1984. Measures of conditional linear dependence and feedback between time series. J. Am. Stat. Assoc. 79, 907-915.） 藉由計量經濟模型與交譜方法研究因果關係（Granger, J., 1969. Investigating causal relations by econometric models and cross-spectral methods. Econometrica 37, 424-438.） 關於同調性與組合同調性之應用與估計與解釋（Halliday, D. M., Rosenberg, J. R., 2000. On the application and estimation and interpretation of coherence and pooled coherence. J. Neurosci. Meth. 100, 173-174.） 用於分析混和時間序列/點過程資訊之架構-生理震顫、單機電單元放電與肌電圖之理論和應用之研究（Halliday, D. M., Rosenberg, J. R., Amjad, A. M., Breeze, P., Conway, B. A., Farmer, S. F., 1995. A framework for the analysis of mixed time series/point process data - Theory and application to the study of physiological tremor, single motor unit discharges and electromyograms. Prog. Biophys. molec. Biol. 64, 237-278.） 使用時變腦波格蘭傑因果關係檢驗神經群之定向相互依賴性（Hesse, W., Möller, E., Arnold, M., Schack, B., 2003. The use of time-variant EEG Granger causality for inspecting directed interdependencies of neural assemblies. J. Neurosci. Meth. 124, 27-44.） 一種描述大腦結構內資訊流動之新方法（Kamiński, M. J., Blinowska, K. J., 1991. A new method of the description of the information flow in the brain structures. Biol. Cybern. 65, 203-210.） 睡眠與覺醒期間腦波之同調性與傳遞之拓樸圖形分析（Kamiński, M. J., Blinowska, K. J., Szelenberger, W., 1997. Topographic analysis of coherence and propagation of EEG activity during sleep and wakefulness. Electroenceph. Clin. Neurophys. 102, 216-227.） 正弦壓力受器刺激與固定頻率呼吸期間心臟循環耦合（Keyl, C., Dambacher, M., Schneider, A., Passino, C., Wegenhorst, U., Bernardi, L., 2000. Cardiocirculatory coupling during sinusoidal baroreceptor stimulation and fixed-frequency breathing. Clinical Science 99, 113-124.） 於大鼠各類型行為期間海馬迴與相關結構間之資訊流動（Korzeniewska, A., Kasicki, S., Kamiński, M. J., Blinowska, K. J., 1997. Information flow between hippocampus and related structures during various types of rat's behavior. J. Neurosci. Meth. 73, 49-60.） 壓縮感之下稀疏腦部網路之恢復（Lee, H., Lee, D. S., Kand, H., Kim, B.-N., Chung, K., 2011. Sparse brain network recovery under compressed sensing. IEEE T. Med. Imaging 30, 1154-1165.） 多時間序列分析之新概述（Lütkepohl, H., 2005. New Introduction to Multiple Time Series Analysis. Springer, pp. 82-87.） 探索自發性心臟週期之方向性與人體收縮壓波動變異性交互作用：評估感壓反射增益影響（Nollo, G., Faes, L., Porta, A., Antolini, R., Ravelli, F., 2005. Exploring directionality in spontaneous heart period and systolic pressure variability interactions in humans: Implications in the evaluation of baroreflex gain. Am. J. Physiol. Heart. Circ. Physiol. 288, 1777-1785.） 於複雜物理系統中穩健估計資訊流動方向（Nolte, G., Ziehe, A., Nikulin, V. V., Schlögl, A., Krämer, N., Brismar, T., Müller, K.-R., 2008. Robustly estimating the flow direction of information in complex physical systems. Phys. Rev. Lett. 100, 234101.） 來自相互作用震盪器相位之耦合方向：一種資訊理論方法（Paluš, M., Stefanovska, A., 2003. Direction of coupling from phases of interacting oscillators: An information-theoretic approach. Phys. Rev. E 67, 055201(R).） 來自雙變量時間序列之耦合方向性：如何避免錯誤因果關係與遺漏之連結（Paluš, M., Vejmelka, M., 2007. Directionality of coupling from bivariate time series: How to avoid false causalities and missed connections. Phys. Rev. E 75, 056211.） 呼吸頻率對於RR間隔與收縮壓波動間關係之影響：一種頻率依賴性現象（Pitzalis, M. V., Mastropasqua, F., Massari, F., Passantino, A., Colombo, R., Mannarini, A., Forleo, C., Rizzon, P., 1998. Effect of respiratory rate on the relationships between RR interval and systolic blood pressure fluctuations: A frequency-dependent phenomenon. Cardiovasc. Res. 38, 332-339.） 使用共同資訊量以測量心臟循環系統中之耦合（Pompe, B., Blidh, P., Hoyer, D., Eiselt, M., 1998. Using mutual information to measure coupling in the cardiorespiratory system. IEEE Eng. Med. Biol. Mag. 17, 32-39.） 譜相分析與時間序列（Priestley, M. B., 1981. Spectral Analysis and Time Series. Academic Press, London, pp 655-668.） 使用所估計之福克－普朗克係數測量耗散動態系統中之相互依賴性（Prusseit, J., Lehnertz, K., 2008. Measuring interdependences in dissipative dynamical systems with estimated Fokker-Planck coefficients. Phys. Rev. E 77, 041914.） 藉由條件概率之循環估計耦合之方向（Romano, M. C., Thiel, M., Kurths, J., Grebogi, C., 2007. Estimation of the direction of the coupling by conditional probabilities of recurrence. Phys. Rev. E 76, 036211.） 耦合方向識別：應用於心臟循環交互作用（Rosenblum, M. G., Cimponeriu, L., Bezerianos, A., Patzak, A., Mrowka, R., 2002. Identification of coupling direction: Application to cardiorespiratory interaction. Phys. Rev. E 65, 041909.） 偵測交互作用震盪器內耦合方向（Rosenblum, M. G., Pikovsky, A. S., 2001. Detecting direction of coupling in interacting oscillators. Phys. Rev. E 64, 045202(R).） 使用部分定向同調性描述神經群交互作用（Sameshima, K., Baccalá, L. A., 1999. Using partial directed coherence to describe neuronal ensemble interactions. J. Neurosci. Meth. 94, 93-103.） 神經網路中偵測直接定向交互作用之方法（Schad, A., Nawrath, J., Jachan, M., Henschel, K., Spindeler, L., Timmer, J., Schelter, B., 2009. Approaches to the detection of direct directed interactions in neuronal networks. In: Velazquez, J. L. P., Wennberg, R. (Eds.), Coordinated Activity in the Brain. Springer, Ch. 3, pp. 43-64.） 使用重整化部分定向同調性評估神經訊號間之定向影響之強度（Schelter, B., Timmer, J., Eichler, M., 2009. Assessing the strength of directed influences among neural signals using renormalized partial directed coherence. J. Neurosci. Meth. 179, 121-130.） 測量資訊轉移（Schreiber, T., 2000. Measuring information transfer. Phys. Rev. Lett. 85, 461-464.） 評估神經訊號間定向影響之強度：一種噪訊資料之方法（Sommerlade, L., Thiel, M., Mader, M., Mader, W., Timmer, J., Platt, B., Schelter, B., Jan 2015. Assessing the strength of directed influences among neural signals: an approach to noisy data. J Neurosci Methods 239, 47-64.） 符號轉置熵（Staniek, M., Lehnertz, K., 2008. Symbolic transfer entropy. Phys. Rev. Lett. 100, 158101.） 探索複雜網路（Strogatz, S. H., Mar 2001. Exploring complex networks. Nature 410 (6825), 268-276.） 從噪訊資料中偵測n:m相位鎖定：應用於腦磁波儀（Tass, P., Rosenblum, M. G., Weule, J., Kurths, J., Pikovsky, A. S., Volkmann, J., Schnitzler, A., Freund, H. J., 1998. Detection of n : m phase locking from noisy data: Application to magnetoencephalography. Phys. Rev. Lett. 81, 3291-3295.） 用於穩健離群值偵測之卡曼濾波器（Ting, J.-A., Theodorou, E., Schaal, S., Oct 2007. A Kalman filter for robust outlier detection. In: Intelligent Robots and Systems, 2007. IROS 2007. IEEE/RSJ International Conference on. pp. 1514-1519.） 藉由條件共同資訊量推論耦合之方向性（Vejmelka, M., Paluš, M., 2008. Inferring the directionality of coupling with conditional mutual information. Phys. Rev. E 77, 026214.）

無

本發明之實施例將參照以下附圖方式進行說明： 圖1係為耦合白噪訊過程之一模擬15維度網路之圖形； 圖2（a）與圖2（b）係分別為使用同調性與部分同調性重建圖1網路之圖形； 圖3係一模擬6維度一階自我迴歸程序之圖形； 圖4係顯示圖3程序其真實係數與使用p=1階真實自我迴歸程序估算之平均係數兩者間差異之絕對值； 圖5係耦合白噪訊過程之一模擬40維度網路之圖形； 圖6（a）與圖6（b）係分別為使用同調性與部分同調性重建圖5網路之圖形； 圖7係為使用p=2階真實程序之定向淨相關分析圖5網路重建之圖形； 圖8係為根據本發明一實施例之方法其步驟之一流程表； 圖9與圖10係分別為使用本發明一實施例之方法估計圖1與圖5之模擬網路圖形； 圖11係顯示圖3程序其真實係數與使用本發明一實施例方法估計之平均係數兩者間差異之絕對值； 圖12係耦合白噪訊過程之一模擬15維度網路之圖形； 圖13係顯示僅使用同調性預測圖12之網路其連結與係數； 圖14係顯示僅使用部分同調性預測圖12之網路其連結與係數； 圖15係顯示使用本發明一實施例方法計算圖12之網路其連結與係數； 圖16係耦合白噪訊過程之一模擬40維度網路之圖形； 圖17顯示僅使用同調性重建圖16之預測網路； 圖18顯示使用部分同調性重建圖16之預測網路； 圖19顯示一健康志願者之腦波圖數據之交互結構； 圖20係以圖表方式描述一健康個體大腦所預期之網路類型； 圖21係以圖表方式描述一具有輕度認知損傷（mild cognitive impairment）個體其大腦所預期之網路類型。

Claims (14)

1. 一種監測一患者腦部功能之電腦執行方法，其包含：藉由一識別方法識別一腦波系統(electroencephalographic(EEG)system)之節點訊號間之網路連結，該腦波系統於一段時間內對於一患者執行一腦波記錄，其中，該識別方法，係於同時產生訊號之互動節點之一網路中，識別所述節點之連結以及估計經識別為相互連結之節點間之連結係數，該識別方法包含以下步驟：於一段預定時間內定期紀錄每一節點之訊號用以形成一資料集；計算該資料集中每一節點組合間之同調性(coherence)與部分同調性(partial coherence)；檢驗每一節點組合，假如該同調性或部分同調性其中之一係低於一第一預定閾值(threshold)，則對於所有後續步驟其對應之連結係數設置為零；一第一估計步驟，估計該資料集中連結係數尚未被設置為零之該節點組合之連結係數；對於經該所述第一估計步驟所估計低於一第二閾值之各連結係數，於所有後續步驟將所述之係數設置為零；以及一第二估計步驟，重新估計該資料集中連結係數尚未被設置為零之該節點組合之連結係數。
2. 如請求項1所述之電腦執行方法，其中，該檢驗步驟包含計算各節點組合其計算之同調性與部分同調性之乘積，以及決定所述乘積是否低於所述之第一預定閾值。
3. 如請求項1所述之電腦執行方法，其中，該所估計之連結係數之該第二閾值係由以下步驟決定： 根據該所估計之係數其歐幾里德距離(Euclidean distance)之平方將該所估計之係數分成兩群組，一包含該些具有高數值之係數之第一群組與一包含該些具有低數值之係數之第二群組；以及設置該所述第二閾值作為大於該所述第二群組中所有係數值之一係數值。
4. 如請求項1所述之電腦執行方法，其中，用於估計該連結係數之該第一與該第二估計步驟估計該資料集之自我迴歸係數(autoregressive coefficients)。
5. 如請求項1所述之電腦執行方法，其中，該網路係已知或經預測為稀疏連結(sparsely connected)。
6. 如請求項1所述之電腦執行方法，其中，該識別方法進一步包括，篩選該資料集用以移除離群值(outliers)之步驟。
7. 如請求項1所述之電腦執行方法，其中，該識別方法進一步包括，過濾該資料集用以移除噪訊(noise)之步驟。
8. 如請求項1所述之電腦執行方法，其包含：使用該識別方法識別腦波圖上節點訊號間之網路連結。
9. 一種用於識別網路連結與估計一腦部活動資料中節點間之連結係數之系統，該系統包含：複數感測器，於一段預定時間內記錄一個體於不同位置之腦部活動以產生一資料集；以及一處理器其組設用以：計算該資料集中各節點組合間之同調性與部分同調性； 檢驗各節點組合，假如該同調性或部分同調性其中之一係低於一第一預定閾值，則對於所有後續步驟所對應之連結係數設置為零；估計資料集中連結係數尚未被設置為零之該節點組合之連結係數；針對各經估計低於一第二閾值之連結係數，對於後續步驟設置該所述係數為零；重新估計資料集中連結係數尚未被設置為零之該節點組合之連結係數。
10. 如請求項9所述之系統，其中，該處理器係經設置用以計算該各節點組合其所計算之同調性與部分同調性之乘積，以及決定該所述之乘積是否低於該所述之第一預定閾值。
11. 如請求項9所述之系統，其中，該處理器係經設置用以藉由以下步驟決定該所述之第二閾值：根據該所估計之係數其歐幾里德距離之平方將該所估計之係數分成兩群組，分別係為包含該些具有高數值之係數之一第一群組，與包含該些具有低數值之係數之一第二群組；以及設置該所述第二閾值作為一大於該所述第二群組中所有係數值之係數值。
12. 如請求項9所述之系統，其中，該處理器係經設置用以估計該連結係數，藉由估計該資料集之自我迴歸係數。
13. 如請求項9所述之系統，其中，該處理器係經設置用以篩選該資料集以移除離群值。
14. 如請求項9所述之系統，其中，該處理器係經設置用以過濾該資料集以移除噪訊。

Images