低氧環(huán)境下腦網(wǎng)絡的建模與分析

王 磊，韋 新，李雪鵬，楊 昱，李海芳

（1.太原理大學信息與計算機學院，山西 太原 030024；2.西安交通大學社會心理學研究所，陜西 西安 710049）

如今，前往高海拔低氧地區(qū)工作、學習或者居住的人數(shù)持續(xù)增長[1]。暴露于高海拔低氧環(huán)境下的人群通常表現(xiàn)出知覺和認知領(lǐng)域的缺陷，例如工作記憶、視覺、注意力、執(zhí)行功能等[2-3]。值得注意的是，有關(guān)高海拔暴露引起認知功能變化的研究大多來自急進海拔[4]，但對長期低氧環(huán)境下人的認知功能腦網(wǎng)絡變化研究尚不完善，對其腦網(wǎng)絡的動態(tài)演化特征更是知之甚少。

在網(wǎng)絡研究中，通過圖論方法來表征復雜網(wǎng)絡的拓撲關(guān)系是研究網(wǎng)絡特性的重要手段。使用功能連接來研究大腦進行認知活動時的腦網(wǎng)絡特征是非常有效的方式。研究腦網(wǎng)絡拓撲屬性有助于研究者更好地考察大腦在進行認知活動時的功能分離和整合等特性，這些特性和大腦功能密切相關(guān)。另外，功能連接可以通過后期的認知訓練進行重組[5-7]，重塑回健康大腦。因此，確定高原暴露對認知功能的影響，從而闡明潛在神經(jīng)機制具有重要意義。

大腦當前的動力狀態(tài)是由特定的大腦編碼狀態(tài)決定的[8]，特定狀態(tài)的編碼將對后續(xù)狀態(tài)具有預測能力[9]，可以通過預測編碼動力學機制中的異常來探究海拔對人腦網(wǎng)絡結(jié)構(gòu)的改變?；诖嘶A提出了腦網(wǎng)絡重組預測模型，該模型可以模擬被試認知任務腦網(wǎng)絡的變化過程，用于揭示腦網(wǎng)絡動力學異常。

該研究采集了3個不同海拔任務態(tài)EEG 數(shù)據(jù)，以通道為節(jié)點，通道間的鎖相值（Phase Locking Value，PLV）為邊，應用滑動窗口技術(shù)構(gòu)建了隨時間變化的動態(tài)PLV 腦網(wǎng)絡，計算了每個被試動態(tài)網(wǎng)絡的拓撲參數(shù)，分析了腦網(wǎng)絡的動態(tài)變化過程，并挑選出海拔組間差異性顯著的時間窗口評估預測模型的準確性。高海拔數(shù)據(jù)集網(wǎng)絡權(quán)值分布不均，為了消除個體權(quán)值的影響和最大程度擬合腦網(wǎng)絡變化過程，設計了6 種不同的權(quán)值更新概率模型，并進行了實驗驗證。α波活動減弱是低氧環(huán)境下的一個重要性表現(xiàn)[10]，因此網(wǎng)絡重組預測模型對α頻段刺激前的腦網(wǎng)絡進行迭代預測，為每個被試挑出擬合度最高的網(wǎng)絡模型。通過分析每個海拔動態(tài)網(wǎng)絡的拓撲參數(shù)差異性和網(wǎng)絡模型選擇的差異性，發(fā)現(xiàn)了高原久居者網(wǎng)絡損傷的異常表現(xiàn)，進而為高原損傷的認知康復訓練提供了重要參考。

1 被試與實驗

1.1 被 試

研究使用的數(shù)據(jù)來自西安交通大學韋新教授團隊采集的高原腦電數(shù)據(jù)[11]。實驗招募了75 名健康男性青年，每組25 人，平均年齡為20.27±2.56 歲，平原出生，進駐高原前從未上過高原，分別移居渭南（400 m）、拉薩（3 700 m）、那曲（4 500 m）。所有受試者在高原居住平均時長為2.12±0.31 年。所有被試均為右利手，無精神障礙。實驗前24 小時內(nèi)未飲用咖啡和酒，無服用精神興奮類藥物，睡眠正常。所有被試均已填寫知情同意書。

1.2 Stroop任務

刺激材料由兩個中文顏色字(紅和綠)和兩種顏色(紅色和綠色)兩兩搭配組成兩種刺激。每個任務中一致刺激為38 個，不一致刺激為102 個，以偽隨機順序呈現(xiàn)。顏色字為48號仿宋體，呈現(xiàn)于計算機14寸顯示器屏幕中央,視角為2.2°。實驗由一個顏色判斷任務和一個字義判斷任務組成。被試通過按鍵作出判斷，正式實驗前讓被試充分練習直到完全理解。記錄準確率和平均命中反應時間進行分析。

1）顏色任務：要求被試克服字義干擾，觀察字體顏色，如果字體顏色為紅色，按“z”鍵，字體顏色為綠色，按“m”鍵。對于每組實驗，首先在屏幕中央呈現(xiàn)一個“+”作為注視點，呈現(xiàn)時間固定為300 ms。顏色刺激呈現(xiàn)時間固定為200 ms，被試在2 000 ms 內(nèi)作出按鍵判斷，接著出現(xiàn)下一個顏色字等待被試按鍵反應，如果被試按鍵錯誤，2 200 ms 后會進入下一輪相同的任務。在任務中要求被試盡可能又快又準地反應。

2）字義任務：與顏色任務相似，要求被試克服字體顏色的干擾，判斷呈現(xiàn)刺激的字義。

圖1 Stroop范式任務示意圖

2 方 法

2.1 EEG數(shù)據(jù)采集與預處理

數(shù)據(jù)采集使用的是Neuroscan 公司的36 導腦電系統(tǒng)，電極分布符合國際10-20 系統(tǒng)標準，參考電極置于耳垂處，記為A1和A2，放置4個記錄電極分別位于眼睛上下部和眼外眥處，用于采集垂直方向眼電和水平方向眼電。采集EEG 數(shù)據(jù)時使用的采樣頻率為1 000 Hz，頭皮與電極之間的阻抗保證低于5 kΩ。

實驗所有的數(shù)據(jù)分析使用的是Matlab2017b 和開源工具箱EEGLAB 進行了預處理，將4 個眼電電極去除。對EEG 數(shù)據(jù)進行了常規(guī)去除噪音預處理，提取所有不一致刺激開始前200 ms 到刺激開始后600 ms 為一個試次，確保與隨后的試驗不重疊，并修正了每個時期的基線。

2.2 腦網(wǎng)絡的構(gòu)建

神經(jīng)科學研究認為，發(fā)生鎖相的神經(jīng)元群體間可以進行有效的信息通信[12]，而PLV可量化在特定頻帶和時間區(qū)域的兩個神經(jīng)信號進入鎖相狀態(tài)的同步程度[13]。因此該研究以頭皮電極為節(jié)點，通道間的PLV值為連邊，構(gòu)建加權(quán)網(wǎng)絡。PLV可通過如下公式求得：

式（1）中，Δφn(t)=φx(t)-φy(t)表示信號x與信號y在t時刻的瞬時相位差，N為時間點。PLV值介于0 到1之間，其值為1 時，表示兩個信號相位完全同步；值為0 時，表示相位不同步，無相位鎖定。

為了評估認知任務期間腦網(wǎng)絡隨時間的變化，研究采用滑動窗口來評估認知任務期間腦網(wǎng)絡隨時間的變化。對每個被試的每個刺激前200 ms 到刺激后600 ms 的數(shù)據(jù)，采用了重疊度為90%，窗口大小為200 ms 的滑動窗口構(gòu)建腦網(wǎng)絡，為了消除特異性，對每個被試的所有刺激網(wǎng)絡進行了平均，并提取了每個被試的網(wǎng)絡屬性[14]。

2.3 網(wǎng)絡屬性

腦功能網(wǎng)絡的拓撲性質(zhì)可以通過一些拓撲結(jié)構(gòu)參數(shù)來表征，該文采用了特征路徑長度、聚類系數(shù)[15]、網(wǎng)絡密度和全局效率[16]對腦網(wǎng)絡進行評估。

特征路徑長度可以表征網(wǎng)絡的集成程度[16]，它定義為網(wǎng)絡中所有節(jié)點對之間的平均最短路徑長度：

式（2）中，dij為節(jié)點i與節(jié)點j之間的最短路徑長度，N為網(wǎng)絡節(jié)點個數(shù)。

平均聚類系數(shù)可以量化網(wǎng)絡的分離度。在加權(quán)網(wǎng)絡的情況下，為了避免主連接權(quán)值的影響，平均聚類系數(shù)可以推廣如下[17]：

式（3）中，Wij表示電極i和j之間的權(quán)值，N為網(wǎng)絡節(jié)點個數(shù)。

網(wǎng)絡密度可以評估網(wǎng)絡的連接強度，其計算公式如下：

式（4）中，L為網(wǎng)絡所有權(quán)值之和，N為網(wǎng)絡節(jié)點個數(shù)。

全局效率可衡量網(wǎng)絡的信息流動快慢[16]，全局效率的定義公式：

式（5）中，dij表示節(jié)點i到節(jié)點j的隨機距離，N為網(wǎng)絡節(jié)點個數(shù)。

2.4 技術(shù)路線

整體實驗技術(shù)路線如圖2 所示，首先選取刺激前200 ms 到刺激后600 ms 的EEG 數(shù)據(jù)，采用窗口大小為200 ms，重疊率為90%的滑動窗口為每個被試構(gòu)建了31 個連續(xù)的動態(tài)PLV 腦網(wǎng)絡，然后分別對每個腦網(wǎng)絡提取網(wǎng)絡屬性，挑選出刺激后3 個海拔組間網(wǎng)絡屬性差異性最大的窗口作為衡量模型預測的準確度。

圖2 技術(shù)路線圖

最后分別對每個被試刺激前200 ms 的EEG 數(shù)據(jù)構(gòu)建PLV 腦網(wǎng)絡，為了確定每個被試的腦網(wǎng)絡重組預測模型，根據(jù)不同的網(wǎng)絡預測模型規(guī)則，對每個被試刺激前的腦網(wǎng)絡權(quán)值進行概率修改，為每個被試挑選出模擬網(wǎng)絡屬性和真實網(wǎng)絡屬性之間均方差（MSE）最小的模型類型。

2.5 網(wǎng)絡預測模型

6 種網(wǎng)絡預測模型如下所示：

1）主連接增強：該模型假設刺激前腦網(wǎng)絡的主要連接（即網(wǎng)絡權(quán)值較大的連接），在認知任務過程中權(quán)值更有可能增加，增加概率為p=P；

2）次連接增強：該模型假設刺激前腦網(wǎng)絡的次要連接（即網(wǎng)絡權(quán)值較低的連接）在認知任務過程中權(quán)值更有可能增加，增加概率為p=(1-P)；

3）隨機增強：該模型假設刺激前腦網(wǎng)絡連接是隨機增加。即所選連接的權(quán)值增加概率p=1；

4）主連接削弱：該模型假設刺激前腦網(wǎng)絡的主要連接（即網(wǎng)絡權(quán)值較大的連接）在認知任務過程中權(quán)值更有可能減小，減小概率為p=P；

5）次連接削弱：該模型假設刺激前腦網(wǎng)絡的次要連接（即網(wǎng)絡權(quán)值較低的連接）在認知任務過程中權(quán)值更有可能減小，減小概率為p=(1-P)；

6）隨機削弱：該模型假設刺激前腦網(wǎng)絡連接是隨機減小。即所選連接的權(quán)值減小概率p=1。

算法步驟總結(jié)如下：

1）隨機選取刺激前腦網(wǎng)絡連通矩陣的一個連接（圖論中通常稱為Wij）；

2）連接值Wij根據(jù)概率p進行權(quán)值修改，增加或減小的大小為（1%Wij）。概率p的值和如何修改取決于當前所考慮的特定模型；

3）計算每個模型修改后腦網(wǎng)絡的網(wǎng)絡屬性；

4）計算3）中腦網(wǎng)絡的網(wǎng)絡屬性與刺激后腦網(wǎng)絡的網(wǎng)絡屬性之間的MSE；

5）重復5 000 次步驟1）、2）、3）和4），并為每次迭代存儲MSE；

6）選擇MSE最小的預測模型。

為了最小化算法內(nèi)在隨機行為產(chǎn)生的偏差，所有實驗重復100 次，將預測模型與真實網(wǎng)絡的屬性之間的MSE最小的作為選定模型。

如圖3 所示，高海拔被試網(wǎng)絡權(quán)值呈正態(tài)分布，并且權(quán)值分布比較集中，被試之間的網(wǎng)絡權(quán)值分布中心也有差異。因此為了消除個體網(wǎng)絡權(quán)值分布集中問題，突出網(wǎng)絡權(quán)值的主次性，并且體現(xiàn)權(quán)值更新的概率性問題，網(wǎng)絡重組預測模型更新概率P的設計尤為重要，該模型設P=num/N，num是所選連接的排序位置，N是網(wǎng)絡總連接數(shù)。此概率設置可以消除網(wǎng)絡權(quán)值分布問題，凸顯各模型之間的差異性。

圖3 被試PLV權(quán)值分布圖

3 實驗結(jié)果與分析

3.1 行為數(shù)據(jù)結(jié)果與分析

對Stroop 任務的正確率和反應時間進行了研究，同時利用信息傳輸率（Information Translate Rate,ITR）量化了被試的大腦信息加工速率，從而評估不同海拔人群大腦認知加工速率的差異性，實驗結(jié)果如表1 示。

從表1中可以看出每個海拔詞義判斷任務的ITR值都高于顏色判斷任務的ITR值，這是因為在Stroop任務中讀詞是行為習慣的自動加工過程，顏色命名是意識控制較強的行為加工過程，自動加工不需要太多注意力的參與，意識控制加工過程則需要有意的控制，相對字義任務，顏色判斷任務的復雜強度更大[18]。

表1 Stroop任務行為結(jié)果統(tǒng)計分析

統(tǒng)計分析發(fā)現(xiàn)，在顏色和詞義任務中，反應時間在各組間存在顯著性差異(P色=0.01，F(xiàn)色=4.912；P詞=0.01，F(xiàn)詞=5.375)，并且反應時間隨著海拔的增高而延長。經(jīng)Bonferroni 校正后，在顏色任務中，ITR 在渭南和那曲組間有顯著性差異(P=0.045)。而正確率在各海拔之間沒有顯著性差異。長期居住在高海拔低壓低氧的環(huán)境下，人的反應速度會相對變慢，而正確率相對平衡[13]，這可能是由于高海拔環(huán)境下大腦長期處于供氧不足，疲憊狀態(tài)，導致認知能力的神經(jīng)元受損[14-15]，用時間上的消耗來補償認知反應。

3.2 網(wǎng)絡動力學結(jié)果與分析

為了評估網(wǎng)絡屬性的時變性，計算了每個動窗口下的網(wǎng)絡屬性。圖4 給出了顏色任務α頻段下不同海拔下網(wǎng)絡屬性的動力學變化。通過特征路徑長度、聚類系數(shù)、網(wǎng)絡密度和全局效率可以觀察到高海拔組相對低海拔組腦網(wǎng)絡表現(xiàn)出更低整合性、更低分離性、更低的連接強度和更低的傳播效率。以往研究表明在300 ms 左右高海拔與低海拔組腦電屬性具有顯著性差異[17-18]，該研究對第20 個窗(200～400 ms)下的腦功能網(wǎng)絡進行了方差檢驗，特征路徑長度(F=5.187,P＜0.01)，聚類系數(shù)(F=5.189,P＜0.01)，網(wǎng)絡密度(F=4.505,P＜0.05)，全局效率(F=5.021,P＜0.01)在組間均具有顯著性差異。此窗口下海拔與特征路徑長度(r=0.336 093,P=0.009)，聚類系數(shù)(r=-0.336 065,P=0.009),網(wǎng)絡密度(r=-0.338 507,P=0.008)和全局效率(r=-0.332 180,P=0.009 514)也具有一定的相關(guān)性。腦功能網(wǎng)絡的集成和連接強度的降低表明了腦網(wǎng)絡之間的通信能力降低，這與先前的研究神經(jīng)活動普遍降低的相結(jié)果一致[2]。功能性腦網(wǎng)絡分離性降低表明大腦局部性溝通較弱或較少分離的神經(jīng)處理[9]。網(wǎng)絡全局效率的降低表明網(wǎng)絡的信息處理能力下降。功能性腦網(wǎng)絡分離性降低表明大腦局部性溝通較弱或較少分離的神經(jīng)處理[9]。網(wǎng)絡全局效率的降低表明網(wǎng)絡的信息處理能力下降。長期的高原暴露會導致人的神經(jīng)活動減弱，進而反映在認知反應速度上的減慢，但是神經(jīng)活動減弱的原因尚不確定是神經(jīng)損傷所致還是供氧不足所致。在以后的研究中可以采集遣返人員腦電信號作進一步分析。

圖4 網(wǎng)絡屬性隨時間變化

3.3 網(wǎng)絡預測模型分析

根據(jù)網(wǎng)絡動力學的結(jié)果，選取刺激后200～400 ms(第20 個窗口)數(shù)據(jù)的網(wǎng)絡屬性作為預測模型的判斷標準。將每個被試刺激前200 ms 數(shù)據(jù)的腦網(wǎng)絡放入預測模型中進行迭代預測，為每一個被試挑選出最能擬合刺激前網(wǎng)絡變化的認知反應模型。如圖5 所示預測模型與認知反應之間的均方誤差(MSE)達到了98.951 2%，這證明了改進模型的準確性。

圖5 模型結(jié)果分析圖

表2 預測模型迭代次數(shù)和準確率

顏色任務下每組的模型選擇分布如圖6(a)所示，高海拔和低海拔被試在認知過程中次連接（即刺激感知前大腦區(qū)域之間低相位同步連接）的增強為主要的腦網(wǎng)絡動力模型，然而，海拔組間模型選擇分布有統(tǒng)計學意義(χ2(3，N=75)=11.04，P＜0.0 5；卡方檢驗)。從刺激前到刺激后網(wǎng)絡參數(shù)的變化量與網(wǎng)絡重組模型的選擇之間存在統(tǒng)計學的顯著關(guān)系((F渭=5.876，P渭=0.009)，(F拉=10.251，P拉=0.001)，(F那=85.903，P那＜0.001)；方差分析)。圖6（b）中給出了次連接增強模型與認知任務中更高的網(wǎng)絡重構(gòu)相關(guān)聯(lián)[19-20]，這可被視為“正常行為”。然而在高海拔組中，部分被試刺激前后網(wǎng)絡參數(shù)差異性較小，網(wǎng)絡重構(gòu)能力更弱，這也證明了長期的高海拔暴露可能會引起認知腦功能網(wǎng)絡的異常。

圖6 模型選擇結(jié)果圖

低壓低氧會誘導灰質(zhì)體積和白質(zhì)纖維束功能的改變，與腦功能區(qū)的異?；顒佑嘘P(guān)，高海拔引起的海馬損傷有隨海拔和暴露時間的增加而加重的趨勢[21]。因此在高海拔組中有更多的被試偏離了次連接增強的模型，這可能是由于腦功能區(qū)的損傷和激活程度的降低，使得全腦網(wǎng)絡差異性減弱[22-23]。

4 結(jié)束語

文中從滑動窗口的角度出發(fā)構(gòu)建腦功能網(wǎng)絡，分析腦網(wǎng)絡屬性的動態(tài)變化過程，找出差異顯著的時間窗口，并利用基于邊概率演化的腦網(wǎng)絡預測模型進行模擬腦網(wǎng)絡變化過程。通過實驗可以得出以下結(jié)論：

1）通過滑動窗口發(fā)現(xiàn)海拔對人腦認知影響主要發(fā)生在300 ms 左右，并且隨著海拔高度的增加，大腦的神經(jīng)活動減弱，腦網(wǎng)絡表現(xiàn)出更低整合性、更低分離性、更低的連接強度和更低的傳播效率。

2）腦網(wǎng)絡預測模型結(jié)果表明長期高海拔暴露會引起大腦網(wǎng)絡動力減弱，同時損害腦功能網(wǎng)絡的重構(gòu)能力，而且大腦的受損程度會隨著海拔高度的增加而加重。