李 強,潘 敏,成 波
(1.安康學院 數學與統(tǒng)計學院,陜西 安康 725000;2.西北大學 醫(yī)學大數據研究中心,陜西 西安 710127)
睡眠是諸多腦功能狀態(tài)中的一種,是指大腦意識相對喪失、自主的肌肉活動消失、受適當刺激可完全恢復到清醒時的生理功能的一種自主行為[1]。睡眠占據了人一生近1/3的時間,在生理功能調節(jié)、記憶鞏固、腦內代謝物清除以及生長發(fā)育等方面起著至關重要的作用。良好的睡眠有助于人們提高工作和學習效率,保持身心健康。根據R-K睡眠評估準則[2],睡眠狀態(tài)主要分為非快速眼動狀態(tài)(non-rapid eye movement,NREM)和快速眼動狀態(tài)(rapid eye movement,REM)。兩種狀態(tài)下大腦的放電活動是截然不同的,表現為腦電節(jié)律之間存在的顯著差異(見表1)。這些腦電節(jié)律的正常發(fā)放有助于維護睡眠穩(wěn)定性[3-5]。
表1 NREM狀態(tài)和REM狀態(tài)分別對應的典型腦電節(jié)律Tab.1 The typical EEG rhythms during NREM and REM
睡眠障礙是指睡眠量和質的異常,或在睡眠過程中發(fā)生某些臨床癥狀(如夢行癥)[6]。失眠是最常見且發(fā)病率極高的一種睡眠障礙,以入睡困難、睡眠維持困難及早醒等為主要特征,且每周至少發(fā)生3次并持續(xù)3個月以上[7]。盡管失眠在短期內不會直接威脅到人們的生命安全,但如果任由病情發(fā)展,會給患者帶來一系列并發(fā)疾病(如精神類疾病、心腦血管疾病、代謝性疾病等),甚至導致死亡。研究表明,臨床上一些失眠患者在睡眠REM階段,會出現明顯的腦電alpha節(jié)律活動增強現象[8]。這一發(fā)現對于通過調節(jié)alpha節(jié)律的發(fā)放來控制失眠患者的病情發(fā)展具有重要的指導意義。然而,關于這一現象發(fā)生的內在機制目前尚不清楚。
作為一類基于生物物理基礎的數學模型,神經計算模型(neural computational model,NCM)為探索、解釋及驗證大腦內各類生理病理現象的內在機制提供了一種有效途徑[9]。NCM主要采用數學理論與方法對大腦內在生理病理機制進行建模,進而通過神經環(huán)路重建的方式來模擬電信號的產生、演變及傳導等神經活動過程。神經集群模型是一類應用極為廣泛的宏觀神經計算模型,它將功能相似、結構相近的神經元看作為一個整體(即神經元集群),并通過描述神經元集群的平均特性來反映其整體放電行為[10-11]。典型的神經集群模型包括Jansen模型[12]、Wendling模型[13]、Liley模型[14]、Suffczynski模型[15]等。其中,Liley模型結構簡單且已成功應用于模擬NREM與REM睡眠狀態(tài)轉遷及對應的腦電節(jié)律[16-17]。
鑒于上述分析,本文將從計算神經科學的角度出發(fā),結合Liley神經計算模型與臨床失眠腦電數據,探索失眠患者在REM期典型放電節(jié)律的發(fā)生機制。首先,設計一種基于馬爾可夫蒙特卡洛的Liley模型關鍵參數的自動選擇與估計方法;其次,結合失眠患者腦電數據,確定影響alpha節(jié)律發(fā)放的模型關鍵參數,并據此給出失眠患者REM期腦電alpha節(jié)律活動增強的機制假設。
2002年,Liley等人通過考慮大腦皮層上神經元之間的興奮性與抑制性連接,提出了一個用于研究哺乳動物大腦電活動中alpha節(jié)律發(fā)放機制的神經計算模型(簡稱為Liley模型)。本節(jié)將對Liley模型的拓撲結構與數學建模過程進行概述。
Liley模型包括兩個神經元集群:興奮性神經元集群(E)與抑制性神經元集群(I),其拓撲結構如圖1所示。其中,集群E與I之間存在相互作用,且均存在自反饋作用以及來自丘腦的外部輸入(分別記為Pee、Pei)。圖1中黑色實線箭頭表示由興奮性神經遞質谷氨酸AMPA受體所介導的興奮性投射作用;黑色實線圓點表示由抑制性神經遞質γ氨基丁酸GABA受體所介導的抑制性投射作用。
圖1 Liley模型的拓撲結構Fig.1 The topological structure of Liley model
Liley模型的數學建模過程主要從以下3個方面來闡述。
1)突觸后膜電位的數學刻畫。
突觸后膜電位I(t)由突觸響應函數h(t)與突觸傳入的興奮性和抑制性脈沖(采用平均脈沖發(fā)放速率Q(t)刻畫)的卷積得到,即
(1)
式中:
h(t)=Γγ·te1-γt,t≥0
(2)
式中:Γ表示突觸后膜電位的最大幅值;γ表示突觸后膜電位的指數衰減時間尺度。
(3)
進一步,可轉化為如下二階微分方程來進行求解
(4)
2)胞體膜電位的數學刻畫。
采用突觸后膜電位I(t)的加權和刻畫胞體膜電位V(t),即
(5)
3)平均脈沖發(fā)放速率的數學刻畫。
由Sigmoid函數S(·)作用于胞體膜電位V(t),則可得平均脈沖發(fā)放速率Q(t),即
(6)
式中:Qmax為集群的最大脈沖發(fā)放速率;θ為發(fā)放閾值;σ為Sigmoid函數的梯度。
Liley模型中電活動的傳導過程主要由兩個計算模塊來完成:“平均脈沖發(fā)放速率-胞體膜電位”計算模塊(記為Q-V)與“胞體膜電位-平均脈沖發(fā)放速率”計算模塊(記為V-Q),如圖2所示。在第1個計算模塊中(藍色虛線),脈沖輸入Q(t)基于式(4)轉換為突觸后膜電位I(t),然后由式(5)得到胞體膜電位V(t);在第2個計算模塊中(紅色虛線),胞體膜電位V(t)基于式(6)轉換為平均脈沖發(fā)放速率Q(t)。
圖2 Liley模型中電活動的傳導過程Fig.2 The detailed transformation of electrical activities in Liley model
Liley模型的數學表達見式(7),模型輸出為興奮性神經元集群的胞體膜電位Ve。Liley模型中所有參數的生理學含義、參考取值及單位見表2[14,18-19]。
Γeγeee(NeeSe(Ve)+Pee(t)),
Γiγiee(NieSi(Vi)+Pie(t)),
Γeγeie(NeiSe(Ve)+Pei(t)),
Γiγiie(NiiSi(Vi)+Pii(t))
(7)
表2 Liley模型的參數取值及其生理學含義Tab.2 The parameters of Liley model and their physiological description
在Liley模型中,對模型行為影響較大的參數稱為模型關鍵參數。如何從眾多模型參數中自動選擇出關鍵參數并對其進行分析,顯然對于理解模型行為具有重要意義。本節(jié)給出一種數據驅動的模型關鍵參數的自動選擇與估計方法,具體步驟可總結為如下算法1。
算法1給定腦電信號S,記θ=(θ1,θ2,…,θN)為模型參數,其中N為參數的總個數。
記π(θj)(j=1,2,…,N)為第j個參數的先驗分布,采用馬爾科夫蒙特卡洛算法估計參數θj的后驗分布,即
(8)
其中
(9)
式(8)~(9)的更多細節(jié)可見文獻[19-20]。
步驟2選擇模型關鍵參數。首先,針對每一參數θj,計算其后驗樣本的均方根偏差
(10)
進一步,為了比較不同參數樣本內部的離散程度,計算參數θj后驗樣本的相對均方根偏差
(11)
式中:θjmax、θjmin分別表示參數θj可取的最大值與最小值。顯然,RRMSDθj的值越小,參數θj后驗分布范圍相對更集中,即從腦電數據S中所獲得的信息越多。
(12)
本文所使用的腦電數據來自于Kermanshah醫(yī)科大學睡眠障礙研究中心臨床采集的失眠腦電數據集[21]。該數據集的具體信息見表3。
上述數據集中,每30 s腦電片段的睡眠狀態(tài)評估結果均由睡眠專家根據AASM準則[22]進行了標注。針對健康受試者與失眠患者,分別隨機選取450個REM睡眠期的30 s單通道腦電片段,并計算其功率譜。所選腦電通道為F4[23]。圖3(a)和圖3(b)分別展示了所選健康受試者和失眠患者的腦電片段對應的頻譜圖。從圖中黑色橢圓圈標記處可以看出,相比于健康受試者,失眠患者腦電alpha節(jié)律所在頻段(8~13 Hz)對應的功率譜值較大,這也驗證了失眠患者REM期腦電alpha節(jié)律活動增強的結論。
圖3 不同受試者腦電片段的頻譜圖Fig.3 The EEG spectrograms of different subjects
表3 Kermanshah醫(yī)科大學睡眠障礙研究中心臨床采集的失眠腦電數據集的具體信息Tab.3 Detail information of EEG data from the Sleep Disorder Research Center in Kermanshah University of Medical Sciences
本小節(jié)基于上述所得腦電片段的功率譜,結合本文所提算法1,確定影響腦電alpha節(jié)律活動的模型關鍵參數及其最優(yōu)取值。
首先,以一個腦電片段為例對以上過程進行說明。圖4(a)展示了來自“失眠患者1”的REM期10 s腦電片段,主要表現為alpha節(jié)律的發(fā)放。采用該腦電片段估計Liley模型中22個參數的后驗分布,并對每個參數后驗分別計算相對均方根偏差(RRMSD)。進而對RRMSD按照從小到大的順序進行排序,并選擇排序前5%(這里將算法1中的參數q%設置為5%)的參數作為模型關鍵參數。圖5展示了所得的參數后驗分布(紅色)、先驗分布(藍色)及參數后驗樣本的RRMSD。從圖中紫色橢圓圈標記處可以看出,參數γi后驗分布范圍相對更小,對應的RRMSD值最小(RRMSDγi=0.012)。由此說明,參數γi從該數據中獲得了更多信息,對于調節(jié)alpha節(jié)律活動具有重要作用,即為模型關鍵參數。
其次,隨機選取80個來自失眠患者的REM期30 s腦電片段,計算所有模型參數后驗樣本的RRMSD值,所得結果如圖6所示??梢钥闯?相比于其他參數,參數γi后驗樣本的RRMSD值整體較小(圖中紅色框標記處所示)。由此可見,參數γi確為對腦電alpha節(jié)律活動影響最大的關鍵參數。結合其生理學含義,可以得出:抑制性突觸活動強度的變化可以影響失眠患者腦電alpha節(jié)律活動。
圖4 腦電片段及功率譜Fig.4 The EEG segments and power spectra
圖5 模型參數的后驗分布(紅色)、先驗分布(藍色)及參數后驗樣本的相對均方根偏差(RRMSD)Fig.5 The posterior (red color) and prior marginal (blue line) distributions of each model parameter and the RRMSD values of parameter posterior samples
圖6 80個腦電片段對應的參數后驗樣本的相對均方根偏差(RRMSD)Fig.6 The RRMSD values of parameter posterior samples obtained from 80 EEG segments
本小節(jié)面向健康受試者與失眠患者,分別根據其腦電片段估計對應的參數γi最優(yōu)取值,并據此分析比較二者之間的差異性,進而給出失眠患者腦電alpha節(jié)律活動增強的一種機制解釋。
隨機選取80個來自健康受試者的REM期30 s腦電片段以及80個來自失眠患者的REM期30 s腦電片段。圖7展示了由這些腦電片段所得到的參數γi最優(yōu)取值的箱線圖。從圖中可以看出,相比于健康受試者,由失眠患者腦電片段得到的參數最優(yōu)值整體偏高,即抑制性突觸后膜電位的衰減速率是增加的,進而抑制性突觸活動的強度是下降的。基于此,本文給出一種機制假設:抑制性突觸活動強度的降低可誘導失眠患者REM期腦電alpha節(jié)律活動的增強。
圖7 參數γi最優(yōu)估計值的箱線圖(由健康受試者與失眠患者腦電數據確定)Fig.7 Box plots of the optimized values of parameter γi(obtained from EEG data of healthy subjects and patients with insomnia)
本文從計算神經科學的角度出發(fā),結合Liley神經計算模型與臨床失眠腦電數據,對失眠患者REM期腦電alpha節(jié)律活動增強現象進行了內在機制上的解釋。具體地,設計了一種基于馬爾可夫蒙特卡洛的Liley模型關鍵參數的自動選擇與估計方法;采用臨床采集的失眠患者腦電數據,通過數值實驗驗證了所提方法的有效性,確定了影響alpha節(jié)律發(fā)放的模型關鍵參數,并據此給出了失眠患者REM期腦電alpha節(jié)律活動增強的機制假設:抑制性突觸活動強度的降低會誘導alpha節(jié)律活動的增強。所得結論對于通過調節(jié)alpha節(jié)律的發(fā)放來改善失眠患者睡眠質量并控制病情發(fā)展,具有一定的指導作用和臨床價值。