基于模型的帕金森病閉環(huán)深部腦刺激效果指標研究

趙德春 陳 歡 沈利豪 焦書洋 蔣宇皓

①(重慶郵電大學(xué)生物信息學(xué)院 重慶 400065)

②(重慶郵電大學(xué)自動化學(xué)院 重慶 400065)

1 引言

帕金森病(Parkinson’s Disease, PD)是一種常見的神經(jīng)退行性疾病，多發(fā)于中老年人。PD的臨床癥狀主要表現(xiàn)為靜止性震顫和運動遲緩等[1]。晚期PD患者生活無法自理，隨著我國老齡化程度的加重，這將給社會帶來沉重的壓力。有研究表明，基底節(jié)為PD的主要致病區(qū)域，但是其潛在的致病機制至今還沒有得到完美的詮釋。目前以基底節(jié)中丘腦底核(SubThalamic Nucleus, STN)為刺激靶點的深部腦刺激(Deep Brain Stimulation, DBS)已經(jīng)成為治療晚期PD的一種常見的外科手術(shù)方式。但存在一定的局限性，比如電池能量消耗過快導(dǎo)致更換電池時會增加患者的開銷與手術(shù)風(fēng)險，患者的癥狀每時每刻都在波動，長期刺激可能會對人體產(chǎn)生如構(gòu)音障礙、言語障礙等副作用[2–4]。如果能隨時監(jiān)測患者癥狀的波動并及時予以刺激，這不僅能節(jié)省電池的壽命，還會在提高刺激效果的同時減少副作用，由此推動了閉環(huán)DBS(Closed-loop DBS,CDBS)的發(fā)展。然而對CDBS算法的探索是一個漫長的過程，需要不斷的研究試錯才能最終應(yīng)用于臨床。但是直接以患者為研究對象是不可行的，同時以非靈長類動物直接進行實驗不僅見效慢還會增加實驗成本。因此建立以基底節(jié)為對象的計算模型，探索其與真實患者之間的內(nèi)在聯(lián)系，以及研究基于模型的CDBS算法有重要意義。

越來越多的證據(jù)表明，PD患者中丘腦底核局部場電位(SubThalamic Nucleus Local Field Potential, STN LFP)的β振蕩(12～35 Hz)可以被持續(xù)檢測到，并且其水平與運動障礙密切相關(guān)，是構(gòu)建閉環(huán)DBS的理想反饋信號[5,6]。Peter Brown團隊通過將β頻段內(nèi)的局部場電位(Local Field Potential, LFP)作為反饋信號，以丘腦底核(Sub-Thalamic Nucleus, STN)為靶點進行單側(cè)或雙側(cè)刺激，成功獲取了不同患者的實驗數(shù)據(jù)[6–8]，證實了CDBS的有效性與可行性。相較于傳統(tǒng)高頻DBS,CDBS刺激效果提升了30%左右，能量消耗大約減少了50%，同時還減少了語言方面的副作用[7,9]。

Terman和Rubin等人[10,11]在2002年時建立了丘腦底核與蒼白球外側(cè)(SubThalamic Nucleus-Globus Pallidus external, STN-GPe)耦合模型以探索神經(jīng)元之間的動態(tài)相互作用，并在2004年時將模型擴展為皮層-基底節(jié)-丘腦環(huán)路，此后的計算模型大多以Rubin和Terman等人建立的模型為基礎(chǔ)進行改進，建立了更加復(fù)雜，更加符合PD患者腦神經(jīng)活動的模型。但是其內(nèi)部數(shù)學(xué)模型卻不盡相同，比如常用的Hodgkin-Huxley方程[4]和Izhikevich方程[12,13]?；谀Ｐ偷腃DBS算法已經(jīng)進行了廣泛的研究，研究者一般采用丘腦中繼準確度[14–16]、β頻帶功率[4,17,18]、能量消耗[4,14,16–19]等作為衡量CDBS的刺激效果，也有采用基底節(jié)中各神經(jīng)元之間的同步性[12,17]作為評估指標。能量消耗是需要達成的目標之一，但在致力達到這個目標之前，需要保證刺激的效果。在臨床上，不管是在藥物治療[5,19–21]，開環(huán)DBS還是為數(shù)不多的閉環(huán)DBS刺激研究中[6–9]，對治療PD患者運動障礙的改善程度主要是通過患者填寫的統(tǒng)一帕金森病評定量表(Unified Parkinson’s Disease Rating Scale, UPDRS)來進行評估。所以為了更直觀，更加貼近臨床評估效果閉環(huán)DBS的刺激效果，建立新的評估指標將UPDRS評分與基于模型的閉環(huán)DBS研究聯(lián)系起來是有必要的。

對健康的非人類靈長類的動物研究中發(fā)現(xiàn)β爆發(fā)的時間動力學(xué)傳遞了重要的生理信息，即β爆發(fā)需要爆發(fā)[22]。在PD患者中也證實了這一點，Tinkhauser等人[6]通過對CDBS刺激后患者的數(shù)據(jù)從時域上進行統(tǒng)計分析，提出在CDBS刺激下持續(xù)時間較長的β爆發(fā)與運動障礙的改善程度呈正相關(guān)，持續(xù)時間較短的β爆發(fā)與運動障礙呈負相關(guān)的結(jié)論。Sun等人[23]也通過對CDBS刺激后PD患者的LFP從時域上進行特征提取，得到PD患者未刺激狀態(tài)下的β頻段有更多的長爆發(fā)而DBS狀態(tài)下有更多的短爆發(fā)的結(jié)論。因此本文基于以基底節(jié)為對象的計算模型，從時域分析CDBS刺激下LFP的爆發(fā)特性，并與PD患者LFP數(shù)據(jù)的時域特性進行對比，研究其臨床相關(guān)性。同時對β爆發(fā)持續(xù)時間的分類進行研究，提出與臨床相關(guān)的指標，探索其作為CDBS刺激后用于估計運動障礙改善的可行性，為今后研究CDBS算法奠定基礎(chǔ)。

2 方法

2.1 皮層-基底節(jié)-丘腦計算模型

基底節(jié)與人體運動控制等功能密切相關(guān)，基底節(jié)的損傷則是運動障礙疾病的生理基礎(chǔ)，比如帕金森病?；坠?jié)接收來自皮層(Cortex, Ctx)的信息，并傳遞給丘腦(Thalamus, Th)。模型中的基底節(jié)(圖1黃色區(qū)域為基底節(jié))主要包括紋狀體(Striatum,Str)、STN、蒼白球外側(cè)(Globus Pallidus external,GPe)、蒼白球內(nèi)側(cè)(Globus Pallidus internal,GPi)。本文使用與Kumaravelu等人[13]相同的基于電導(dǎo)的皮層-基底節(jié)-丘腦模型，如圖1所示，每類核團均包含10個神經(jīng)元，每個核團之間的連接包含興奮性連接與抑制性連接，根據(jù)網(wǎng)絡(luò)之間的連接可以得到各神經(jīng)元的膜電位方程，具體內(nèi)容見文獻[13]。該模型的正常狀態(tài)，患病狀態(tài)與STN DBS狀態(tài)下基底節(jié)各神經(jīng)元的神經(jīng)振蕩活動與平均放電率已經(jīng)經(jīng)過6-羥基多巴胺(6-HydroxyDopamine, 6-OHDA)損毀的大鼠實驗進行了驗證。

圖1 皮層-基底節(jié)-丘腦網(wǎng)絡(luò)模型

構(gòu)建模型的主要目的是獲取STN核團處的局部場電位(Local Field Potential, LFP)以進行下一步的研究，LFP計算公式為[24,25]

其中，R= 1 ?·m 表 示均勻的細胞外電阻率，ISTNi為第i個STN神經(jīng)元的總輸入電流，Dci為 第i個STN神經(jīng)元與中心記錄電極之前的距離，假設(shè)每個神經(jīng)元與刺激電極之間的距離相等Dci= 1 m m。

2.2 實驗設(shè)置

在PD狀態(tài)時，由于多巴胺能神經(jīng)元的缺失，紋狀體中M型鉀離子電流減少，紋狀體對接收皮層輸入的敏感性降低，同時GPe內(nèi)軸突側(cè)支的突觸強度加強，所以改變對應(yīng)的參數(shù)即可模擬PD狀態(tài)[13]。為了避免模型的偶然性對后續(xù)研究造成影響，通過改變參數(shù)的初始值模擬10名PD“患者”。有證據(jù)表明STN LFP中β頻帶的過度振蕩與帕金森病的運動障礙密切相關(guān)[6]，所以在進行后續(xù)研究之前先驗證模型的可行性。如圖2為模擬的10名PD患者中Patient 2在LFP上β頻帶的功率，相較于正常狀態(tài)，患病狀態(tài)下的功率表現(xiàn)出了明顯的急劇上升，特別在峰值范圍內(nèi)(19±4 Hz)。這很好地模擬了β頻帶的劇烈振蕩，同時具有較好的魯棒性。

圖2 LFP中β帶功率譜密度

目前報道CDBS改善PD患者臨床癥狀的研究者主要有Peter Brown團隊，他們廣泛采用的閉環(huán)控制策略為閾值刺激法[6–9]，接下來為了使研究更具有實際應(yīng)用意義本文將以閾值刺激法為例進行后續(xù)研究。以往的研究表明，刺激閾值至少大于β振幅的50%時將具有良好的臨床效果[6]，并且具體數(shù)值與患者的狀況相關(guān)，因此本文隨機選取大于β振幅52%分位處的數(shù)值作為刺激閾值。獲取PD患者在未刺激狀態(tài)下的LFP，并圍繞β頻段內(nèi)的峰值活動頻率范圍進行帶通濾波，然后包絡(luò)以產(chǎn)生β振幅，通過計算獲得PD患者在分位處的閾值大小為在線刺激提供參考。基于閾值控制策略的閉環(huán)DBS刺激如圖3所示，首先獲取未刺激狀態(tài)下的LFP，然后在β頻帶的峰值頻率附近進行帶通濾波，接下來對β爆發(fā)進行包絡(luò)以獲得β振幅，結(jié)合在未刺激狀態(tài)下得到的刺激閾值對β振幅的大小進行篩選，當振幅大于閾值時通過DBS波形發(fā)生器產(chǎn)生刺激波形對STN進行刺激直到振幅小于閾值才停止刺激。在開始刺激之前延遲時間為30 ms，同時為了避免神經(jīng)組織受損，在真正的DBS外科手術(shù)中常采用電荷平衡雙相脈沖[26]。所以本文采用脈沖寬度為100 μs，先陽極后陰極的130 Hz電荷平衡雙相對稱脈沖[8]。

獲取每種狀態(tài)下200 s的數(shù)據(jù)，然后使用Matlab2017a進行處理。對數(shù)據(jù)進行分析時為了提高計算效率，將數(shù)據(jù)降采樣至200 Hz。在頻域上，理論上DBS刺激應(yīng)該能消除β頻段的劇烈振蕩，因此與圖2類似，通過對比刺激與未刺激這兩種狀態(tài)下β功率的大小來評估刺激效果。在時域上，主要分析β爆發(fā)的持續(xù)時間與平均幅值的時域特性。首先將β爆發(fā)定義為從β振幅大于刺激閾值開始一直持續(xù)到小于閾值的時間段，如圖3的紅色橢圓框所示，而平均幅值定義為這段持續(xù)時間所對應(yīng)的幅值的平均值。

圖3 閉環(huán)DBS刺激流程

在閉環(huán)DBS刺激期間存在頻繁的刺激打開與停止的情況，刺激后的信號存在更大的方差[6]。對時域數(shù)據(jù)進行分析時，在相同百分位的閾值下β爆發(fā)的持續(xù)時間長短受不同狀態(tài)下閾值大小的影響，所以為了避免低估刺激狀態(tài)下的持續(xù)時間將閾值設(shè)定為兩種狀態(tài)下閾值的平均值，并將該閾值應(yīng)用于所有狀態(tài)。同樣地，不同百分位的平均閾值也能影響爆發(fā)持續(xù)時間的長短，因此為了尋求不同閾值下數(shù)據(jù)分布的相似性，定義了不同百分位(55%, 60%,65%, 70%, 75%)的閾值。以65%分位下的閾值為例，圖4為Patient 2在兩種狀態(tài)下持續(xù)時間與幅值分布的散點圖，從中可以觀察出β爆發(fā)持續(xù)時間的大小主要集中在500 ms以內(nèi)。為了避免峰值頻段以外的噪聲影響，只保留大于100 ms的數(shù)據(jù)。但保留大于500 ms的數(shù)據(jù)也是有必要的，因為持續(xù)時間更長的突發(fā)變得越來越不頻繁，并且未刺激狀態(tài)下?lián)碛懈喔L的持續(xù)時間，其可能與PD癥狀的致病相關(guān)[6]。為了更加具體地研究不同爆發(fā)持續(xù)時間的分布情況，將β爆發(fā)劃分為9個時間窗口以保證每個時間窗內(nèi)都能擁有足夠的β爆發(fā)個數(shù)。同時根據(jù)分布情況進一步將β爆發(fā)的持續(xù)時間劃分為2個時間窗口以觀察長爆發(fā)與短爆發(fā)的時域分布情況。

圖4 65%分位閾值下β爆發(fā)時域圖(Patient 2)

由于每名模擬患者β爆發(fā)的總數(shù)上都存在差別，所有數(shù)據(jù)均以平均值±標準差的形式出現(xiàn)，使用單因素重復(fù)測量方差分析數(shù)據(jù)之間的顯著性，在此之前采用Kolmogorov Smirnov測試對數(shù)據(jù)進行正態(tài)分布檢測，不滿足正態(tài)分布的采用Wilcoxon秩和檢驗。

3 結(jié)果分析

不管是外科手術(shù)治療[3]還是藥物治療[5]都會降低PD患者的β頻帶功率，因此在閉環(huán)DBS的算法研究中廣泛采用β帶功率的變化作為刺激效果的評估指標。而本文基于模型的閾值刺激策略也能顯著降低β頻段的功率(F=8.77, p<0.01) (圖5)，所以本文應(yīng)用的閉環(huán)刺激算法是可行的。

圖5 不同狀態(tài)下β頻帶的功率變化

圖6為200 s的時間長度下兩種狀態(tài)β爆發(fā)的分布情況，圖6(a)表明刺激狀態(tài)下位于0.10～0.15 s的時間窗內(nèi)β爆發(fā)數(shù)量顯著多于未刺激狀態(tài)(F=11.32,p<0.01)，相反的是在0.20～0.25 s, 0.30～0.35 s和大于0.5 s的時間窗內(nèi)未刺激狀態(tài)下β爆發(fā)的數(shù)量數(shù)明顯多于刺激狀態(tài)下的數(shù)量(F=7.14, p<0.05;F=7.27, p<0.05; F=5.28, p<0.05)，其中?表示p<0.05，??表 示p<0.01，???表示p<0.001。結(jié)果表明與未刺激狀態(tài)相比，DBS狀態(tài)下短爆發(fā)(0.1～0.15 s)的數(shù)量更多，而長爆發(fā)(0.20～0.25 s,0.30～0.35 s和大于0.5 s)的數(shù)量更少。為了分析兩種狀態(tài)下爆發(fā)分布的具體差異，同時探尋不同閾值下β爆發(fā)分布的一致性，在相同條件下本文分析了所有閾值下β爆發(fā)的分布情況并總結(jié)得到了不同時間窗內(nèi)β爆發(fā)的差值分布(圖6(b))。這里?表示所有閾值下存在顯著性差異的時間窗，由于剛好在0.20～0.25 s的時間窗內(nèi)刺激狀態(tài)下的數(shù)量顯著多于未刺激狀態(tài)，所以本文進一步將時間窗二分類為0.10～0.25 s的短爆發(fā)與大于0.25 s的長爆發(fā)。

圖6 β爆發(fā)的分布情況

如圖7(a)，兩種狀態(tài)在二分類下的爆發(fā)分布具有顯著性差異(F=11.32, p<0.01)，與未刺激狀態(tài)相比DBS刺激下具有更多的短爆發(fā)和更少的長爆發(fā)，但是在總的爆發(fā)數(shù)量上兩者卻沒有顯著性差異(F=0.21, p>0.05) (圖7(b))。同樣在平均持續(xù)時間上也能得到相似的分布(圖7(c)和圖7(d))，即在總的持續(xù)時間沒有顯著性差異的情況下(F=4.4,p>0.05)，刺激狀態(tài)下長爆發(fā)對應(yīng)的總持續(xù)時間顯著大于未刺激狀態(tài)，而短爆發(fā)下則顯著少于未刺激狀態(tài)(F=10.58, p<0.01; F=8.74, p<0.01)，并且所有不同閾值下得到的結(jié)論一致。根據(jù)以往的研究經(jīng)驗，短爆發(fā)可能是人體產(chǎn)生的正常爆發(fā)，而長爆發(fā)可能是所謂的“致病”爆發(fā)[27]，PD癥狀的產(chǎn)生可能與此相關(guān)。結(jié)論表明CDBS能把未刺激狀態(tài)下較長的爆發(fā)轉(zhuǎn)移為更多較短的爆發(fā)，所以刺激狀態(tài)下?lián)碛懈嗟亩瘫l(fā)和更大的總持續(xù)時間，這有可能是改善PD癥狀的機制之一。

圖7 65%分位下的β爆發(fā)分布

結(jié)合以往的研究經(jīng)驗，開環(huán)DBS和閉環(huán)DBS都能在一定程度上降低β帶功率，雖然開環(huán)DBS降低的功率更多，但閉環(huán)DBS的改善效果更好[7]，所以在評估閉環(huán)DBS的刺激效果時單一地使用β功率或許不夠全面。Tinkhauser等人[6]從時域分析得出β爆發(fā)的持續(xù)時間與UPRDS評分之間存在相關(guān)性，即持續(xù)時間較短爆發(fā)與運動障礙呈負相關(guān)，而持續(xù)時間較長的爆發(fā)與運動障礙呈正相關(guān)。也就是說較長持續(xù)時間的β爆發(fā)所占的比例越多，運動障礙相關(guān)的癥狀越嚴重，因此本文旨在研究相關(guān)的指標用于估計UPDRS評分。

根據(jù)UPDRS經(jīng)過DBS刺激后評分可能會降低的特性，本文將估計指標定義為長爆發(fā)與短爆發(fā)之間的比值，并命名為UPDRS類估計(Similar to UPDRS Estimates, SUE)，即

其中，Ldur為持續(xù)時間較長的β爆發(fā)的數(shù)量，Sdur為持續(xù)時間較短的β爆發(fā)數(shù)量，這里長爆發(fā)與短爆發(fā)見圖7(a)。

為了確認未刺激時SUE與持續(xù)時間之間是否相關(guān)，同時確認不同閾值之間的一致性，本文根據(jù)式(2)對所有閾值下的不同時間窗與SUE進行相關(guān)性分析。圖8(a)表明在所有的閾值下SUE與持續(xù)時間較短的爆發(fā)呈負相關(guān)，而與持續(xù)時間較長的爆發(fā)呈正相關(guān)，示例如圖8(b)和圖8(c)所示。這和運動障礙與爆發(fā)持續(xù)時間之間的臨床相關(guān)性所得的結(jié)論一致，因此初步認定SUE可用于估計運動障礙的程度(UPDRS-III第20, 22, 23項)。

圖8 持續(xù)時間與SUE之間的相關(guān)性

如圖9所示，經(jīng)過DBS刺激后SUE值顯著低于未刺激狀態(tài)下(F=18.42, p<0.001)，并且其改善效果較好，平均改善率幾乎都能達到30%以上，這與實際情況相符，即經(jīng)過有效的DBS刺激之后PD患者的UPDRS-III評分降低。然而有趣的是，從SUE的結(jié)果上來看模擬的10名患者中并不是所有的患者都對DBS的刺激有積極的響應(yīng)，其中Patient 9的癥狀不僅沒有得到改善，反而加重。產(chǎn)生該現(xiàn)象的原因可能是由于患者之間的差異性，如文獻[6]描述的那樣，每名患者達到較好的刺激效果時，所采用的刺激閾值是不同的。說明不適當?shù)拇碳ら撝禃又鼗颊叩陌Y狀，從神經(jīng)元的振蕩方面來解釋就是會增加長β爆發(fā)的產(chǎn)生。因此SUE不僅能很好地反映DBS刺激效果較好時UPDRS-III的變化，同時還能反映效果不好時的變化，表征患者之間的差異性。

圖9 不同狀態(tài)下的SUE

4 結(jié)束語

本文以皮層-基底節(jié)-丘腦網(wǎng)絡(luò)模型為平臺，閾值刺激為閉環(huán)控制策略研究能夠估計臨床UPDRS評分的指標。首先通過頻域內(nèi)LFP β頻帶的功率變化驗證了計算模型與閉環(huán)控制算法的可行性，然后在此基礎(chǔ)上對比分析了LFP的時域特性，并提出了能估計在不同情形下動態(tài)變化的UPDRS-III指標。并且計算模型得到的β爆發(fā)持續(xù)時間更短，因此每個時間窗口的長度取較小值(50 ms)以保證充足的數(shù)據(jù)量。為了避免數(shù)據(jù)長度對結(jié)果產(chǎn)生的影響，分析了25 s和100 s的數(shù)據(jù)長度，得到了與200 s數(shù)據(jù)長度下相似的結(jié)論，閾值范圍為55%～75%，并且在不同百分位閾值下數(shù)據(jù)分布一致性較高，能顯著性地區(qū)分長爆發(fā)與短爆發(fā)。

在以往的研究中對計算模型有效性的驗證往往是根據(jù)信號的放電率，神經(jīng)元之間的同步性，或者丘腦中繼準確度等指標。Sun等人[23]通過特征提取的方式驗證了真實PD患者的刺激狀態(tài)比未刺激狀態(tài)擁有更多的短爆發(fā)，Tinkhauser等人[6]也經(jīng)過統(tǒng)計分析得到了相似結(jié)論。本文通過在計算模型上進行的時域分析得到了相同的結(jié)論，該結(jié)論在頻域的基礎(chǔ)上更加符合DBS刺激狀態(tài)下的真實特性，因此今后在對閉環(huán)DBS研究之前可以通過持續(xù)時間的顯著性分布來驗證計算模型與閉環(huán)刺激策略的可行性。

根據(jù)Peter Brown團隊的研究結(jié)論，LFP β爆發(fā)的持續(xù)時間與患者的臨床損害相關(guān)。該現(xiàn)象直接將人體信號的內(nèi)部特性與外部表現(xiàn)聯(lián)系了起來，而本文所提出的SUE指標能在一定程度上描述UPDRS-III在不同狀態(tài)下的變化，將模型研究與臨床效果相關(guān)，使得今后基于模型的閉環(huán)DBS算法研究更具有可靠性。此外通過時域特性分析表明閾值刺激法能通過抑制其可能“致病”的長爆發(fā)，誘導(dǎo)短爆發(fā)的產(chǎn)生來達到改善PD癥狀的目的。這或許是閉環(huán)DBS的治療機制之一，因此基于模型的時域分析法也為今后探索閉環(huán)DBS的潛在治療機制奠定基礎(chǔ)。