功能連接方法及其在fMRI數(shù)據(jù)中的應(yīng)用

文章來源于微信公眾號(茗創(chuàng)科技)尽狠,歡迎有興趣的朋友搜索關(guān)注响鹃。

導(dǎo)讀

強大的非侵入性神經(jīng)成像技術(shù)的可用性引發(fā)了各種旨在繪制人類大腦的研究油狂。這些研究不僅聚焦于發(fā)現(xiàn)大腦激活信號趟济,還聚焦于理解大腦網(wǎng)絡(luò)中功能性交流的整體結(jié)構(gòu)戴差∷桶郑基于不同的大腦區(qū)域在功能上相互連接并不斷共享信息的原理,文獻中提出了尋找這些功能網(wǎng)絡(luò)的各種方法。在這篇文章中袭厂,研究者概述了在fMRI數(shù)據(jù)中估計和描述功能連接的最常用方法墨吓。研究者用人類連接組項目的靜息態(tài)fMRI數(shù)據(jù)說明了這些方法,闡述了這些方法的實施細節(jié)和對結(jié)果解釋的見解纹磺。本研究的目標是通過提供必要的工具來測量和描述大腦回路帖烘,從而指導(dǎo)神經(jīng)成像領(lǐng)域的新研究人員。


1 引言

功能磁共振成像(fMRI)技術(shù)已成為描述人類大腦連接及其與健康橄杨、行為和生活方式關(guān)系的有力工具秘症。fMRI測量包括基于血氧水平依賴(BOLD)對比的間接的和非侵入性的大腦活動測量。與正電子發(fā)射斷層掃描(PET)和腦電圖(EEG)等替代性腦成像方式相比式矫,fMRI具有無創(chuàng)性和高空間分辨率的特點乡摹,這使得其成為大量腦成像研究的熱門選擇。例如采转,人類連接項目聪廉,該項目旨在通過描述健康成人大腦的連接模式來理解大腦的潛在功能。

這類研究主要有兩個目標:第一故慈,識別大腦中對外部刺激做出反應(yīng)的位置信號板熊;第二,識別大腦在靜息態(tài)或任務(wù)態(tài)時出現(xiàn)的大腦時-空關(guān)聯(lián)模式察绷。這些關(guān)聯(lián)模式是解剖學(xué)上不同腦區(qū)的功能連接時間序列中共同激活的度量干签,稱為功能連接。有證據(jù)表明拆撼,這些連接模式中的個體差異可以解釋個體在認知和行為功能上的重要差異容劳。因此,了解這些模式可以在預(yù)測神經(jīng)退行性疾病的早期發(fā)病以及監(jiān)測疾病護理和治療方面發(fā)揮重要作用情萤。

功能性MRI數(shù)據(jù)通常是高維的鸭蛙,由一段時間內(nèi)收集的3D腦體積圖像組成。在典型的fMRI研究中筋岛,體素Nv的數(shù)量是數(shù)十萬級的娶视,時間點T的數(shù)量也是成百上千的。估算體素空間連接性的Nv×Nv相關(guān)矩陣是一項具有挑戰(zhàn)性的任務(wù)睁宰,因此在具體研究中進行估算是需要一些策略和前期的假設(shè)肪获。一種簡單的方法是首先預(yù)先選定稱為“種子”的腦區(qū),然后計算種子間的相互關(guān)系和大腦中其他體素的功能時間序列柒傻。這種基于種子的方法計算簡單并且容易解釋孝赫,所以受到研究者的廣泛歡迎。種子的大小可以不同红符,大小可以和單個體素一樣小青柄。如果種子是一個區(qū)域伐债,則通常將某段時間進程內(nèi)的信號進行平均,并將其與響應(yīng)時間段內(nèi)所有其他體素相關(guān)致开。為了提高可測量性峰锁,通常首先將大腦分成幾個小的腦區(qū),并使用這些區(qū)域的平均時間序列來對腦網(wǎng)絡(luò)進行估計双戳。當比較精神疾病患者和正常人大腦的運作模式時虹蒋,基于種子的方法是非常有用的辦法。例如飒货,Wang等人(2006)使用基于種子的方法表明魄衅,與對照組被試相比,海馬體和其他腦區(qū)之間的連接會隨著阿爾茨海默病的早期癥狀而改變塘辅。盡管基于種子的方法很受歡迎晃虫,但是研究者對這種方法仍有各種各樣的批評。首先扣墩,因為預(yù)先確定了特定的種子傲茄,所以其他不同腦區(qū)出現(xiàn)的潛在模式會被忽略。其次沮榜,該方法是對某段時間序列內(nèi)的ROI進行平均而忽略了體素之間的可變性。第三喻粹,該方法只計算了節(jié)點對之間的相關(guān)性蟆融,而忽略了整個網(wǎng)絡(luò)中其他節(jié)點的潛在影響。

與預(yù)先指定特定腦區(qū)相比守呜,降維法則通過使用少量潛在成分來表示數(shù)據(jù)型酥,從而描述空間和/或時間連接模式。主成分分析(PCA)和獨立成分分析(ICA)是兩種最常用的方法查乒。這兩種方法都將高維成像數(shù)據(jù)投射到低維子空間弥喉。在PCA中,這種投射由正交成分組成玛迄,這些正交成分使投射到低維子空間的數(shù)據(jù)的方差最大化由境。在ICA中,投射由盡可能在空間上獨立的成分組成蓖议。然后虏杰,這些成分被組合成具有體素值的地形圖,每個體素值表示特定體素的相對數(shù)量勒虾,受該成分的激活調(diào)節(jié)纺阔。與基于種子的方法相比,PCA和ICA都具有自動化提取成分的優(yōu)勢修然,無需預(yù)先指定種子區(qū)域笛钝,也就是說這些方法是數(shù)據(jù)驅(qū)動的质况。Smith等人(2009)使用ICA將大腦網(wǎng)絡(luò)分解為空間子網(wǎng)絡(luò),這些子網(wǎng)絡(luò)在靜息態(tài)和任務(wù)態(tài)fMRI數(shù)據(jù)中具有相似的功能玻靡。

還有一些方法將基于種子的方法中分割腦區(qū)的策略與降維方法相結(jié)合结榄,來描述大腦神經(jīng)元回路。Ting等人(2020)使用解剖圖預(yù)先確定興趣區(qū)(ROIs)啃奴,然后通過主成分分析的方法從每個聚類中提取特征潭陪。然后使用多個提取的成分來估計這些興趣區(qū)之間的連接性,這種連接性使用RV系數(shù)來表示最蕾,是對特征集之間相關(guān)性的度量依溯。

除了用于估計功能連接的方法外,研究者通常還使用圖論工具來描述功能網(wǎng)絡(luò)瘟则。在圖論中黎炉,大腦網(wǎng)絡(luò)被視為由邊連接的節(jié)點的集合。通常將腦區(qū)之間的功能連接定義為連邊醋拧。對于指定的節(jié)點慷嗜,通過對功能連接矩陣閾值化,會得到一個二元矩陣丹壕。然后庆械,這個二元圖譜被用來計算描述大腦網(wǎng)絡(luò)性質(zhì)的各種圖形參數(shù)。這些參數(shù)代表了大腦網(wǎng)絡(luò)的一些關(guān)鍵特征菌赖,通常包括有助于確定圖節(jié)點是以隨機秩序連接還是以小世界秩序連接的數(shù)量缭乘。隨機網(wǎng)絡(luò)的連接模式更具有全局性,而小世界網(wǎng)絡(luò)則表現(xiàn)出較高的局部秩序性琉用。統(tǒng)計學(xué)的網(wǎng)絡(luò)模型將這些圖論測量作為預(yù)測具有多個個體特征的全腦網(wǎng)絡(luò)的輸入堕绩。

本文的目的是概述在靜息態(tài)功能磁共振成像數(shù)據(jù)中估計和描述功能連接性最常用的方法。研究者通過分析人類連接組項目中單個被試的數(shù)據(jù)來說明這些方法邑时。雖然研究者并不打算詳盡地介紹該領(lǐng)域快速發(fā)展的方法奴紧,但希望本文提供的信息能夠指導(dǎo)神經(jīng)成像領(lǐng)域的新研究人員探索這些數(shù)據(jù)。

接下來晶丘,本文的內(nèi)容包括以下幾個方面:

在第2節(jié)中黍氮,研究者描述了評估功能連通性的不同方法,主要針對單個被試的數(shù)據(jù)浅浮。在第3節(jié)中滤钱,使用第2節(jié)中描述的方法,來測量人類連接組項目中單個被試靜息態(tài)數(shù)據(jù)的功能連接脑题。在第4節(jié)中件缸,介紹了幾種針對多個被試功能連接估計的方法。在第5節(jié)中叔遂,描述了統(tǒng)計學(xué)網(wǎng)絡(luò)模型他炊。最后争剿,在第6節(jié)中對功能連接進行總結(jié)。


2 功能連接的方法

在本節(jié)中痊末,研究者將回顧估計單個被試腦功能連接的不同方法蚕苇。并在第4節(jié)中討論群體腦功能連接性。對于所有計算凿叠,研究者假設(shè)矩陣Y是大小為T×Nv的矩陣涩笤,由表示被試在每個體素v=1,….., Nv上的BOLD信號的Nv時間進程組成。這里盒件,為了簡單起見蹬碧,研究者通過將每個體素數(shù)據(jù)(列)減去其相應(yīng)時間過程的平均值來對矩陣Y中心化。無論是在靜息態(tài)還是任務(wù)態(tài)中功能連接分析的目標是描述大腦各區(qū)域如何相互作用炒刁。

2.1 基于種子的分析

計算大量體素之間所有兩兩成對的相關(guān)性在計算上代價是非常大的恩沽,因為這需要

次分析∠枋迹基于種子的分析(SBA)依賴于估計感興趣區(qū)域(ROI)之間或種子區(qū)域與大腦中所有其他體素之間的成對相關(guān)性罗心。

為了估計ROI之間的相關(guān)性,一種常見的方法是根據(jù)解剖模板劃分大腦區(qū)域城瞎,通常稱為腦圖譜渤闷。目前比較常用的腦圖譜有AAL、Tailarach圖譜和MN1-152圖譜脖镀。為了估計種子區(qū)域和其他體素之間的相關(guān)性肤晓,研究者通常會根據(jù)專家意見來選擇種子或選擇在功能磁共振成像實驗中表現(xiàn)出更大激活的體素作為種子。后一種方法在涉及任務(wù)的實驗中更常見认然。選擇種子后,通過量化相關(guān)性來估計連接性漫萄。以往的文獻中提出了各種量化的方法卷员。研究者在附錄A中提供了關(guān)于這些量化方法的更多細節(jié)√谖瘢基于種子的方法匯總?cè)缦拢?/p>

①選擇種子腦區(qū)或體素毕骡;

②將腦區(qū)或體素的時間序列與大腦中的所有其他體素計算相關(guān)。如果種子是一個腦區(qū)岩瘦,則在將該腦區(qū)的時間序列與大腦中的所有其他體素計算相關(guān)之前未巫,對該腦區(qū)的時間序列求平均值。使用附錄A中所述的方法之一启昧;

③呈現(xiàn)相關(guān)性的3D像或顯示相關(guān)性閾值(僅顯示達到顯著性水平的相關(guān))叙凡。注:為了降低假陽的概率,確定顯著性密末,因而需要考慮多重比較握爷。Bonferroni和FDR是兩種廣泛使用的多重比較校正的方法跛璧。

或者,在使用腦圖譜將大腦劃分為不同的ROI后新啼,可以通過對體素進行平均或計算第一主成分來總括該腦區(qū)的時間序列追城。接下來,研究者使用這些總括的時間序列來與其他所有區(qū)域進行相關(guān)燥撞。研究者在第3節(jié)中對這兩種選擇都進行了說明座柱。

2.2 分解方法

雖然基于種子的方法在解釋時比較簡單,但在種子選擇技術(shù)上容易存在偏差物舒,所以在使用時需要謹慎色洞。主成分分析和獨立成分分析旨在通過提供功能連接的數(shù)據(jù)驅(qū)動方法來解決這個問題。這些分解方法在功能性神經(jīng)成像應(yīng)用中發(fā)揮著重要作用茶鉴。它們在預(yù)處理步驟中可以用于去除數(shù)據(jù)偽影和降低數(shù)據(jù)維度锋玲,并且它們可能會不止出現(xiàn)在一個估計不同人群中的功能連接性步驟中。在本節(jié)中涵叮,研究者將重點介紹它們在估計單個被試功能磁共振成像數(shù)據(jù)中功能連接方法中的作用惭蹂,而在第4節(jié)中,將探討它們在多被試分析中的貢獻割粮。

作為基于種子分析的替代方法盾碗,分解方法的目標是將體素域表示為空間分量的較小子集。每個空間分量都有一個單獨的時間過程舀瓢,代表許多體素的功能磁共振信號的同時變化廷雅。在本節(jié)中,研究者假設(shè)矩陣Y的每一列數(shù)據(jù)都減去了平均值京髓。

2.2.1 主成分分析(PCA)

PCA是一種常用的方法航缀,可以降低數(shù)據(jù)維數(shù),同時最大限度地減少數(shù)據(jù)信息的損失堰怨,最大限度地提高數(shù)據(jù)的可變性芥玉。主成分可以通過兩種方式來獲得:1.通過樣本協(xié)方差矩陣YTY的特征分解。2.通過使用奇異值分解(SVD)理論尋找數(shù)據(jù)矩陣Y的特征向量备图。

數(shù)據(jù)矩陣的秩是r=min{T灿巧,Nv}(通常是T<Nv和r=T),因此可以通過分解找到r個主成分

其中T×r矩陣U包含一個正交左奇異向量uk∈RT揽涮、r×Nv矩陣W包含正交右奇異向量wk∈RNv抠藕,且r×r對角矩陣∑包含有序奇異值。注意蒋困,YTY的特征分解定義為WT∑2W盾似。r×Nv矩陣W的正交行被稱為特征圖像,形成大腦地形圖雪标,每個地形圖代表由該成分激活所調(diào)節(jié)的給定體素的相對數(shù)量颜说。

另一種方法是將原始的功能磁共振成像數(shù)據(jù)投射到第一個p主成分所跨越的空間中购岗,其中p的選擇基于該成分能夠解釋數(shù)據(jù)可變異性指標。投射數(shù)據(jù)矩陣YW由這個新子空間中的區(qū)域時間序列組成门粪。Zhou等人(2009)利用這一想法降低某些ROI中功能磁共振成像數(shù)據(jù)的維數(shù)喊积,然后對兩個大腦區(qū)域的區(qū)組時間序列進行格蘭杰因果關(guān)系分析,以推斷方向連接玄妈。雖然可以使用這些投射數(shù)據(jù)點的時間序列計算相關(guān)性乾吻,但由于新子空間中的每個空間區(qū)域代表原始數(shù)據(jù)空間中不同體素的線性組合,因此結(jié)果沒有明確的解釋拟蜻。

2.2.2 獨立成分分析(ICA)

ICA旨在使用獨立因素的潛在表征來代表大腦數(shù)據(jù)绎签。與PCA不同,其目標是將Y分解為混合矩陣和空間獨立成分(IC)組合的結(jié)果酝锅。

式中诡必,M是帶有mk列的T×K混合矩陣,K×Nv矩陣C是帶有ck行的獨立成分矩陣搔扁,其中每個ck包含對應(yīng)于成分K的大腦網(wǎng)絡(luò)爸舒,總共有K個獨立成分。這些成分代表各種功能的網(wǎng)絡(luò)稿蹲,如運動扭勉、視覺、聽覺等苛聘。矩陣E是獨立的高斯噪聲涂炎。

假設(shè)成分映射(component maps),ck设哗,k=1, ....., K代表可能的重疊和統(tǒng)計學(xué)上的相關(guān)信號勾徽,但單個成分圖譜分布是獨立的嘱根,即面殖,如果P(ck)表示第K成分圖譜中體素值的概率分布峻贮,則

每個獨立成分ck是一個Nv大小的向量,可以形成大腦地形圖澎粟。與PCA一樣,這些地形圖表示給定體素的相對數(shù)量欢瞪,該體素由該成分的激活調(diào)節(jié)活烙。

2.3 計算方面

在腦成像應(yīng)用中,估計主成分需要對Nv×Nv矩陣YTY進行分解遣鼓,這通常是不可行的啸盏。文獻中提出了許多算法來有效地估計這種高維環(huán)境中的成分。Zou等人(2006)開發(fā)了稀疏PCA(SPCA)骑祟,它基于使用lasso懲罰回歸優(yōu)化估計過程(這種lasso懲罰是指高維環(huán)境下的多級函數(shù)主成分)通過迭代閾值算法估計主成分的稀疏集回懦。其中一些使用的一些工具包很容易獲取气笙,可以在作者的網(wǎng)站上下載。

類似地怯晕,估計獨立成分并不簡單潜圃,對成分進行排序也很困難,因為IC通常不是正交的舟茶,每個成分解釋的方差之和不能與原始數(shù)據(jù)的方差之和相加谭期。最早的算法之一是Infomax,其目的是最大化適當變換的成分映射的聯(lián)合熵吧凉。最近隧出,更現(xiàn)代的算法側(cè)重于從包含大量特征的數(shù)據(jù)矩陣中提取稀疏的特征集。例如Le等人(2011)提出的具有重建成本(RICA)的獨立成分分析(ICA)阀捅,在Matlab中有相應(yīng)的工具箱胀瞪。

2.4 一種混合方法

Ting等人(2020)給出了估算功能連接的不同方法。作者提出了一個基于多拓撲尺度網(wǎng)絡(luò)的多尺度模型饲鄙,從體素水平到由體素聚類組成的腦區(qū)凄诞,以及由這些腦區(qū)的集合組成的更大網(wǎng)絡(luò)。在實際應(yīng)用中傍妒,這些體素集合是預(yù)先指定的幔摸,通常被視為解剖學(xué)上的ROI。然后颤练,這些ROI可以組合成ROI聚類既忆。他們的方法包括一個維度縮減步驟,直到每個ROI中包含一個因子模型嗦玖。設(shè)r表示第r個ROI患雇,Yr為T×pr數(shù)據(jù)矩陣,由屬于第r個ROI的體素時間序列組成(包含總共的pr體素宇挫,其中

pr=Nv苛吱,R為ROI總數(shù))。然后器瘪,可以得到這樣一個公式

其中翠储,Yr(t)是一個大小為pr的列向量,fr(t)是一個潛在公共因子的mr×1向量橡疼,其中有許多因子mr<<pr援所,Qr是一個mr×pr因子加載矩陣,通過fr的混合定義pr體素之間的依賴關(guān)系欣除,而Er(t)=[er1(t)住拭,…,erpr(t)]'是白噪聲的pr×1向量,E(Er(t))=0滔岳,

=Cov(Er(t))=diag(

杠娱,…,

)谱煤。然后將這些因子模型連接起來進行定義

其中Y(t)是大小為

pr=Nv的列向量摊求,Q=diag(Q1,…趴俘,QR)是一個

pr×

mr塊對角混合矩陣睹簇,f(t)=[fr(t),…寥闪,fR(t)]'是一個

mr×1聚合潛在因素向量太惠。

這些不同拓撲尺度下的網(wǎng)絡(luò)協(xié)方差矩陣使用低秩矩陣按以下方式進行估計。設(shè)

為ROI r內(nèi)的協(xié)方差矩陣疲憋。模型(4)包含以下分解

類似地凿渊,從模型(5)中,可以得到

低維因子協(xié)方差矩陣∑ff是一個區(qū)組矩陣缚柳,用于估計ROI之間的零滯后相關(guān)(lag-zero dependency)埃脏,如下所示。

對角線塊∑fr fr秋忙,r=1,....R是對角協(xié)方差矩陣彩掐,用于捕獲每個ROI內(nèi)因素的總方差。非對角塊∑fk fj灰追,j≠j是因子之間的交叉協(xié)方差矩陣堵幽,總結(jié)了聚類j和k之間的相關(guān)性。

作者使用RV系數(shù)總結(jié)了ROI之間的相關(guān)性弹澎,RV系數(shù)是平方相關(guān)系數(shù)的多元推廣朴下。聚類j和k中因子之間的RV系數(shù)定義如下:

其中

在實際估計中苦蒿,作者應(yīng)用該模型來估計靜息態(tài)腦網(wǎng)絡(luò)殴胧。他們使用PCA估計因子fr和矩陣Qr,并根據(jù)解剖學(xué)上的腦圖譜預(yù)先指定ROI佩迟。Miranda等人(2021)利用人類連接組項目中工作記憶任務(wù)的數(shù)據(jù)团滥,使用這種方法來估計ROI之間的功能連接。

2.5 腦網(wǎng)絡(luò)

使用圖論中的工具來表示大腦是很常見的报强。在這個框架中灸姊,我們可以將功能連接視為一個由圖表示的網(wǎng)絡(luò),其中空間單元是節(jié)點躺涝,它們之間的連接是邊厨钻。網(wǎng)絡(luò)被視為由連邊連接的節(jié)點(頂點)的集合。圖(網(wǎng)絡(luò))表示為G=(V坚嗜,E)對夯膀,其中V和E分別是頂點和邊的集合。此外苍蔬,圖形可以加權(quán)诱建,在這種情況下,將用三重G=(V碟绑,E俺猿,W)表示,W(E)表示每條邊的權(quán)重格仲。

首先要做的是選擇網(wǎng)絡(luò)的節(jié)點押袍。與基于種子的估計連接類似,這些節(jié)點由體素或腦圖譜給出的ROI細分定義凯肋。根據(jù)節(jié)點的規(guī)格谊惭,必須確定它們的邊(連接)。這些邊量化了這些不同節(jié)點之間的關(guān)聯(lián)強度侮东,即功能連接圈盔。前面討論的功能連接和附錄A中描述測量方法用來量化邊的強度。

圖論分析的大多數(shù)標準工具都是由二進制網(wǎng)絡(luò)開發(fā)的悄雅,其中每條邊要么存在要么不存在驱敲。通常稱為鄰接矩陣的二元矩陣是通過對連通矩陣進行閾值化得到的。雖然使用標準圖論機制對加權(quán)圖進行閾值設(shè)置很方便宽闲,但在這個過程中众眨,原始信號的信息會有損失。此外便锨,在大多數(shù)情況下围辙,閾值的選擇不是唯一的,很難說明這樣的操作是合理的放案。一種策略是使用大規(guī)模單變量方法姚建,在這種方法中,對網(wǎng)絡(luò)中的每個可能的連邊進行統(tǒng)計測試吱殉,然后使用標準技術(shù)(如Bonferroni校正或錯誤發(fā)現(xiàn)率(FDR))進行多重比較校正掸冤。

在對腦網(wǎng)絡(luò)進行估計后,計算了一些描述性的度量友雳,作為描述拓撲圖性質(zhì)的手段稿湿。在估計腦網(wǎng)絡(luò)是常用的指標有特征路徑長度、集群系數(shù)和度分布押赊。有關(guān)神經(jīng)成像中使用的綜合拓撲度量的列表饺藤,請參見Rubinov等人(2010)的研究。

典型路徑長度(Characteristic path length)。路徑是不同節(jié)點的序列涕俗,代表路徑較短的成對大腦區(qū)域之間的潛在信息流罗丰,這意味著整合的潛力更強。路徑的長度估計了大腦區(qū)域之間功能整合的可能性再姑。功能整合最常用的度量之一是網(wǎng)絡(luò)中所有節(jié)點對之間的平均最短路徑長度萌抵,它被定義為典型路徑長度。斷開連接的節(jié)點之間的路徑被定義為具有無限長元镀,這是計算該度量時的一個問題绍填,尤其是在稀疏網(wǎng)絡(luò)中,例如在函數(shù)連接性中栖疑。實際上讨永,研究者只在現(xiàn)有路徑之間取平均值,這可能存在一定的問題遇革。關(guān)于這個問題的討論住闯,請參考Rubinov等人(2010)的研究。

節(jié)點度分布澳淑。作為中心性的一種度量比原,單個節(jié)點的度等于連接到該節(jié)點的連邊數(shù),即該節(jié)點的鄰居數(shù)杠巡。因此量窘,節(jié)點度分布就是網(wǎng)絡(luò)中所有節(jié)點的度分布。在功能連接中氢拥,高度節(jié)點度與網(wǎng)絡(luò)中的許多其他節(jié)點進行功能性交互蚌铜,稱為節(jié)點中心。

集群系數(shù)嫩海。作為分離的一種度量冬殃,集群系數(shù)是節(jié)點的鄰居中彼此也是鄰居的部分,在圖論中是單個節(jié)點周圍三角形的部分叁怪。功能網(wǎng)絡(luò)中的集群系數(shù)表明审葬,一種統(tǒng)計相關(guān)性的組織,表示分離的功能性神經(jīng)加工奕谭,即在緊密相連的大腦區(qū)域群中進行特定加工的能力涣觉。平均而言,網(wǎng)絡(luò)的平均集群系數(shù)反映了單個節(jié)點周圍的聚集連接的普遍程度血柳。平均集群系數(shù)是對每個節(jié)點分別進行歸一化官册,并且會受到度值較低節(jié)點的不成比例的影響。

許多其他網(wǎng)絡(luò)指標受基本網(wǎng)絡(luò)特征的影響很大难捌,例如節(jié)點和連邊的數(shù)量以及本節(jié)介紹的分布程度膝宁。


3 真實數(shù)據(jù)示例

研究者分析了來自人類連接組項目(HCP)的靜息態(tài)數(shù)據(jù)鸦难。選擇使用FIX pipeline處理已降噪的數(shù)據(jù)。使用高通時間濾波器员淫,執(zhí)行運動回歸(即24個運動參數(shù)的回歸:六個剛體運動參數(shù)明刷、它們的后向時間導(dǎo)數(shù)和這12個時間序列的平方),并應(yīng)用基于ICA的回歸满粗,以消除與信號成分正交的噪聲成分的方差。對于單被試分析愚争,研究者考慮了從示例的右-左階段(受試者100307)收集的像映皆。每720ms獲得一個fMRI像。每個像由大小為91×109×91的圖像組成轰枝,共有1200個采樣點捅彻。

3.1 單被試示例

3.1.1 基于ROI的功能連接分析

使用MNI152空間的AAL Atlas版本的腦圖譜劃分ROI鞍陨。數(shù)據(jù)共涉及166個ROI步淹,并使用以下方法估計功能鏈接:(a)每個ROI中平均時間序列的交叉相關(guān);(b)每個ROI中平均時間序列的偏相關(guān)诚撵;(c)將ROI數(shù)據(jù)的時間序列投射到其第一主成分空間的交叉缭裆;(d)將ROI數(shù)據(jù)的時間序列投射到其第一主成分空間的偏相關(guān);(e)對于每個ROI寿烟,研究者考慮占ROI變異性20%的主要成分澈驼,并計算方程(8)中描述的RV系數(shù)。

估計的功能連接結(jié)果如圖1所示筛武。觀察圖1可以發(fā)現(xiàn)缝其,圖片(a)和(c)中的交叉相關(guān)比其相應(yīng)的偏相關(guān)(圖(b)和(d))捕捉到更大的相關(guān)性。圖(e)中使用的RV系數(shù)方法似乎能夠捕捉到大量ROI之間的低相關(guān)徘六。在從圖中得出任何結(jié)論之前内边,應(yīng)該首先檢驗這些值是否顯著。對于前四個矩陣待锈,在進行檢驗時漠其,首先將這些值轉(zhuǎn)換為z分數(shù),然后對它們進行閾值化竿音,以確定重要(顯著)的相關(guān)辉懒。對于圖(e)中的RV系數(shù),顯著性是基于Ting等人(2020)研究中詳述的系數(shù)的漸近分布來確定的谍失。

圖1.基于AAL分區(qū)的ROI的功能連接眶俩。圖(a)描述了ROI的平均時間序列的交叉相關(guān),圖(b)描述了ROI的平均時間序列的偏相關(guān)快鱼,圖(c)描述了與第一偏相關(guān)投射的ROI數(shù)據(jù)的時間序列的交叉相關(guān)颠印,圖(d)描述了與第一偏相關(guān)投射的ROI數(shù)據(jù)的時間序列的偏相關(guān)纲岭,圖(e)表示RV系數(shù),每個ROI保留了能夠解釋其20%變異性的主成分线罕。

接下來止潮,研究者使用這些功能連接矩陣可以獲得一個二元圖,其連邊是根據(jù)從相關(guān)系數(shù)的z分數(shù)中獲得的p值確定的钞楼,如附錄a中方程(A2)所述喇闸。p值根據(jù)Bonferroni校正進行了閾值化,顯著性為5%询件。對于圖(e)中的RV系數(shù)燃乍,使用系數(shù)的漸近分布將值轉(zhuǎn)換為z分數(shù),并基于Bonferroni校正進行閾值化宛琅,以找到顯著性為5%的標準正態(tài)分布分位數(shù)刻蟹。考慮到這個標準嘿辟,研究者計算了如圖2所示的鄰接矩陣舆瘪。

觀察圖2可以發(fā)現(xiàn)(a)和(c)圖都存在大量連邊。這表明不同興趣區(qū)之間存在高水平的相互作用红伦。在圖(b)和(d)中沒有發(fā)現(xiàn)這種高水平的相互作用英古。在圖(e)中,可以觀察到一些ROI與許多其他ROI之間存在中等程度的相互作用昙读,而在靜息態(tài)實驗期間哺呜,一些區(qū)域是安靜的。

圖2.從圖1的閾值連接矩陣得到的二元圖箕戳。

3.1.2 腦網(wǎng)絡(luò)總結(jié)性測量

研究者利用上面獲得的二元圖某残,使用第2.5節(jié)中所述的圖論方法來估計一些總結(jié)性測量。各計算中使用的公式詳見附錄B陵吸。表1顯示了計算結(jié)果玻墅。CPL是排除網(wǎng)絡(luò)中所有無限路徑的特征路徑長度,DG是腦網(wǎng)絡(luò)的平均度壮虫,其中度表示每個節(jié)點中的連邊數(shù)澳厢,CC是腦網(wǎng)絡(luò)的平均集聚系數(shù),Inf是網(wǎng)絡(luò)中無限路徑的數(shù)量囚似。表1中的數(shù)值反映了研究者在圖2中觀察到的情況剩拢。度表示興趣區(qū)之間的連接數(shù)。正如前面所注意到的饶唤,圖(a)和(c)中的圖說明ROI之間存在很多交互作用徐伐,表明度值很高。RV系數(shù)的特征路徑長度(CPL)表明募狂,平均的腦網(wǎng)絡(luò)路徑長度較短办素,其值與圖2(a)和(c)中的腦網(wǎng)絡(luò)相當角雷。這表明,盡管連接的區(qū)域很少性穿,但連接的區(qū)域彼此相鄰勺三。

表1 腦網(wǎng)絡(luò)總結(jié)性測量


3.1.3 基于體素的功能連接

基于種子的分析。對于基于種子的分析需曾,研究者選擇左側(cè)額下回島蓋部作為種子點吗坚。獲取該腦區(qū)的平均時間序列,并計算與剩余體素的交叉相關(guān)呆万。研究者進行了Bonferroni校正商源,α=0.05,來確定相關(guān)性的閾值桑嘶。圖3顯示了生成的大腦地形圖。研究者將大于125的簇(聚類類別)視為重要體素躬充,將他們覆蓋在標準腦上逃顶,標準腦由這里使用的示例被試數(shù)據(jù)的平均時間點組成。

圖3.左側(cè)額下回島蓋部的基于種子分析的功能連接充甚。圖中顯示了帶有體素的矢狀位切片以政,這些體素與紅色顯示的種子ROI有顯著的聯(lián)系。

分解方法伴找。首先獲得數(shù)據(jù)矩陣Y的主成分盈蛮。需要注意的是,需要大量的主成分來表示數(shù)據(jù)的可變異性技矮,而傳統(tǒng)的主成分具有第2.3節(jié)中討論的問題抖誉。對于這一特定數(shù)據(jù),需要150個成分來解釋20%的數(shù)據(jù)可變性衰倦,需要463個成分來解釋50%的數(shù)據(jù)可變性袒炉。研究者在圖4中展示了前五個按特征值(即載荷)縮放(降維)的成分。

接下來樊零,為了估計獨立成分我磁,研究者使用Beckmann等人(2004)研究中提出的概率獨立成分分析,并通過FSL中的MELODIC(多元探索線性優(yōu)化分解為獨立成分)函數(shù)實現(xiàn)驻襟。圖5描述了結(jié)果夺艰。

為了便于說明,研究者在不設(shè)置閾值的情況下呈現(xiàn)成分沉衣。在多被試方法中郁副,更常見的是使用單個成分的映射作為輸入,然后在群組成分中執(zhí)行閾值化豌习,以識別群組腦網(wǎng)絡(luò)霞势。

圖4.從第一個(頂部)到第五個(底部)的有序主成分圖的矢狀位視圖烹植。
圖5.根據(jù)從第一個(頂部)到第五個(底部)不斷增加的唯一解釋方差對獨立成分的地形圖進行排序的矢狀位視圖。


4 多被試的功能連接

在對功能性MRI建模時愕贡,一個重要的目標是通過利用跨被試的共享結(jié)構(gòu)來識別多被試數(shù)據(jù)中的功能連接結(jié)構(gòu)草雕。多被試功能連接模型的范圍可以從約束張量分解模型(constrained tensor decomposition models)(例如PARAFAC)到更靈活的方法,其中首先估計特定被試的連接矩陣或PCA和ICA模型固以,并將其串聯(lián)結(jié)果(concatenated results)用作基于組水平估計的輸入墩虹。最佳模型將取決于哪一級別的靈活性水平最能捕獲組水平的功能連接特征。

在多被試的ICA模型中憨琳,一個簡單的過程是使用預(yù)先指定的ROI(如第2節(jié)中描述的基于種子的方法)估計單被試的功能連接诫钓,然后將這些結(jié)果聚合到單個矩陣中,隨后使用主成分進一步對該矩陣進行分解(decomposing)篙螟。然后菌湃,可以映射主成分以估計基于組水平的功能連接。

Caffo等人(2010)的研究中采用了一種多階段方法來比較阿爾茨海默病高家族風(fēng)險被試(臨床無癥狀)與匹配對照組之間的功能連接遍略。該方法遵循四個步驟惧所,包括特定被試的SVD、匯總特定被試特征向量的組水平降維绪杏、將單被試水平的數(shù)據(jù)投射到總體特征向量上以獲得特定被試的因子載荷下愈,以及在功能性邏輯回歸模型中使用特定被試的因子載荷。

許多研究提出了一種組水平ICA方法蕾久,其中势似,fMRI數(shù)據(jù)要么跨被試暫時性連接,要么作為多維數(shù)組僧著。FMRIB軟件庫(FSL)是一個包含各種成像數(shù)據(jù)的圖像分析和統(tǒng)計工具的軟件庫履因,在其MELODIC函數(shù)中提供了組ICA和張量ICA(tensorial ICA)。本節(jié)將重點介紹這兩種方法盹愚。

4.1 組水平ICA

Calhoun等人(2001)的研究中首次提出了對一組被試的fMRI數(shù)據(jù)進行ICA的方法搓逾。假設(shè)我們觀察n名被試的功能磁共振成像數(shù)據(jù)。假設(shè)Yi是一個大小為T×Nv的矩陣杯拐,由Nv個時間序列組成霞篡,代表每個體素v=1,...端逼,Nv處的BLOD信號朗兵,被試i=1,...顶滩,n余掖。他們的模型包括以下多個階段。

①單個被試水平數(shù)據(jù)的歸約(reduction)礁鲁。在這一步中盐欺,在時域中應(yīng)用歸約(reduction)赁豆。對于每個被試i=1,...冗美,n魔种。歸約的數(shù)據(jù)如下所示:

其中

是L×T歸約矩陣,Xi是代表歸約數(shù)據(jù)的L×Nv矩陣粉洼。實際上节预,F(xiàn)?1通過PCA分解(decomposition)得到;

②匯總單被試水平數(shù)據(jù)的歸約(reduction)属韧。對NL×Nv矩陣[

安拟,…,

]T進行歸約宵喂,以獲得K×Nv矩陣X=G?1[

糠赦,…,

]T锅棕,其中K是要獲得的組水平成分數(shù)量拙泽,G?1是一個K×NL歸約化矩陣,實際上是由主成分得到的哲戚;

③獨立成分的估計奔滑。ICA分解應(yīng)用于矩陣X艾岂,如第2.2.2節(jié)所述顺少。

其中M是K×K混合矩陣,C是K×Nv成分映射矩陣王浴。得到的組水平ICA成分可以通過首先將其轉(zhuǎn)換為Z分數(shù)來設(shè)定閾值脆炎。

可以通過對單個被試劃分GM(其中G=(G?1)T)獲得個體水平的地形圖煞茫,并按照如下所示返回前面的步驟烟央。

基于這些步驟,矩陣GMC是單個映射的大小為NL×Nv的矩陣枚碗,可以進行分區(qū)钞啸,因此GiMiCi=

Yi几蜻,Ci包含個體水平的地形圖。

4.2 張量ICA

張量ICA是基于張量分解的体斩,它可以得到多維數(shù)組的低秩表達梭稚。PARAFAC是一種常見的分解方法。設(shè)X∈RT×Nv×N是一個陣列絮吵,分別包含功能磁共振成像數(shù)據(jù)和時間弧烤、體素和被試。三維陣列X可以按以下方式分解為R外積之和

其中

蹬敲。這種分解意味著數(shù)組X中的每個元素都可以寫成暇昂。

在分解中的向量可以用矩陣表示莺戒,例如A= [a1a2 ...aR],同樣地急波,也可以得到矩陣B和C从铲。三維數(shù)組可以在一個稱為矩陣化的過程中展開為矩陣。展開可以發(fā)生在三個維度中的任何一個幔崖。在第二維度上食店,X(2)∈RNv×NT是X的二模矩陣化。類似地赏寇,我們可以延伸生成二模和三模矩陣吉嫩。

關(guān)于PARAFAC分解和矩陣化的詳細信息,請參考Kolda等人(2009)的研究嗅定。根據(jù)這些定義自娩,可以將公式寫為:

,其中☉表示Katri–Rao乘積(product)渠退。Beckmann等人(2005)的研究中提出了方程的ICA分解

忙迁,其中

、混合矩陣

以及成分矩陣BT的估計方法詳見Beckmann等人(2004)的研究碎乃。


5 統(tǒng)計學(xué)網(wǎng)絡(luò)模型

在本節(jié)中姊扔,研究者按照Solo等人(2018)研究中的符號來描述統(tǒng)計網(wǎng)絡(luò)模型,目的是描述大腦網(wǎng)絡(luò)梅誓。在這些模型中恰梢,首先使用第2節(jié)中的技術(shù)估計單個被試的功能連接。估計完單個被試的數(shù)據(jù)之后梗掰,考慮多個感興趣的變量和拓撲網(wǎng)絡(luò)特征對整體網(wǎng)絡(luò)結(jié)構(gòu)的影響嵌言。

假設(shè)(Yi,Xi)分別代表被試i的腦網(wǎng)絡(luò)和協(xié)變量及穗。給定協(xié)變量的網(wǎng)絡(luò)的概率密度函數(shù)用P(Yi | Xi摧茴,θi)表示,其中θi描述Yi和Xi的關(guān)系埂陆。這些協(xié)變量可以是特定于節(jié)點的協(xié)變量苛白,如大腦位置,也可以是網(wǎng)絡(luò)的功能焚虱,如路徑長度或第2.5節(jié)中描述的其他指標购裙。常用的密度函數(shù)建模方法包括指數(shù)隨機圖模型(ERGM)和混合模型。

在ERGM中著摔,研究者認為二元圖和模型適用于每個被試缓窜,如下所示。設(shè)Yi是一個由R×R節(jié)點組成的網(wǎng)絡(luò)。然后禾锤,如果節(jié)點j和k之間存在連接私股,則Yijk=1,否則Yijk=0恩掷。概率質(zhì)量函數(shù)具有正則指數(shù)族的形式:

參數(shù)θ的估計使用MCMC ML方法完成倡鲸。在Simpson等人(2012)的研究中確定了每個受試者網(wǎng)絡(luò)中最重要的解釋性指標g(yi)。接下來黄娘,為所有被試創(chuàng)建一個基于組的擬合參數(shù)值θ的匯總測量峭状。他們使用這些基于組水平的解釋性度量和參數(shù),通過ERGM來擬合基于組水平的代表性網(wǎng)絡(luò)逼争。

目前的ERGM估計方法的一個局限性是其可伸縮性(scalability)优床。主要問題不是ROI本身的數(shù)量,而是腦網(wǎng)絡(luò)的連邊結(jié)構(gòu)誓焦,這可能會導(dǎo)致沒有辦法收斂的問題胆敞。此外,大多數(shù)模型都是為二元圖開發(fā)的杂伟,不太適合對連邊水平上的問題進行檢驗移层。作為ERGM的替代方案,混合模型允許連邊水平的檢查以及多被試的比較赫粥。該框架定義了一個由兩部分組成的混合效應(yīng)观话,對連接存在或不存在的概率以及連接強度(如果存在)進行建模。假設(shè)Yi表示附錄a中列出的一個相關(guān)度量給出的功能連接強度越平,讓Rijk表示j和k之間是否存在連接频蛔。那么條件概率是

其中,βr是固定效應(yīng)的向量喧笔,該向量與每個參與者和節(jié)點對的協(xié)變量Xijk相關(guān)帽驯,而bri是代表特定被試和節(jié)點的參數(shù)的隨機效應(yīng)龟再。

假設(shè)Zijk為與隨機效應(yīng)bri相關(guān)的設(shè)計矩陣书闸;模型分為兩部分。模型的第一部分使用logit鏈接函數(shù)將節(jié)點j和k之間的連接概率與協(xié)變量聯(lián)系起來利凑,如下所示浆劲。

第二部分通過將節(jié)點i和j之間的相關(guān)系數(shù)的Fisher Z變換與協(xié)變量線性關(guān)聯(lián),在存在連接的情況下哀澈,對節(jié)點j和k之間的連接強度進行建模牌借。讓Sijk=Yijk | Rijk=1,則

其中割按,βr是一個總體參數(shù)向量膨报,它與每個被試和節(jié)點對的同一組協(xié)變量Xijk的連接強度相關(guān),bsi是一個被試和節(jié)點特定參數(shù)向量,它捕捉了這種關(guān)系如何隨總體平均βs變化现柠,eijk是主體i和節(jié)點j和k的隨機噪聲院领。Simpson等人(2015)給出了兩部分建模方法的詳細信息。

這些模型存在一個問題够吩,即使考慮了多重比較比然,基于連接權(quán)重的閾值選擇將影響網(wǎng)絡(luò)拓撲。Solo等人認為周循,持久同源性提供了一個解決閾值問題的多尺度層次框架强法。該方法是一種計算拓撲技術(shù),為比較網(wǎng)絡(luò)提供了一致的數(shù)學(xué)框架湾笛。持久同源性記錄了拓撲網(wǎng)絡(luò)特征在多分辨率和多尺度上的變化饮怯,而不是以固定閾值觀察網(wǎng)絡(luò)。通過這樣做嚎研,它揭示了對噪聲具有魯棒性的特征硕淑,即最“持久”的拓撲特征。


6 總結(jié)

本文回顧了fMRI數(shù)據(jù)中最常用的功能連接估計的方法嘉赎。對于單被試數(shù)據(jù)置媳,可以通過直接量化感興趣區(qū)域和/或種子區(qū)域之間的相關(guān)性,或通過找到一組代表同時活動的潛在成分來進行估計公条。前一種方法的在解釋上比較簡單,而后一種方法的解釋則不那么清楚靶橱。在提供的示例中寥袭,表示數(shù)據(jù)可變異性所需的成分圖數(shù)量非常多,因此关霸,僅對少數(shù)成分的調(diào)查可能無法反映大腦網(wǎng)絡(luò)的整體情況传黄。第2節(jié)中獲得的結(jié)果表明,即使以等效的方式定義區(qū)域队寇,不同的連接估計程序也會導(dǎo)致對腦網(wǎng)絡(luò)的不同解釋膘掰。因此,意識到每種方法的局限性是非常重要的佳遣,尤其是在解釋單個被試數(shù)據(jù)的結(jié)果時识埋。

盡管單被試分析面臨挑戰(zhàn),但單被試分析的方法也可以成功地應(yīng)用于多被試功能連接分析零渐。尤其是如果這個方法不需要多個階段窒舟,而是執(zhí)行張量ICA框架中的聯(lián)合估計。其他新興的多被試腦網(wǎng)絡(luò)的方法诵盼,如持久同源性(persistent homology)惠豺,是估計大腦回路的一種很有前途的方法银还,尤其是可伸縮性(scalability)可以實現(xiàn)的話。

原文:Functional Connectivity Methods and Their Applications in fMRI Data.

https://doi.org/10.3390/e24030390

文章來源于微信公眾號(茗創(chuàng)科技)洁墙,歡迎有興趣的朋友搜索關(guān)注见剩。

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 序言:七十年代末,一起剝皮案震驚了整個濱河市扫俺,隨后出現(xiàn)的幾起案子苍苞,更是在濱河造成了極大的恐慌,老刑警劉巖狼纬,帶你破解...
    沈念sama閱讀 206,126評論 6 481
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件羹呵,死亡現(xiàn)場離奇詭異,居然都是意外死亡疗琉,警方通過查閱死者的電腦和手機冈欢,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 88,254評論 2 382
  • 文/潘曉璐 我一進店門,熙熙樓的掌柜王于貴愁眉苦臉地迎上來盈简,“玉大人凑耻,你說我怎么就攤上這事∧停” “怎么了香浩?”我有些...
    開封第一講書人閱讀 152,445評論 0 341
  • 文/不壞的土叔 我叫張陵,是天一觀的道長臼勉。 經(jīng)常有香客問我邻吭,道長,這世上最難降的妖魔是什么宴霸? 我笑而不...
    開封第一講書人閱讀 55,185評論 1 278
  • 正文 為了忘掉前任囱晴,我火速辦了婚禮,結(jié)果婚禮上瓢谢,老公的妹妹穿的比我還像新娘畸写。我一直安慰自己,他們只是感情好氓扛,可當我...
    茶點故事閱讀 64,178評論 5 371
  • 文/花漫 我一把揭開白布枯芬。 她就那樣靜靜地躺著,像睡著了一般幢尚。 火紅的嫁衣襯著肌膚如雪破停。 梳的紋絲不亂的頭發(fā)上翅楼,一...
    開封第一講書人閱讀 48,970評論 1 284
  • 那天尉剩,我揣著相機與錄音,去河邊找鬼毅臊。 笑死理茎,一個胖子當著我的面吹牛,可吹牛的內(nèi)容都是我干的。 我是一名探鬼主播皂林,決...
    沈念sama閱讀 38,276評論 3 399
  • 文/蒼蘭香墨 我猛地睜開眼朗鸠,長吁一口氣:“原來是場噩夢啊……” “哼!你這毒婦竟也來了础倍?” 一聲冷哼從身側(cè)響起烛占,我...
    開封第一講書人閱讀 36,927評論 0 259
  • 序言:老撾萬榮一對情侶失蹤,失蹤者是張志新(化名)和其女友劉穎沟启,沒想到半個月后忆家,有當?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體,經(jīng)...
    沈念sama閱讀 43,400評論 1 300
  • 正文 獨居荒郊野嶺守林人離奇死亡德迹,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點故事閱讀 35,883評論 2 323
  • 正文 我和宋清朗相戀三年芽卿,在試婚紗的時候發(fā)現(xiàn)自己被綠了。 大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片胳搞。...
    茶點故事閱讀 37,997評論 1 333
  • 序言:一個原本活蹦亂跳的男人離奇死亡卸例,死狀恐怖,靈堂內(nèi)的尸體忽然破棺而出肌毅,到底是詐尸還是另有隱情筷转,我是刑警寧澤,帶...
    沈念sama閱讀 33,646評論 4 322
  • 正文 年R本政府宣布悬而,位于F島的核電站旦装,受9級特大地震影響,放射性物質(zhì)發(fā)生泄漏摊滔。R本人自食惡果不足惜阴绢,卻給世界環(huán)境...
    茶點故事閱讀 39,213評論 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一處隱蔽的房頂上張望艰躺。 院中可真熱鬧呻袭,春花似錦、人聲如沸腺兴。這莊子的主人今日做“春日...
    開封第一講書人閱讀 30,204評論 0 19
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽页响。三九已至篓足,卻和暖如春,著一層夾襖步出監(jiān)牢的瞬間闰蚕,已是汗流浹背栈拖。 一陣腳步聲響...
    開封第一講書人閱讀 31,423評論 1 260
  • 我被黑心中介騙來泰國打工, 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留没陡,地道東北人涩哟。 一個月前我還...
    沈念sama閱讀 45,423評論 2 352
  • 正文 我出身青樓索赏,卻偏偏與公主長得像,于是被迫代替她去往敵國和親贴彼。 傳聞我的和親對象是個殘疾皇子潜腻,可洞房花燭夜當晚...
    茶點故事閱讀 42,722評論 2 345

推薦閱讀更多精彩內(nèi)容