創傷后應激障礙（PTSD）臨床表現復雜多樣，僅依靠臨床評估進行診斷缺乏準確性和客觀性。因此，迫切需要建立可靠客觀的輔助診斷模型，為PTSD患者提供有效的診斷。目前，利用圖神經網絡表征PTSD受限于已有模型的表達能力，并不能取得較好的分類診斷結果。為此，本文提出多圖譜多核圖卷積模型對PTSD數據進行分類。首先，采用不同圖譜構建同一受試者不同尺度下的功能連接矩陣，并通過K最近鄰算法構建為圖。其次，對同一受試者不同尺度的腦圖結構，提出多圖譜多核圖卷積方法來提高模型特征的提取能力。最后，將提取到的不同尺度下多種受試者特征融合后進行分類，并使用圖類激活映射確定Top 10的分類貢獻腦區。在地震后PTSD數據上的實驗結果表明，本文模型在地震后PTSD患者和未患病受試者分類任務上獲得了84.75%的準確率、84.02%的特異性及85%的曲線下面積。本文結果為地震后PTSD輔助診斷提供了有力證據，并有望在其他PTSD疾病分類診斷中也能可靠地識別特異性的大腦區域，為臨床醫生提供參考。
引用本文： 周麗君, 朱鴻儒, 劉云飛, 莫賢, 袁軍, 羅昌宇, 張俊然. 基于多圖譜多核圖卷積的創傷后應激障礙分類研究. 生物醫學工程學雜志, 2024, 41(6): 1110-1118. doi: 10.7507/1001-5515.202407031 復制
版權信息： ?四川大學華西醫院華西期刊社《生物醫學工程學雜志》版權所有，未經授權不得轉載、改編
0 引言
創傷后應激障礙（post-traumatic stress disorder，PTSD）是一種在創傷性事件之后較常見的精神障礙疾病，常見的創傷源為戰爭、自然災害（如地震）、車禍等。其特征是對創傷性事件的持續侵入性記憶、過度覺醒、回避與創傷相關的線索，以及思維和情緒的負面變化[1-2]。這種慢性精神疾病嚴重影響著患者的生活質量，并可造成相當大的經濟負擔[3]。由于PTSD的臨床表現復雜多樣，目前仍缺乏可靠的客觀生物學標志，單純依靠量表和醫生的主觀經驗來進行臨床評估和診斷可能造成準確性的偏差和一致性缺乏[4]，同時，由于客觀診斷手段（如影像標記）的可重復性差，亦會導致大量PTSD患者未被發現或誤診[5]。因此，迫切需要建立可靠、客觀、指標可重復的輔助診斷模型，通過客觀手段重現PTSD的臨床測量來彌合主觀評估和客觀測量之間的差距[6]，并促進未來診斷方法的改進。
目前，研究者主要基于功能磁共振成像（functional magnetic resonance imaging，fMRI）數據來建立大腦的功能連接網絡（functional connectivity network，FCN）等特征映射[7]。隨后，利用機器學習或深度學習方法，學習和選擇精神障礙疾病特征映射中相較于傳統方法（如稀疏學習等）更敏感、更有效的特異性指標，以提高精神障礙疾病診斷的準確性和可靠性[8-10]，并將參與分類的主要貢獻特征作為與疾病高度相關的影像標記。Campbell等[11]結合多個尺度特征（即區域均勻性、低頻波動的局部幅度和基于功能連接的網絡），訓練多層感知機對有意識清醒和麻醉誘導的無意識狀態進行分類，用于從靜息態fMRI數據中識別麻醉誘導的無意識狀態的穩健影像標記，取得了0.84的平均曲線下面積（area under the curve，AUC）。Chang等[12]使用深度堆疊自編碼器提取低頻振幅的低維特征識別主要精神疾病及其亞型，并推定了新的診斷影像標記，為精神疾病的精準診斷提供了啟發。值得注意的是，雖然已有大量研究將影像數據與深度學習方法結合起來分類精神分裂癥[13]、重度抑郁癥[14]、自閉癥[15]、阿爾茨海默病[16]等精神疾病的不同群體，但利用深度學習方法對PTSD進行分類診斷的研究還很有限。Zhu等[17]將深度學習技術結合基于圖論的分析指標識別PTSD患者獲得了80%的平均準確率，高于傳統支持向量機（57%）。Zhu等[18]使用自編碼器用于PTSD多站點數據的分類，取得了75%的AUC值。盡管上述方法在捕獲和提取與PTSD相關的大腦功能連接特征方面取得了一定的成果，能夠從高維fMRI數據中提取復雜的空間和時間特征，減少了對人工提取特征的依賴[6]，但是，目前使用深度學習等方法對PTSD的分類準確率還有待提高[17]，同時，這類方法對精神障礙疾病類別歸屬的可解釋性也有限[18]，這可能降低臨床醫生的信任度和接受度。
最近，圖神經網絡（graph neural network，GNN）已經成為一種強大的深度學習工具，被廣泛用于解決醫學fMRI圖像中非歐氏空間圖數據問題[19]。GNN還能夠將重要特征可視化，以反駁深度學習是“黑匣子”的典型批判，能夠用于揭示精神疾病的影像特征[20]。GNN通過將大腦的功能連接網絡表示為圖來建模大腦的復雜網絡拓撲結構，其中，節點對應腦模板定義的感興趣區域（regions of interest，ROI），邊量化了ROI之間的功能連接。GNN已在之前的研究中已成功應用于自閉癥譜系障礙[21]、阿爾茨海默病[22]、重度抑郁癥[23]等精神障礙疾病的表征，取得了顯著的成果，各項分類指標提高的同時增強了對精神疾病的影像標記的理解。然而，之前的研究發現，使用單一圖譜劃分ROI，忽略了不同空間尺度下FCN潛在的互補拓撲信息[24]，單一圖譜并不能完全捕捉大腦區域連接的復雜性和異質性[25-26]，這種限制可能會削弱大腦網絡的多尺度特征刻畫能力。此外，傳統圖卷積網絡（graph convolutional network，GCN）模型[27]在所有層中使用恒定大小的卷積核，導致了每個節點的特征以固定的跳數學習鄰居的特征，不利于獲取類可分離的輸出特征[28]，從而限制了模型的特征提取能力，導致分類性能受到影響。
鑒于上述問題，本文提出多圖譜多核圖卷積（multi-kernel graph convolution neural network，MK-GCN）模型用于PTSD的分類診斷。主要在兩個方面做出了改進：首先，本文采用多個腦圖譜劃分ROI，以構建多尺度的FCNs進行表征學習，從輸入層面獲取多個圖譜FCN之間的互補拓撲信息；其次，通過改變GCN層中卷積核的大小來捕獲FCN中類可分離的輸出特征。最后，本文將從臨床角度出發，分析PTSD患者和創傷后未患PTSD的健康受試者之間的分類貢獻腦區，并在不同腦圖譜的模板下進行重疊腦區的確認，期望更可靠地識別地震后PTSD患者異常的大腦區域，進而為臨床醫生提供參考。
1 材料和方法
本文方法主要分為多尺度圖構建、多核圖卷積特征提取、分類輸出等幾個部分，模型整體框架如圖1所示：主要包括不同圖譜功能連接矩陣構建、K最近鄰算法選擇邊構建圖、多核圖卷積提取特征、特征融合后全連接分類輸出。首先，本文采用先驗圖譜自動解剖標記（automated anatomical labeling，AAL116）和基于數據驅動的腦圖譜（craddock200，CC200）分別對預處理后的信號構建功能連接矩陣。然后，使用K最近鄰算法選擇邊后構建為圖數據。接著，將不同尺度的圖數據送入多核圖卷積網絡中進行模型訓練和特征提取，捕獲更優特征。最后，將每個尺度提取到的特征融合后輸出到全連接層完成二分類，實現地震后PTSD組和地震后未患PTSD（Non-PTSD）的健康對照組的識別。
圖1
地震后應激障礙預測整體框圖
Figure1.
Post-earthquake traumatic stress disorder classification and prediction framework
1.1 數據采集及預處理
1.1.1 受試者信息
本研究招募的138名參與者均為2008年5月12日汶川地震后的幸存者，最后納入2015年至2017年來自同一地震災區的受試，年齡18～60歲，均為右利手，未接受任何常規治療（藥物或心理治療）。研究中兩名精神科醫生根據《精神疾病診斷和統計手冊》中創傷后應激障礙量表（Clinically Administered PTSD Scale，CAPS）評估PTSD癥狀，核磁共振掃描和心理測試在入組后兩天內進行評估。受試者排除標準包括：有任何神經系統疾病史、嚴重頭部損傷且意識喪失時間超過10 min、雙相情感障礙或其他精神障礙病史的參與者；智力低下、懷孕、酗酒和/或其他物質依賴/濫用；起搏器或手術夾等金屬植入物的干擾，以及其他磁共振成像掃描禁忌等情況。本研究經四川大學華西醫院醫學倫理委員會批準，并獲得每位受試者的書面知情同意書。對采集后的圖像質量進行了篩選，最后保留了110例受試者的靜息態fMRI數據，受試者信息如表1所示。
表1
110名受試者的人口統計學信息
Table1.
Demographic information of the 110 subjects
最終受試者分為兩組：① 59名PTSD患者，其中男性16名，女性43名。② 51名地震后未患PTSD的健康受試者，其中男性16名，女性35名，作為對照組。
1.1.2 靜息態fMRI數據采集
本研究110名受試者的fMRI數據采集均在3.0T MRI掃描儀上使用12通道相控陣頭部線圈進行，使用梯度回波平面成像序列獲取fMRI血氧水平依賴性圖像，其參數設置如下：重復時間/回波時間=2 000 ms/30 ms；翻轉角=90°；切片厚度=5 mm；視野=240 mm×240 mm；體素大小=3.75 mm×3.75 mm×5 mm；采集矩陣=64×64。所有受試者的掃描時長均接近7 min，對應205個時間點。
1.1.3 數據預處理
使用神經影像腦連接組分析軟件（graph theoretical network analysis，GRETNA）對采集數據進行預處理。處理步驟包括：① 考慮到初始靜息態fMRI信號不穩定和參與者前期的適應時間，去除每個功能數據的前10個時間點，剩余195個時間點用于后續處理。② 剩余的數據首先進行時間層校正，然后重新對齊以去除頭動偽影的影響。將頭動位移超過2.5 mm或者頭動旋轉角度超過2.5°的受試者剔除。本研究中所有受試者的頭動位移和頭動旋轉角度均未超過上述標準。③ 回歸協變量后，使用蒙特利爾神經學研究所（Montreal Neurological Institute，MNI）的模板將得到的圖像配準至3 mm體素的標準空間。④ 最后對每個體素時間序列進行帶通濾波（0.01～0.08 Hz），以減少低頻漂移和高頻生理噪聲的影響。
1.2 模型描述
同一受試者預處理后的圖像使用AAL116和CC200進行不同尺度腦區劃分；同時分別建立功能連接矩陣后對兩個尺度的功能連接（functional connectivity，FC）矩陣構建圖并進行多核卷積圖表征學習，最后在特征融合后進行分類輸出。
1.2.1 多尺度圖構建
本文使用三元組
表示具有N個節點的無向圖。其中
是N個腦區的集合，邊的集合
表示FC矩陣通過K最近鄰算法稀疏化后保留的連邊。
表示節點特征矩陣，其中F表示節點特征維度，矩陣中的每個元素表示腦區i和j之間的皮爾遜相關系數值。通過上述操作，對于同一受試者得到兩個特證矩陣：
和
，同時獲得兩個鄰接矩陣
和
，最后，基于兩個空間尺度為每個受試構造兩個圖
和
。
1.2.2 多核卷積圖表征學習
使用譜圖卷積定義整個卷積過程來分析大腦的FCN，并提出多核卷積的圖表征方法診斷PTSD。圖上的傅立葉變換依賴于拉普拉斯矩陣的特征向量，歸一化后的圖拉普拉斯矩陣L可以表示為式(1)。
 |
其中，
表示單位矩陣，
是度矩陣，在這里表示腦區連接的數量，
表示由傅里葉基組成的特征矩陣，
表示特征值的對角矩陣。則對于每個腦區的信號x的傅里葉變換就可以表示為
，其逆變換就可以表示為
。譜圖卷積定義為信號x與傅里葉域中可學習的濾波器
，為了防止信號x的傅立葉變換矩陣乘法所需的巨大計算量，使用切比雪夫多項式參數化重新定義濾波器表示為式(2)：
 |
其中，
是k次切比雪夫系數的向量，同時，定義
為一個函數
，對信號進行頻譜濾波，以腦區i為中心的濾波器
的腦區j的值可以表示為式(3)：
 |
受文獻[28]啟發，當k大于0且為整數時，代表了k跳鄰居，則對于任意的腦區節點i和j則有式(4)成立。
 |
其中，
是
和
之間的最短路徑距離，
是從
到
的最短路徑上的所有邊權重的總和。式（3）中拉普拉斯算子的k階多項式表示的譜濾波器恰好反映其k跳鄰域局部化表示。對于完整特征集X，使用濾波器
的輸出
可以表示為式(5)。
 |
為了全面考慮腦區特征之間的緊密程度及其鄰域的情況，本文使用多個不同大小的濾波器
（多核卷積）獲取對應的輸出特征
，其中，
表示單跳鄰居（最小）和6跳鄰居（最大）的內核大小，每個內核卷積后的輸出特征聚合后輸出至下一層。同時為了保證節點特征不被平滑化，選用3層GCN模型提取特征如式(6)所示。
 |
其中，
表示第l層神經網絡的權重矩陣，
，
，
表示歸一化鄰接矩陣和單位矩陣的和，
是非線性激活函數。輸出特征
，d表示最終輸出的特征維度。
最后，采用Lee等[29]提出的全局平均池化操作聚合節點特征并讀出同一受試的最終圖級表征，可表示為式(7)：
 |
其中，
表示卷積運算得到的第i個腦區的特征向量。整個多核卷積和圖讀出過程用圖2進行描述。
圖2
多核圖卷積網絡詳細描述
Figure2.
Detailed description of multi-kernel graph convolutional networks
1.2.3 特征融合與分類輸出
通過卷積操作，特征提取后最終獲取了每個尺度的128維的圖表示向量，分別表示為
和
，M表示受試者個數。如圖1所示，本文使用拼接操作進行特征融合，同時采用自適應權重策略自動獲取每個尺度的特征融合比例。最后，通過特征融合得到每個受試的兩個尺度特征更新后的256維圖表示向量，將此特征進行分類得出預測結果，并用圖類激活映射找出分類貢獻腦區（Top 10）。
1.3 實驗設置
本文模型基于PyTorch和PyG實現。3個圖卷積層的神經元個數分別為128、64、32，使用ReLU非線性激活和dropout操作防止模型過擬合，dropout設置為0.3。采用余弦退火學習率動態調整學習率，其初始值為0.001。同時采用帶標簽平滑的交叉熵損失作為優化函數，模型參數使用網格搜索進行優化，訓練迭代次數為100次。
實驗中進行了5折交叉驗證。評價指標包括分類精度、靈敏度、特異性。同時，繪制了受試者工作特性（receiver operating characteristic，ROC）曲線并計算了AUC，曲線的橫坐標為假陽性率，縱坐標為真陽性率。ROC曲線的AUC可以定量評價PTSD組和Non-PTSD組分類的準確性。本文采用雙樣本t檢驗進行統計學分析，所有統計分析的檢驗水準設定為0.05。顯著性檢驗結果的解釋以P值為準，P值小于0.05被認為具有統計學意義。數據分析使用Python 3.8.10進行。
1.4 對比方法
本文提出的MK-GCN方法與傳統機器學習方法——支持向量機（support vector machine，SVM）和邏輯回歸（logistic regression，LR），以及經典GCN模型[27]進行了對比實驗，對比模型具體實驗設置及參數參見附件1。
2 實驗結果
2.1 分類結果
在PTSD數據上進行對比實驗，對比結果如表2所示。同時，圖3直觀地顯示了這幾種方法的ROC曲線。從表2和圖3可看出以下結果：首先，在四個評價指標和ROC曲線上，本文提出的MK-GCN_F方法在PTSD和Non-PTSD分類中都取得了最優結果；其次，與單尺度結果相比，MK-GCN_F在四項評價指標上提高了至少2%；此外，對于單一尺度，使用CC200腦圖譜劃分ROI比AAL116腦圖譜劃分ROI取得了更好的結果，同時，使用GCN的方法比機器學習方法獲得了更優的分類性能，而多核圖卷積MK-GCN獲得了比傳統GCN方法至少高4%的AUC。這些結果表明，與使用單一腦圖譜劃分ROI相比，使用多尺度腦圖譜有助于提高PTSD分類性能，其原因可能是不同尺度FCNs特征包含互補的拓撲信息，這些互補信息的協同可以改善分類結果；此外，通過改變圖卷積過程濾波器的大小，可以改善GCN的特征提取能力，更好地捕捉類可分離的輸出特征，從而進一步提高PTSD分類性能。為了進一步驗證本文模型的的泛化能力，在公開的自閉癥譜系障礙（autism spectrum disorder，ASD）數據集ABIDE上取了最大站點NYU進行了驗證，此站點共有受試184人，ASD 79人，健康受試者105人；其中ASD患者女性11人、男性68人，健康受試者女性26人、男性79人，年齡、性別、教育時間均無顯著差異。實驗結果表明，本文模型在ASD中的分類結果相較于對比模型各項指標也有明顯提升，分類精度、靈敏度、特異性至少提升了4%，取得了86.95%的AUC值，完整數據參見附件2。
表2
不同方法在PTSD數據中評價指標對比結果（%）
Table2.
Comparison of evaluation metrics across different methods for PTSD data (%)
圖3
PTSD分類受試者工作特征曲線
Figure3.
Receiver operating characteristic curve for PTSD classification
2.2 不同腦圖譜分類結果
本研究還采用Power264（基于元分析將皮質區域映射到整個皮層的更精細的功能劃分圖譜）和Zalesky980（基于全腦解剖連接建模的更精細的圖譜）來觀察不同腦圖譜對分類結果的影響。實驗結果發現，與AAL116和CC200腦圖譜相比，更精細的腦圖譜對PTSD分類性能反而有所下降。為了更全面地評估多圖譜融合對分類結果的影響，本文還研究了圖譜兩兩組合、3個圖譜組合、4個圖譜組合對PTSD分類結果的影響，實驗結果見表3。
表3
不同腦圖譜及融合后對比結果（%）
Table3.
Different brain atlases and post-fusion comparison results (%)
2.3 不同大小卷積核分類結果
為了研究不同卷積核
的大小對每個尺度輸入圖性能的影響，設定在相同的輸入（即相同的節點特征和圖結構表示）下，ks1和ks2
[1, 6]的分類性能，其中，ks1表示第二層GCN的內核大小，ks2表示第三層GCN的內核大小，ks = 1和ks = 6分別表示了單跳鄰居（最小）和6跳鄰居（最大）的內核大小。對比結果見圖4（以CC200為例）。
圖4
不同層卷積核ks內核大小對分類結果影響熱圖
Figure4.
Heatmap of the effect of different layers of convolutional kernel kernel size of ks on classification results
2.4 分類貢獻Top 10腦區及可視化結果
本文使用圖類激活映射方法[30]找出不同圖譜（AAL116和CC200）對PTSD分類貢獻Top 10的腦區及對應腦區的MNI坐標（完整坐標數據參見附件3）。從Top 10腦區的MNI坐標可以發現，AAL116和CC200腦圖譜中有一個重疊腦區顯示為左顳極的左顳中回（Temporal_Pole_Mid_L，MTG），為了更直觀地說明，本文使用腦可視化工具BrainNet Viewer[31]分別顯示了兩個圖譜中Top 10的腦區，結果如圖5所示。其中，重疊腦區用紅色虛線框標出。
圖5
不同圖譜分類貢獻Top 10腦區
紅色虛線框表示不同圖譜重疊腦區
Figure5. Top 10 brain regions contributing to classification across different atlasesred dashed boxes indicate overlapping brain regions for different atlases
3 討論
由表3結果可見，AAL116圖譜和CC200圖譜的融合達到了最優的分類結果，這可能是因為AAL116腦圖譜是根據腦結構劃分的功能模版，與數據驅動的腦圖譜劃分方法CC200之間的FCN拓撲信息的互補是最佳的，這與文獻[24]報道是一致的。同時，添加更多腦圖譜并不能提高分類性能，甚至導致分類性能下降，這可能是因為多個腦圖譜中獲取的高維節點特征包含了更多冗余信息和噪聲，反而造成性能下降，這與文獻[25]報道結果一致。
從圖4結果發現，頻譜卷積對
的選取比較敏感，AUC值隨
的增大會呈現整體增加的趨勢。當ks1=ks2 = 4的組合時AUC結果取得最優，這可能是因為前一層和后一層的卷積在這里取得了輸入特征全局和局部信息的最佳平衡（圖4中黑框顯示）。同時，當ks1 = ks2 = 1的組合時，AUC結果并不理想（圖4中白框顯示），這與文獻[28]研究結果一致，進而說明了多核卷積的有效性。
從附件3結果發現，Top 10的腦區中，AAL116中有3個腦區（分別為：Heschl_L、Heschl_R和MTG）以及CC200中有5個腦區（分別為：Paracingulate Gyrus_R、Cingulate Gyrus_R、Precuneous Cortex_R、Planum Temporale_R和Superior Temporal Gyrus_L）都屬于默認模式網絡（default mode network，DMN），這與Zhu等[32]研究結果中DMN網絡與地震后PTSD高度相關是一致的。其次，剩余Top 10腦區分別位于背側注意網絡、軀體運動網絡、中央執行網絡、視覺網絡和邊緣網絡，也與文獻[17]發現一致。此外，本文發現的Top 10的腦區多分布于額葉和顳葉皮層，這與Zhang等[33]的研究一致。Geuze等[34]研究發現，與健康對照者相比，在編碼階段，PTSD患者表現出額葉皮層激活不足和顳葉皮層過度激活，同時，在研究過程中觀察到患有PTSD的退伍軍人的右額葉皮層、雙側顳中回和左海馬后/海馬旁回的激活減少，這種差異表明額葉顳葉網絡可能參與了PTSD中的認知缺陷，這些發現也進一步證實了本文發現額葉和顳葉皮層對PTSD分類具有較大貢獻這一結果。
此外，從圖5可以直觀地發現AAL116和CC200腦圖譜中的重疊腦區顯示為左顳極的MTG。華西醫院龔啟勇團隊2016年發表的論文[35]中發現，左腦島、紋狀體、杏仁核和顳中回的異常可能是長期創傷導致PTSD的原因。同時，在Herringa等[36]的研究中發現在持續行動中遭受戰斗創傷的PTSD患者的MTG觀察到灰質體積減少。本文絕大多數PTSD患者在經歷地震創傷之前沒有罹患PTSD的經驗，因此沒有意識到自己具有PTSD的癥狀，是通過反復的被動篩查后納入受試組；另外，此數據是在2008年汶川地震后7～9年采集，其結果可能與地震后長期創傷有關。綜上，MTG異常有望作為中長期PTSD患者較為客觀的影像標記，從而為臨床醫生提供進一步的參考。
4 結論
本文提出一種多圖譜MK-GCN的框架，基于靜息態fMRI數據進行PTSD的分類診斷。通過多個腦圖譜構建不同尺度的FCNs進行表征學習，從輸入層面獲取了多個圖譜FCN之間的互補信息。同時，通過改變傳統GCN層中卷積核的大小來捕獲FCN中類可分離的輸出特征，從而改善不同尺度下的FCNs的特征學習能力，使GCN模型可以選擇更優特征來解決圖分類問題。此外，本文從臨床角度出發，分析了PTSD患者和創傷后未患PTSD的健康受試之間的分類貢獻Top 10的腦區，并在不同腦圖譜的功能連接模式下發現了重疊腦區，有望作為診斷中長期PTSD更為可靠客觀的影像標記，為臨床醫生提供參考。但是，本文提出的多圖譜MK-GCN方法在識別PTSD患者的問題上雖取得了較好的實驗結果，仍存在一些局限性。首先，本研究在基于圖神經網絡的框架中只研究了110名受試者樣本，樣本量較小。其次，盡管本研究在單站點小樣本數據中取得了較好的實驗結果，但還沒有在多中心的腦影像數據中進行驗證。在針對PTSD患者的神經影像學研究中，單個中心較難獲得大量樣本。最近的研究趨勢傾向于提高神經影像數據的共享水平，通過多中心數據來驗證研究的可重復性。在未來的工作中，基于本文提出的診斷模型有望推廣至其他類型PTSD的研究中，并進一步探索PTSD相關的神經機制，從而開展更具臨床應用價值的研究。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：周麗君主要負責實驗流程、數據記錄與分析、論文編寫以及算法程序設計；劉云飛主要負責程序指導；莫賢主要負責影像分析和論文審閱；袁軍、羅昌宇主要負責整理數據和實驗結果；張俊然、朱鴻儒主要負責實驗指導、數據分析指導、論文審閱修訂。
倫理聲明：本研究通過了四川大學華西醫院醫學倫理委員會的審批［2018年 審（236）號］。
本文附件見本刊網站的電子版本（
0 引言
創傷后應激障礙（post-traumatic stress disorder，PTSD）是一種在創傷性事件之后較常見的精神障礙疾病，常見的創傷源為戰爭、自然災害（如地震）、車禍等。其特征是對創傷性事件的持續侵入性記憶、過度覺醒、回避與創傷相關的線索，以及思維和情緒的負面變化[1-2]。這種慢性精神疾病嚴重影響著患者的生活質量，并可造成相當大的經濟負擔[3]。由于PTSD的臨床表現復雜多樣，目前仍缺乏可靠的客觀生物學標志，單純依靠量表和醫生的主觀經驗來進行臨床評估和診斷可能造成準確性的偏差和一致性缺乏[4]，同時，由于客觀診斷手段（如影像標記）的可重復性差，亦會導致大量PTSD患者未被發現或誤診[5]。因此，迫切需要建立可靠、客觀、指標可重復的輔助診斷模型，通過客觀手段重現PTSD的臨床測量來彌合主觀評估和客觀測量之間的差距[6]，并促進未來診斷方法的改進。
目前，研究者主要基于功能磁共振成像（functional magnetic resonance imaging，fMRI）數據來建立大腦的功能連接網絡（functional connectivity network，FCN）等特征映射[7]。隨后，利用機器學習或深度學習方法，學習和選擇精神障礙疾病特征映射中相較于傳統方法（如稀疏學習等）更敏感、更有效的特異性指標，以提高精神障礙疾病診斷的準確性和可靠性[8-10]，并將參與分類的主要貢獻特征作為與疾病高度相關的影像標記。Campbell等[11]結合多個尺度特征（即區域均勻性、低頻波動的局部幅度和基于功能連接的網絡），訓練多層感知機對有意識清醒和麻醉誘導的無意識狀態進行分類，用于從靜息態fMRI數據中識別麻醉誘導的無意識狀態的穩健影像標記，取得了0.84的平均曲線下面積（area under the curve，AUC）。Chang等[12]使用深度堆疊自編碼器提取低頻振幅的低維特征識別主要精神疾病及其亞型，并推定了新的診斷影像標記，為精神疾病的精準診斷提供了啟發。值得注意的是，雖然已有大量研究將影像數據與深度學習方法結合起來分類精神分裂癥[13]、重度抑郁癥[14]、自閉癥[15]、阿爾茨海默病[16]等精神疾病的不同群體，但利用深度學習方法對PTSD進行分類診斷的研究還很有限。Zhu等[17]將深度學習技術結合基于圖論的分析指標識別PTSD患者獲得了80%的平均準確率，高于傳統支持向量機（57%）。Zhu等[18]使用自編碼器用于PTSD多站點數據的分類，取得了75%的AUC值。盡管上述方法在捕獲和提取與PTSD相關的大腦功能連接特征方面取得了一定的成果，能夠從高維fMRI數據中提取復雜的空間和時間特征，減少了對人工提取特征的依賴[6]，但是，目前使用深度學習等方法對PTSD的分類準確率還有待提高[17]，同時，這類方法對精神障礙疾病類別歸屬的可解釋性也有限[18]，這可能降低臨床醫生的信任度和接受度。
最近，圖神經網絡（graph neural network，GNN）已經成為一種強大的深度學習工具，被廣泛用于解決醫學fMRI圖像中非歐氏空間圖數據問題[19]。GNN還能夠將重要特征可視化，以反駁深度學習是“黑匣子”的典型批判，能夠用于揭示精神疾病的影像特征[20]。GNN通過將大腦的功能連接網絡表示為圖來建模大腦的復雜網絡拓撲結構，其中，節點對應腦模板定義的感興趣區域（regions of interest，ROI），邊量化了ROI之間的功能連接。GNN已在之前的研究中已成功應用于自閉癥譜系障礙[21]、阿爾茨海默病[22]、重度抑郁癥[23]等精神障礙疾病的表征，取得了顯著的成果，各項分類指標提高的同時增強了對精神疾病的影像標記的理解。然而，之前的研究發現，使用單一圖譜劃分ROI，忽略了不同空間尺度下FCN潛在的互補拓撲信息[24]，單一圖譜并不能完全捕捉大腦區域連接的復雜性和異質性[25-26]，這種限制可能會削弱大腦網絡的多尺度特征刻畫能力。此外，傳統圖卷積網絡（graph convolutional network，GCN）模型[27]在所有層中使用恒定大小的卷積核，導致了每個節點的特征以固定的跳數學習鄰居的特征，不利于獲取類可分離的輸出特征[28]，從而限制了模型的特征提取能力，導致分類性能受到影響。
鑒于上述問題，本文提出多圖譜多核圖卷積（multi-kernel graph convolution neural network，MK-GCN）模型用于PTSD的分類診斷。主要在兩個方面做出了改進：首先，本文采用多個腦圖譜劃分ROI，以構建多尺度的FCNs進行表征學習，從輸入層面獲取多個圖譜FCN之間的互補拓撲信息；其次，通過改變GCN層中卷積核的大小來捕獲FCN中類可分離的輸出特征。最后，本文將從臨床角度出發，分析PTSD患者和創傷后未患PTSD的健康受試者之間的分類貢獻腦區，并在不同腦圖譜的模板下進行重疊腦區的確認，期望更可靠地識別地震后PTSD患者異常的大腦區域，進而為臨床醫生提供參考。
1 材料和方法
本文方法主要分為多尺度圖構建、多核圖卷積特征提取、分類輸出等幾個部分，模型整體框架如圖1所示：主要包括不同圖譜功能連接矩陣構建、K最近鄰算法選擇邊構建圖、多核圖卷積提取特征、特征融合后全連接分類輸出。首先，本文采用先驗圖譜自動解剖標記（automated anatomical labeling，AAL116）和基于數據驅動的腦圖譜（craddock200，CC200）分別對預處理后的信號構建功能連接矩陣。然后，使用K最近鄰算法選擇邊后構建為圖數據。接著，將不同尺度的圖數據送入多核圖卷積網絡中進行模型訓練和特征提取，捕獲更優特征。最后，將每個尺度提取到的特征融合后輸出到全連接層完成二分類，實現地震后PTSD組和地震后未患PTSD（Non-PTSD）的健康對照組的識別。
圖1
地震后應激障礙預測整體框圖
Figure1.
Post-earthquake traumatic stress disorder classification and prediction framework
1.1 數據采集及預處理
1.1.1 受試者信息
本研究招募的138名參與者均為2008年5月12日汶川地震后的幸存者，最后納入2015年至2017年來自同一地震災區的受試，年齡18～60歲，均為右利手，未接受任何常規治療（藥物或心理治療）。研究中兩名精神科醫生根據《精神疾病診斷和統計手冊》中創傷后應激障礙量表（Clinically Administered PTSD Scale，CAPS）評估PTSD癥狀，核磁共振掃描和心理測試在入組后兩天內進行評估。受試者排除標準包括：有任何神經系統疾病史、嚴重頭部損傷且意識喪失時間超過10 min、雙相情感障礙或其他精神障礙病史的參與者；智力低下、懷孕、酗酒和/或其他物質依賴/濫用；起搏器或手術夾等金屬植入物的干擾，以及其他磁共振成像掃描禁忌等情況。本研究經四川大學華西醫院醫學倫理委員會批準，并獲得每位受試者的書面知情同意書。對采集后的圖像質量進行了篩選，最后保留了110例受試者的靜息態fMRI數據，受試者信息如表1所示。
表1
110名受試者的人口統計學信息
Table1.
Demographic information of the 110 subjects
最終受試者分為兩組：① 59名PTSD患者，其中男性16名，女性43名。② 51名地震后未患PTSD的健康受試者，其中男性16名，女性35名，作為對照組。
1.1.2 靜息態fMRI數據采集
本研究110名受試者的fMRI數據采集均在3.0T MRI掃描儀上使用12通道相控陣頭部線圈進行，使用梯度回波平面成像序列獲取fMRI血氧水平依賴性圖像，其參數設置如下：重復時間/回波時間=2 000 ms/30 ms；翻轉角=90°；切片厚度=5 mm；視野=240 mm×240 mm；體素大小=3.75 mm×3.75 mm×5 mm；采集矩陣=64×64。所有受試者的掃描時長均接近7 min，對應205個時間點。
1.1.3 數據預處理
使用神經影像腦連接組分析軟件（graph theoretical network analysis，GRETNA）對采集數據進行預處理。處理步驟包括：① 考慮到初始靜息態fMRI信號不穩定和參與者前期的適應時間，去除每個功能數據的前10個時間點，剩余195個時間點用于后續處理。② 剩余的數據首先進行時間層校正，然后重新對齊以去除頭動偽影的影響。將頭動位移超過2.5 mm或者頭動旋轉角度超過2.5°的受試者剔除。本研究中所有受試者的頭動位移和頭動旋轉角度均未超過上述標準。③ 回歸協變量后，使用蒙特利爾神經學研究所（Montreal Neurological Institute，MNI）的模板將得到的圖像配準至3 mm體素的標準空間。④ 最后對每個體素時間序列進行帶通濾波（0.01～0.08 Hz），以減少低頻漂移和高頻生理噪聲的影響。
1.2 模型描述
同一受試者預處理后的圖像使用AAL116和CC200進行不同尺度腦區劃分；同時分別建立功能連接矩陣后對兩個尺度的功能連接（functional connectivity，FC）矩陣構建圖并進行多核卷積圖表征學習，最后在特征融合后進行分類輸出。
1.2.1 多尺度圖構建
本文使用三元組表示具有N個節點的無向圖。其中是N個腦區的集合，邊的集合表示FC矩陣通過K最近鄰算法稀疏化后保留的連邊。表示節點特征矩陣，其中F表示節點特征維度，矩陣中的每個元素表示腦區i和j之間的皮爾遜相關系數值。通過上述操作，對于同一受試者得到兩個特證矩陣：和，同時獲得兩個鄰接矩陣和，最后，基于兩個空間尺度為每個受試構造兩個圖和。
1.2.2 多核卷積圖表征學習
使用譜圖卷積定義整個卷積過程來分析大腦的FCN，并提出多核卷積的圖表征方法診斷PTSD。圖上的傅立葉變換依賴于拉普拉斯矩陣的特征向量，歸一化后的圖拉普拉斯矩陣L可以表示為式(1)。
|
其中，表示單位矩陣，是度矩陣，在這里表示腦區連接的數量，表示由傅里葉基組成的特征矩陣，表示特征值的對角矩陣。則對于每個腦區的信號x的傅里葉變換就可以表示為，其逆變換就可以表示為。譜圖卷積定義為信號x與傅里葉域中可學習的濾波器，為了防止信號x的傅立葉變換矩陣乘法所需的巨大計算量，使用切比雪夫多項式參數化重新定義濾波器表示為式(2)：
|
其中，是k次切比雪夫系數的向量，同時，定義為一個函數，對信號進行頻譜濾波，以腦區i為中心的濾波器的腦區j的值可以表示為式(3)：
|
受文獻[28]啟發，當k大于0且為整數時，代表了k跳鄰居，則對于任意的腦區節點i和j則有式(4)成立。
|
其中，是和之間的最短路徑距離，是從到的最短路徑上的所有邊權重的總和。式（3）中拉普拉斯算子的k階多項式表示的譜濾波器恰好反映其k跳鄰域局部化表示。對于完整特征集X，使用濾波器的輸出可以表示為式(5)。
|
為了全面考慮腦區特征之間的緊密程度及其鄰域的情況，本文使用多個不同大小的濾波器（多核卷積）獲取對應的輸出特征，其中，表示單跳鄰居（最小）和6跳鄰居（最大）的內核大小，每個內核卷積后的輸出特征聚合后輸出至下一層。同時為了保證節點特征不被平滑化，選用3層GCN模型提取特征如式(6)所示。
|
其中，表示第l層神經網絡的權重矩陣，，，表示歸一化鄰接矩陣和單位矩陣的和，是非線性激活函數。輸出特征，d表示最終輸出的特征維度。
最后，采用Lee等[29]提出的全局平均池化操作聚合節點特征并讀出同一受試的最終圖級表征，可表示為式(7)：
|
其中，表示卷積運算得到的第i個腦區的特征向量。整個多核卷積和圖讀出過程用圖2進行描述。
圖2
多核圖卷積網絡詳細描述
Figure2.
Detailed description of multi-kernel graph convolutional networks
1.2.3 特征融合與分類輸出
通過卷積操作，特征提取后最終獲取了每個尺度的128維的圖表示向量，分別表示為和，M表示受試者個數。如圖1所示，本文使用拼接操作進行特征融合，同時采用自適應權重策略自動獲取每個尺度的特征融合比例。最后，通過特征融合得到每個受試的兩個尺度特征更新后的256維圖表示向量，將此特征進行分類得出預測結果，并用圖類激活映射找出分類貢獻腦區（Top 10）。
1.3 實驗設置
本文模型基于PyTorch和PyG實現。3個圖卷積層的神經元個數分別為128、64、32，使用ReLU非線性激活和dropout操作防止模型過擬合，dropout設置為0.3。采用余弦退火學習率動態調整學習率，其初始值為0.001。同時采用帶標簽平滑的交叉熵損失作為優化函數，模型參數使用網格搜索進行優化，訓練迭代次數為100次。
實驗中進行了5折交叉驗證。評價指標包括分類精度、靈敏度、特異性。同時，繪制了受試者工作特性（receiver operating characteristic，ROC）曲線并計算了AUC，曲線的橫坐標為假陽性率，縱坐標為真陽性率。ROC曲線的AUC可以定量評價PTSD組和Non-PTSD組分類的準確性。本文采用雙樣本t檢驗進行統計學分析，所有統計分析的檢驗水準設定為0.05。顯著性檢驗結果的解釋以P值為準，P值小于0.05被認為具有統計學意義。數據分析使用Python 3.8.10進行。
1.4 對比方法
本文提出的MK-GCN方法與傳統機器學習方法——支持向量機（support vector machine，SVM）和邏輯回歸（logistic regression，LR），以及經典GCN模型[27]進行了對比實驗，對比模型具體實驗設置及參數參見附件1。
2 實驗結果
2.1 分類結果
在PTSD數據上進行對比實驗，對比結果如表2所示。同時，圖3直觀地顯示了這幾種方法的ROC曲線。從表2和圖3可看出以下結果：首先，在四個評價指標和ROC曲線上，本文提出的MK-GCN_F方法在PTSD和Non-PTSD分類中都取得了最優結果；其次，與單尺度結果相比，MK-GCN_F在四項評價指標上提高了至少2%；此外，對于單一尺度，使用CC200腦圖譜劃分ROI比AAL116腦圖譜劃分ROI取得了更好的結果，同時，使用GCN的方法比機器學習方法獲得了更優的分類性能，而多核圖卷積MK-GCN獲得了比傳統GCN方法至少高4%的AUC。這些結果表明，與使用單一腦圖譜劃分ROI相比，使用多尺度腦圖譜有助于提高PTSD分類性能，其原因可能是不同尺度FCNs特征包含互補的拓撲信息，這些互補信息的協同可以改善分類結果；此外，通過改變圖卷積過程濾波器的大小，可以改善GCN的特征提取能力，更好地捕捉類可分離的輸出特征，從而進一步提高PTSD分類性能。為了進一步驗證本文模型的的泛化能力，在公開的自閉癥譜系障礙（autism spectrum disorder，ASD）數據集ABIDE上取了最大站點NYU進行了驗證，此站點共有受試184人，ASD 79人，健康受試者105人；其中ASD患者女性11人、男性68人，健康受試者女性26人、男性79人，年齡、性別、教育時間均無顯著差異。實驗結果表明，本文模型在ASD中的分類結果相較于對比模型各項指標也有明顯提升，分類精度、靈敏度、特異性至少提升了4%，取得了86.95%的AUC值，完整數據參見附件2。
表2
不同方法在PTSD數據中評價指標對比結果（%）
Table2.
Comparison of evaluation metrics across different methods for PTSD data (%)
圖3
PTSD分類受試者工作特征曲線
Figure3.
Receiver operating characteristic curve for PTSD classification
2.2 不同腦圖譜分類結果
本研究還采用Power264（基于元分析將皮質區域映射到整個皮層的更精細的功能劃分圖譜）和Zalesky980（基于全腦解剖連接建模的更精細的圖譜）來觀察不同腦圖譜對分類結果的影響。實驗結果發現，與AAL116和CC200腦圖譜相比，更精細的腦圖譜對PTSD分類性能反而有所下降。為了更全面地評估多圖譜融合對分類結果的影響，本文還研究了圖譜兩兩組合、3個圖譜組合、4個圖譜組合對PTSD分類結果的影響，實驗結果見表3。
表3
不同腦圖譜及融合后對比結果（%）
Table3.
Different brain atlases and post-fusion comparison results (%)
2.3 不同大小卷積核分類結果
為了研究不同卷積核的大小對每個尺度輸入圖性能的影響，設定在相同的輸入（即相同的節點特征和圖結構表示）下，ks1和ks2[1, 6]的分類性能，其中，ks1表示第二層GCN的內核大小，ks2表示第三層GCN的內核大小，ks = 1和ks = 6分別表示了單跳鄰居（最小）和6跳鄰居（最大）的內核大小。對比結果見圖4（以CC200為例）。
圖4
不同層卷積核ks內核大小對分類結果影響熱圖
Figure4.
Heatmap of the effect of different layers of convolutional kernel kernel size of ks on classification results
2.4 分類貢獻Top 10腦區及可視化結果
本文使用圖類激活映射方法[30]找出不同圖譜（AAL116和CC200）對PTSD分類貢獻Top 10的腦區及對應腦區的MNI坐標（完整坐標數據參見附件3）。從Top 10腦區的MNI坐標可以發現，AAL116和CC200腦圖譜中有一個重疊腦區顯示為左顳極的左顳中回（Temporal_Pole_Mid_L，MTG），為了更直觀地說明，本文使用腦可視化工具BrainNet Viewer[31]分別顯示了兩個圖譜中Top 10的腦區，結果如圖5所示。其中，重疊腦區用紅色虛線框標出。
圖5
不同圖譜分類貢獻Top 10腦區
紅色虛線框表示不同圖譜重疊腦區
Figure5. Top 10 brain regions contributing to classification across different atlasesred dashed boxes indicate overlapping brain regions for different atlases
3 討論
由表3結果可見，AAL116圖譜和CC200圖譜的融合達到了最優的分類結果，這可能是因為AAL116腦圖譜是根據腦結構劃分的功能模版，與數據驅動的腦圖譜劃分方法CC200之間的FCN拓撲信息的互補是最佳的，這與文獻[24]報道是一致的。同時，添加更多腦圖譜并不能提高分類性能，甚至導致分類性能下降，這可能是因為多個腦圖譜中獲取的高維節點特征包含了更多冗余信息和噪聲，反而造成性能下降，這與文獻[25]報道結果一致。
從圖4結果發現，頻譜卷積對的選取比較敏感，AUC值隨的增大會呈現整體增加的趨勢。當ks1=ks2 = 4的組合時AUC結果取得最優，這可能是因為前一層和后一層的卷積在這里取得了輸入特征全局和局部信息的最佳平衡（圖4中黑框顯示）。同時，當ks1 = ks2 = 1的組合時，AUC結果并不理想（圖4中白框顯示），這與文獻[28]研究結果一致，進而說明了多核卷積的有效性。
從附件3結果發現，Top 10的腦區中，AAL116中有3個腦區（分別為：Heschl_L、Heschl_R和MTG）以及CC200中有5個腦區（分別為：Paracingulate Gyrus_R、Cingulate Gyrus_R、Precuneous Cortex_R、Planum Temporale_R和Superior Temporal Gyrus_L）都屬于默認模式網絡（default mode network，DMN），這與Zhu等[32]研究結果中DMN網絡與地震后PTSD高度相關是一致的。其次，剩余Top 10腦區分別位于背側注意網絡、軀體運動網絡、中央執行網絡、視覺網絡和邊緣網絡，也與文獻[17]發現一致。此外，本文發現的Top 10的腦區多分布于額葉和顳葉皮層，這與Zhang等[33]的研究一致。Geuze等[34]研究發現，與健康對照者相比，在編碼階段，PTSD患者表現出額葉皮層激活不足和顳葉皮層過度激活，同時，在研究過程中觀察到患有PTSD的退伍軍人的右額葉皮層、雙側顳中回和左海馬后/海馬旁回的激活減少，這種差異表明額葉顳葉網絡可能參與了PTSD中的認知缺陷，這些發現也進一步證實了本文發現額葉和顳葉皮層對PTSD分類具有較大貢獻這一結果。
此外，從圖5可以直觀地發現AAL116和CC200腦圖譜中的重疊腦區顯示為左顳極的MTG。華西醫院龔啟勇團隊2016年發表的論文[35]中發現，左腦島、紋狀體、杏仁核和顳中回的異常可能是長期創傷導致PTSD的原因。同時，在Herringa等[36]的研究中發現在持續行動中遭受戰斗創傷的PTSD患者的MTG觀察到灰質體積減少。本文絕大多數PTSD患者在經歷地震創傷之前沒有罹患PTSD的經驗，因此沒有意識到自己具有PTSD的癥狀，是通過反復的被動篩查后納入受試組；另外，此數據是在2008年汶川地震后7～9年采集，其結果可能與地震后長期創傷有關。綜上，MTG異常有望作為中長期PTSD患者較為客觀的影像標記，從而為臨床醫生提供進一步的參考。
4 結論
本文提出一種多圖譜MK-GCN的框架，基于靜息態fMRI數據進行PTSD的分類診斷。通過多個腦圖譜構建不同尺度的FCNs進行表征學習，從輸入層面獲取了多個圖譜FCN之間的互補信息。同時，通過改變傳統GCN層中卷積核的大小來捕獲FCN中類可分離的輸出特征，從而改善不同尺度下的FCNs的特征學習能力，使GCN模型可以選擇更優特征來解決圖分類問題。此外，本文從臨床角度出發，分析了PTSD患者和創傷后未患PTSD的健康受試之間的分類貢獻Top 10的腦區，并在不同腦圖譜的功能連接模式下發現了重疊腦區，有望作為診斷中長期PTSD更為可靠客觀的影像標記，為臨床醫生提供參考。但是，本文提出的多圖譜MK-GCN方法在識別PTSD患者的問題上雖取得了較好的實驗結果，仍存在一些局限性。首先，本研究在基于圖神經網絡的框架中只研究了110名受試者樣本，樣本量較小。其次，盡管本研究在單站點小樣本數據中取得了較好的實驗結果，但還沒有在多中心的腦影像數據中進行驗證。在針對PTSD患者的神經影像學研究中，單個中心較難獲得大量樣本。最近的研究趨勢傾向于提高神經影像數據的共享水平，通過多中心數據來驗證研究的可重復性。在未來的工作中，基于本文提出的診斷模型有望推廣至其他類型PTSD的研究中，并進一步探索PTSD相關的神經機制，從而開展更具臨床應用價值的研究。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：周麗君主要負責實驗流程、數據記錄與分析、論文編寫以及算法程序設計；劉云飛主要負責程序指導；莫賢主要負責影像分析和論文審閱；袁軍、羅昌宇主要負責整理數據和實驗結果；張俊然、朱鴻儒主要負責實驗指導、數據分析指導、論文審閱修訂。
倫理聲明：本研究通過了四川大學華西醫院醫學倫理委員會的審批［2018年 審（236）號］。
本文附件見本刊網站的電子版本（
