大腦神經元過度的興奮可以導致癲癇的發作, 維持興奮-抑制平衡是癲癇控制的關鍵。提出了癇性指數來描述癲癇發作程度, 并用作PID控制器的被控參數來對癲癇發作進行控制。把神經群模型作為平臺來仿真興奮性增加導致的癲癇發作程度變化, 進而對兩種癲癇控制策略進行了仿真。實驗結果表明興奮強度增加而保持抑制強度不變會導致癇性指數的大幅度增加, 導致癲癇。而用PID控制器分別降低興奮強度或增加抑制強度都可以維持興奮-抑制平衡, 并緩解癲癇的發作。癇性指數可以描述腦電信號的線性和非線性特性, PID控制器簡單而且不依賴于潛在的生理結構, 為本方法應用于臨床打下了基礎。
引用本文： 馬震. 基于神經群模型的致癇興奮性控制研究. 生物醫學工程學雜志, 2016, 33(2): 255-259. doi: 10.7507/1001-5515.20160044 復制
引言
癲癇是一種常見的慢性腦部疾病，嚴重危害人類健康。腦電圖(electroencephalogram, EEG)是腦神經細胞電活動在大腦皮層或頭皮表面的總體反映，是研究癲癇發作特征的重要工具[1]。現在根據EEG進行的癲癇診斷、治療往往是醫療工作者依據經驗人工完成，是一項復雜而又耗時的工作。因此，計算機輔助癲癇控制是一種極有價值的輔助治療手段[2]。
當神經群的興奮性增加，打破正常的興奮-抑制平衡關系，將會引起神經元同步放電，進而誘發癲癇。所以如何控制誘發癲癇的過度興奮成為癲癇控制的關鍵。神經群模型包括能夠與癲癇重要機制直接聯系的參數，例如興奮、抑制或者連接[3]，可以很好地揭示癲癇發作的生理機制，為研究過度興奮的控制提供了平臺。皮層柱模型是Jansen等[4-5]提出的，被用來研究自發EEG信號和視覺皮層的誘發電位。這個模型被后來的研究者改進[3, 6-10]，并應用在不同領域的研究中。雖然改進的模型可以產生更加豐富的波形，但是這些模型中的某些參數并不具有明確的生理學意義。為了更好地探索癲癇發作的生理學機制，我們采用了模型參數具有生理學意義的Jansen模型。
Wang等[11]提出了一種神經群模型線性化的方法，經過線性化的神經群模型可以用一個傳輸函數進行表示。癲癇發作被看作神經群模型不穩定所導致，并給出了一種PD控制器來保持神經群的穩定性，抑制癲癇發作[11]。Ma等[12]針對同步行為的可傳遞性，提出了面向神經群族的控制方案，可以有效地通過打破神經群族的同步性來防止癲癇發作。
在本文中，我們提出了一種致癇興奮性控制方法。首先，給出癇性指數的定義，并用來描述癲癇發作的程度。依據腦電信號的癇性指數來產生控制信號，從而對癲癇興奮性進行控制，阻止癲癇發作。
1 基于神經生理學的神經群數學模型
神經群模型的結構如圖 1所示。在該模型中，錐體神經元和中間神經元通過突觸互相連接。圖中，實線表示的是大腦內部信號流，虛線表示外部控制信號流。
圖1
神經群模型
Figure1.
Neural mass model
本地神經群模型中，一個線性系統可以將興奮或者抑制興奮神經子群的平均行動電位發放率m(t)轉換為平均后突觸電位v(t)，也就是v(t)=he, i(t)*m(t)，可以表示為：
| $ {{h}_{\text{e, i}}}\left(t \right)=\frac{{{A}_{\text{e, i}}}}{{{\tau }_{\text{e, i}}}}t{{e}^{-\frac{1}{{{\tau }_{\text{e, i}}}}}}, $ |
式(1)中下標e和i分別表示興奮和抑制神經子群。Ae, i描述了興奮/抑制子群后突觸電位的特性，可以分別用來調節興奮和抑制突觸的強度，τe, i為時間常數。本文中，Ae, i由神經群本身的興奮/抑制程度和兩個外部控制器(EPIDC和IPIDC)的影響共同決定，也就是
| $ {{A}_{\text{e, i}}}=A_{\text{e, i}}^{n}-A_{\text{e, i}}^{c}, $ |
其中，Ae, in為神經群自身的興奮/抑制強度，Ae, ic為來自于控制器EPIDC或IPIDC控制信號。
非線性環節S(v)將一個神經元子群的平均膜電位變換為動作電位的平均發放率，這個變換用如下形式的Sigmoid函數來描述，
| $ S\left(v \right)=2{{e}_{0}}/\left(1+{{e}^{r\left({{v}_{0}}-v \right)}} \right), $ |
這里2e0是最大發放率，v0是對應于發放率e0的后突觸電位，r是Sigmoid函數的陡峭度。外部興奮輸入P(t)代表了附近群的影響，在本文中通過正均值的高斯白噪聲來表示，其均值為90，方差為30。錐體神經元和中間神經元之間興奮和抑制突觸的數目是由C1~C4四個連接參數來描述[4-5, 12]。
2 癲癇控制
2.1 特征提取
Pincus[13]在研究心率變化時提出了近似熵(approximate entropy，ApEn)的算法，越復雜的時間序列對應的近似熵越大。而當癲癇發作時，腦電信號的復雜度會降低，因此，本文采用Olson等[14]的方法來計算近似熵，并作為描述腦電信號復雜性的非線性指標。
此外，我們采用小波變換技術對腦電信號進行多尺度分解，提取對應頻段在3~29 Hz范圍內各尺度上重構信號的波動指數[15]作為腦電數據的線性特征。通常情況下，癲癇發作時腦電信號的波動會比未發作的間歇期劇烈，波動指數可以用來衡量信號變化的強度，其定義如下:
| $ {{F}_{i}}=\sum\limits_{j=1}^{M-1}{\left| {{a}_{i}}\left(j+1 \right)-{{a}_{i}}\left(j \right) \right|}, $ |
式中: ai為第i尺度上小波重構信號的幅度，M為信號的長度。
我們綜合二者，提出癇性指數(epileptiform index, Ie)來描述腦電信號的發作程度，定義如下：
| $ {{I}_{c}}=\sum\limits_{i\in N}{{{F}_{i}}/{{E}_{a}}}, $ |
這里，Fi為小波變換后第i尺度的的波動指數，N為納入考慮的尺度集合，在本文中，選取對應頻段在3~29 Hz范圍內的第3、4、5尺度，即N={3, 4, 5}，并分別計算其波動指數。Ea為癲癇信號的近似熵。
2.2 控制方法
我們所提出的控制策略的總體目標是通過降低興奮強度或者增加抑制強度來使得興奮和抑制處于平衡狀態，從而抑制發作。所以，我們引入兩個PID控制器EPIDC和IPIDC，其中EPIDC用來控制興奮線性模塊，產生控制信號Aec；IPIDC用來控制抑制線性模塊，產生控制信號Aic。兩個控制器的控制信號是癇性指數的函數。
首先，計算腦電信號的近似熵Ea，并進行小波變換，計算各尺度的波動指數Fi；然后，計算癇性指數Ie，如果大于閾值TI則認為興奮性過高，足以引發癲癇，控制器開始工作；最后，根據計算癇性指數Ie確定控制量，控制量初始置為0，如果需要進行調節，則更新控制量，否則保持控制量不變。控制量為：
| $ A_{\text{e, i}}^{c}\left(t \right)={{k}_{\text{p}}}*e\left(t \right)+{{k}_{\text{i}}}*\int{e\left(t \right)+{{k}_{\text{d}}}}*\frac{\text{d}e\left(t \right)}{\text{d}t}, $ |
其中，e(t)=Ie-TI。對于EPIDC，kp、ki、kd三個參數為正數，而對于IPIDC，這三個參數為負數。
3 實驗結果
構建包括一個神經群及兩個PID控制器的實驗系統，調節神經群的參數保持在標準值[5]，控制器保持在待命狀態以模擬正常的腦電行為。然后，通過改變Aen使得神經群的興奮性產生變化，來驗證本文所提出的控制器是否能有效地控制致癇興奮性。
3.1 無控制器時的發作情況
為了研究興奮性的增加對癇性指數的影響，保持Ain不變，Aen從2增加到11，步長為0.01，對應的癇性指數變化曲線如圖 2所示。從圖 2可以看到，在沒有控制器介入的情況下，癇性指數會隨Aen的增加而增大。當Aen=7.71時，癇性指數達到最大值34.350 7。當Aen在2~3.45范圍內變化，癇性指數一直保持增加，從0.684 5增加到1.350 3，但增加幅度不大。而當Aen=3.46時，癇性指數突然增加到4.225 1。隨著Aen繼續增加，癇性指數一直保持在較高值。當在Aen=7.71達到最大值后，癇性指數有一定程度下降，但是依然保持在25左右。
圖2
無控制器時癇性指數隨興奮性增加的變化
Figure2.
Epileptiform index as a function of excitatory strength in the case of no controller
3.2 有控制器時的發作情況
為了檢驗本文提出的控制器是否能較好地調節過度興奮引發的癲癇發作，通過增加興奮強度Aen來仿真神經群興奮性增加的情況，然后應用本文的控制器來對興奮性的增加進行控制。具體為Ain保持標準值，Aen從2增加到8，然后由8減少到2，步長為0.1。仿真時間為122 s，每秒鐘計算一次癇性指數Ie，如果超過閾值TI則控制器介入進行調節，TI設置為1.5，這里閾值取不同的值會使得控制器有不同的靈敏性，介入的時刻也不同。在仿真甚至實際應用中，TI可以根據患者的情況進行設置，使控制器可以更好地適應不同的病情。為了分別討論兩種控制策略(降低興奮強度和增加抑制強度)對癲癇發作的調節效果，EPIDC和IPIDC分別應用于神經群模型。EPIDC采用的kp、ki、kd參數值分別為0.11、0.003 5和0.002 4，IPIDC的kp、ki、kd參數值分別為-1.82、-0.122和0.116。
只有EPIDC工作的時候，Aen和Ace如圖 3(a)所示。從圖中可以看出，通過PID控制器，Ace能根據Aen的變化而有效地進行改變，給出的控制反饋可以抵消過度興奮所造成的影響，使興奮強度保持在正常水平。從圖 3(b)可以發現，采用這種方法保持了興奮-抑制平衡，癇性指數均值為1.961 1，與無控制器的情況相比有大幅度下降，保持在正常范圍，避免了發作。在A通過臨界點3.5的時候，無論是增加過程還是下降過程，都會引起癇性指數的突變，但控制器都可以很快地隨之產生相應的控制信號，從而有效地降低興奮強度。
圖3
只有EPIDC起作用時的控制結果
(a)
(a)
只有IPIDC工作的時候，Aen和Aic如圖 4(a)所示。在圖中，因為Aen和Aic數值差距較大，為了對比更加明顯，將Aen放大十倍進行顯示。從圖中可以看出，與EPIDC通過產生控制信號直接抵消興奮強度的增加不同，IPIDC根據癇性指數的變化產生控制信號來增加抑制強度，從而維持興奮-抑制平衡。從圖 4(b)可以看出，癇性指數并未像無控制器情況下隨興奮強度增加而迅速增加，均值為1.760 3。在仿真過程中，在興奮強度增加通過臨界點A=3.5的時候，會引發癇性指數的突然增加，控制器可以很快地產生增加抑制強度的控制信號，有效地控制了癲癇發作。這種策略雖然不通過直接降低興奮強度來維持興奮-抑制平衡，但是增加抑制強度也可以抵消興奮強度增加帶來的影響，保持二者之間的平衡。
圖4
只有IPIDC起作用時的控制結果
(a)
(a)
4 結論
當癲癇發作時，腦電信號從高復雜度過渡到低復雜度，本文用近似熵來描述腦電信號的復雜度，并作為表征癲癇發作程度的非線性指標。波動指數可以反映腦電信號的波動劇烈程度，而在癲癇發作時，腦電信號波動加劇，所以本文采用波動指數作為線性指標來描述癲癇發作程度。進而，提出可以反映腦電信號線性和非線性特征的癇性指數，并用來作為PID控制器的被控參數來維持大腦的興奮-抑制平衡。
采用神經群模型作為研究平臺，對興奮性增加引起的癇性指數的變化進行了討論。如果抑制強度保持不變，隨著興奮強度的增加癇性指數也會增加。當興奮強度超過某一臨界值，癇性指數會突然大幅度增加，導致癲癇發作。本文采用兩個PID控制器來模擬兩種控制興奮性增加的策略——降低興奮強度和增加抑制強度。仿真結果證明這兩種策略都可以很好地控制興奮性增加導致的癇性發作。本文所提出的控制方法為設計外部發作控制設備提供了理論的基礎。PID控制器結構簡單而且不依賴于潛在的生理結構，有利于將其應用于臨床。將此控制方法應用于動物模型將是我們下一步研究的重點。
引言
癲癇是一種常見的慢性腦部疾病，嚴重危害人類健康。腦電圖(electroencephalogram, EEG)是腦神經細胞電活動在大腦皮層或頭皮表面的總體反映，是研究癲癇發作特征的重要工具[1]。現在根據EEG進行的癲癇診斷、治療往往是醫療工作者依據經驗人工完成，是一項復雜而又耗時的工作。因此，計算機輔助癲癇控制是一種極有價值的輔助治療手段[2]。
當神經群的興奮性增加，打破正常的興奮-抑制平衡關系，將會引起神經元同步放電，進而誘發癲癇。所以如何控制誘發癲癇的過度興奮成為癲癇控制的關鍵。神經群模型包括能夠與癲癇重要機制直接聯系的參數，例如興奮、抑制或者連接[3]，可以很好地揭示癲癇發作的生理機制，為研究過度興奮的控制提供了平臺。皮層柱模型是Jansen等[4-5]提出的，被用來研究自發EEG信號和視覺皮層的誘發電位。這個模型被后來的研究者改進[3, 6-10]，并應用在不同領域的研究中。雖然改進的模型可以產生更加豐富的波形，但是這些模型中的某些參數并不具有明確的生理學意義。為了更好地探索癲癇發作的生理學機制，我們采用了模型參數具有生理學意義的Jansen模型。
Wang等[11]提出了一種神經群模型線性化的方法，經過線性化的神經群模型可以用一個傳輸函數進行表示。癲癇發作被看作神經群模型不穩定所導致，并給出了一種PD控制器來保持神經群的穩定性，抑制癲癇發作[11]。Ma等[12]針對同步行為的可傳遞性，提出了面向神經群族的控制方案，可以有效地通過打破神經群族的同步性來防止癲癇發作。
在本文中，我們提出了一種致癇興奮性控制方法。首先，給出癇性指數的定義，并用來描述癲癇發作的程度。依據腦電信號的癇性指數來產生控制信號，從而對癲癇興奮性進行控制，阻止癲癇發作。
1 基于神經生理學的神經群數學模型
神經群模型的結構如圖 1所示。在該模型中，錐體神經元和中間神經元通過突觸互相連接。圖中，實線表示的是大腦內部信號流，虛線表示外部控制信號流。
圖1
神經群模型
Figure1.
Neural mass model
本地神經群模型中，一個線性系統可以將興奮或者抑制興奮神經子群的平均行動電位發放率m(t)轉換為平均后突觸電位v(t)，也就是v(t)=he, i(t)*m(t)，可以表示為：
| $ {{h}_{\text{e, i}}}\left(t \right)=\frac{{{A}_{\text{e, i}}}}{{{\tau }_{\text{e, i}}}}t{{e}^{-\frac{1}{{{\tau }_{\text{e, i}}}}}}, $ |
式(1)中下標e和i分別表示興奮和抑制神經子群。Ae, i描述了興奮/抑制子群后突觸電位的特性，可以分別用來調節興奮和抑制突觸的強度，τe, i為時間常數。本文中，Ae, i由神經群本身的興奮/抑制程度和兩個外部控制器(EPIDC和IPIDC)的影響共同決定，也就是
| $ {{A}_{\text{e, i}}}=A_{\text{e, i}}^{n}-A_{\text{e, i}}^{c}, $ |
其中，Ae, in為神經群自身的興奮/抑制強度，Ae, ic為來自于控制器EPIDC或IPIDC控制信號。
非線性環節S(v)將一個神經元子群的平均膜電位變換為動作電位的平均發放率，這個變換用如下形式的Sigmoid函數來描述，
| $ S\left(v \right)=2{{e}_{0}}/\left(1+{{e}^{r\left({{v}_{0}}-v \right)}} \right), $ |
這里2e0是最大發放率，v0是對應于發放率e0的后突觸電位，r是Sigmoid函數的陡峭度。外部興奮輸入P(t)代表了附近群的影響，在本文中通過正均值的高斯白噪聲來表示，其均值為90，方差為30。錐體神經元和中間神經元之間興奮和抑制突觸的數目是由C1~C4四個連接參數來描述[4-5, 12]。
2 癲癇控制
2.1 特征提取
Pincus[13]在研究心率變化時提出了近似熵(approximate entropy，ApEn)的算法，越復雜的時間序列對應的近似熵越大。而當癲癇發作時，腦電信號的復雜度會降低，因此，本文采用Olson等[14]的方法來計算近似熵，并作為描述腦電信號復雜性的非線性指標。
此外，我們采用小波變換技術對腦電信號進行多尺度分解，提取對應頻段在3~29 Hz范圍內各尺度上重構信號的波動指數[15]作為腦電數據的線性特征。通常情況下，癲癇發作時腦電信號的波動會比未發作的間歇期劇烈，波動指數可以用來衡量信號變化的強度，其定義如下:
| $ {{F}_{i}}=\sum\limits_{j=1}^{M-1}{\left| {{a}_{i}}\left(j+1 \right)-{{a}_{i}}\left(j \right) \right|}, $ |
式中: ai為第i尺度上小波重構信號的幅度，M為信號的長度。
我們綜合二者，提出癇性指數(epileptiform index, Ie)來描述腦電信號的發作程度，定義如下：
| $ {{I}_{c}}=\sum\limits_{i\in N}{{{F}_{i}}/{{E}_{a}}}, $ |
這里，Fi為小波變換后第i尺度的的波動指數，N為納入考慮的尺度集合，在本文中，選取對應頻段在3~29 Hz范圍內的第3、4、5尺度，即N={3, 4, 5}，并分別計算其波動指數。Ea為癲癇信號的近似熵。
2.2 控制方法
我們所提出的控制策略的總體目標是通過降低興奮強度或者增加抑制強度來使得興奮和抑制處于平衡狀態，從而抑制發作。所以，我們引入兩個PID控制器EPIDC和IPIDC，其中EPIDC用來控制興奮線性模塊，產生控制信號Aec；IPIDC用來控制抑制線性模塊，產生控制信號Aic。兩個控制器的控制信號是癇性指數的函數。
首先，計算腦電信號的近似熵Ea，并進行小波變換，計算各尺度的波動指數Fi；然后，計算癇性指數Ie，如果大于閾值TI則認為興奮性過高，足以引發癲癇，控制器開始工作；最后，根據計算癇性指數Ie確定控制量，控制量初始置為0，如果需要進行調節，則更新控制量，否則保持控制量不變。控制量為：
| $ A_{\text{e, i}}^{c}\left(t \right)={{k}_{\text{p}}}*e\left(t \right)+{{k}_{\text{i}}}*\int{e\left(t \right)+{{k}_{\text{d}}}}*\frac{\text{d}e\left(t \right)}{\text{d}t}, $ |
其中，e(t)=Ie-TI。對于EPIDC，kp、ki、kd三個參數為正數，而對于IPIDC，這三個參數為負數。
3 實驗結果
構建包括一個神經群及兩個PID控制器的實驗系統，調節神經群的參數保持在標準值[5]，控制器保持在待命狀態以模擬正常的腦電行為。然后，通過改變Aen使得神經群的興奮性產生變化，來驗證本文所提出的控制器是否能有效地控制致癇興奮性。
3.1 無控制器時的發作情況
為了研究興奮性的增加對癇性指數的影響，保持Ain不變，Aen從2增加到11，步長為0.01，對應的癇性指數變化曲線如圖 2所示。從圖 2可以看到，在沒有控制器介入的情況下，癇性指數會隨Aen的增加而增大。當Aen=7.71時，癇性指數達到最大值34.350 7。當Aen在2~3.45范圍內變化，癇性指數一直保持增加，從0.684 5增加到1.350 3，但增加幅度不大。而當Aen=3.46時，癇性指數突然增加到4.225 1。隨著Aen繼續增加，癇性指數一直保持在較高值。當在Aen=7.71達到最大值后，癇性指數有一定程度下降，但是依然保持在25左右。
圖2
無控制器時癇性指數隨興奮性增加的變化
Figure2.
Epileptiform index as a function of excitatory strength in the case of no controller
3.2 有控制器時的發作情況
為了檢驗本文提出的控制器是否能較好地調節過度興奮引發的癲癇發作，通過增加興奮強度Aen來仿真神經群興奮性增加的情況，然后應用本文的控制器來對興奮性的增加進行控制。具體為Ain保持標準值，Aen從2增加到8，然后由8減少到2，步長為0.1。仿真時間為122 s，每秒鐘計算一次癇性指數Ie，如果超過閾值TI則控制器介入進行調節，TI設置為1.5，這里閾值取不同的值會使得控制器有不同的靈敏性，介入的時刻也不同。在仿真甚至實際應用中，TI可以根據患者的情況進行設置，使控制器可以更好地適應不同的病情。為了分別討論兩種控制策略(降低興奮強度和增加抑制強度)對癲癇發作的調節效果，EPIDC和IPIDC分別應用于神經群模型。EPIDC采用的kp、ki、kd參數值分別為0.11、0.003 5和0.002 4，IPIDC的kp、ki、kd參數值分別為-1.82、-0.122和0.116。
只有EPIDC工作的時候，Aen和Ace如圖 3(a)所示。從圖中可以看出，通過PID控制器，Ace能根據Aen的變化而有效地進行改變，給出的控制反饋可以抵消過度興奮所造成的影響，使興奮強度保持在正常水平。從圖 3(b)可以發現，采用這種方法保持了興奮-抑制平衡，癇性指數均值為1.961 1，與無控制器的情況相比有大幅度下降，保持在正常范圍，避免了發作。在A通過臨界點3.5的時候，無論是增加過程還是下降過程，都會引起癇性指數的突變，但控制器都可以很快地隨之產生相應的控制信號，從而有效地降低興奮強度。
圖3
只有EPIDC起作用時的控制結果
(a)
(a)
只有IPIDC工作的時候，Aen和Aic如圖 4(a)所示。在圖中，因為Aen和Aic數值差距較大，為了對比更加明顯，將Aen放大十倍進行顯示。從圖中可以看出，與EPIDC通過產生控制信號直接抵消興奮強度的增加不同，IPIDC根據癇性指數的變化產生控制信號來增加抑制強度，從而維持興奮-抑制平衡。從圖 4(b)可以看出，癇性指數并未像無控制器情況下隨興奮強度增加而迅速增加，均值為1.760 3。在仿真過程中，在興奮強度增加通過臨界點A=3.5的時候，會引發癇性指數的突然增加，控制器可以很快地產生增加抑制強度的控制信號，有效地控制了癲癇發作。這種策略雖然不通過直接降低興奮強度來維持興奮-抑制平衡，但是增加抑制強度也可以抵消興奮強度增加帶來的影響，保持二者之間的平衡。
圖4
只有IPIDC起作用時的控制結果
(a)
(a)
4 結論
當癲癇發作時，腦電信號從高復雜度過渡到低復雜度，本文用近似熵來描述腦電信號的復雜度，并作為表征癲癇發作程度的非線性指標。波動指數可以反映腦電信號的波動劇烈程度，而在癲癇發作時，腦電信號波動加劇，所以本文采用波動指數作為線性指標來描述癲癇發作程度。進而，提出可以反映腦電信號線性和非線性特征的癇性指數，并用來作為PID控制器的被控參數來維持大腦的興奮-抑制平衡。
采用神經群模型作為研究平臺，對興奮性增加引起的癇性指數的變化進行了討論。如果抑制強度保持不變，隨著興奮強度的增加癇性指數也會增加。當興奮強度超過某一臨界值，癇性指數會突然大幅度增加，導致癲癇發作。本文采用兩個PID控制器來模擬兩種控制興奮性增加的策略——降低興奮強度和增加抑制強度。仿真結果證明這兩種策略都可以很好地控制興奮性增加導致的癇性發作。本文所提出的控制方法為設計外部發作控制設備提供了理論的基礎。PID控制器結構簡單而且不依賴于潛在的生理結構，有利于將其應用于臨床。將此控制方法應用于動物模型將是我們下一步研究的重點。
