A Mixed Model Approach for Estimating Regional Functional Connectivity from Voxel-level BOLD Signals

Technische Samenvatting: Een Mixed-Model Benadering voor het Schatten van Regionale Functionele Connectiviteit vanuit Voxel-niveau BOLD-signalen

Probleemstelling
Rustende-toestand functionele connectiviteit (FC) kwantificeert de synchronie tussen activiteitspatronen van anatomisch gescheiden hersenregio's. In functionele magnetische resonantiefotografie (fMRI) worden gegevens verkregen als bloedzuurstofgehalte-afhankelijke (BOLD) signalen over talrijke ruimtelijke voxels. De gangbare methode voor het schatten van inter-regionale connectiviteit is de Correlation of Averages (CA) estimator, die voxel-niveau signalen binnen een regio aggregeert (meestal via ruimtelijke middeling) en de Pearson-correlatie tussen deze regionale gemiddelden berekent.

Hoewel de CA-estimator computationeel eenvoudig is, lijdt deze aan een aantoonbare bias en een gebrek aan statistisch gefundeerde rechtvaardiging. Zoals aangetoond door Achard et al. (2023), convergeert de CA-estimator naar een waarde die niet alleen afhangt van de ware inter-regionale correlatie, maar ook van nuisance-parameters, waaronder intra-regionale ruimtelijke correlatie, signaal-ruisverhoudingen en het aantal voxels. Bijgevolg kunnen CA-schattingen ernstig vertekend zijn, met name wanneer ruimtelijke correlaties hoog zijn of de signaalsterkte laag is. Bovendien hebben bestaande voxel-niveau modelleringsbenaderingen (bijv. Kang et al., 2012; Castruccio et al., 2018) zich primair gericht op taak-gebaseerde activatie of groep-niveau inferentie, waarbij ze vaak falen om intra-regionale van inter-regionale variabiliteit duidelijk te onderscheiden of om robuuste individuele-niveau netwerkschattingen te bieden.

Methodologie
De auteurs stellen een rigoureus statistisch kader voor op basis van lineaire mixed-effects modellen om inter-regionale functionele connectiviteit direct te schatten vanuit voxel-niveau BOLD-signalen.

1. Spatiotemporeel Mixed-Effects Model:
   Het geobserveerde BOLD-signaal bij voxel $v_{jl}$ in regio $R_j$ en wavelet-coëfficiënt $m$ wordt gemodelleerd als:
   $$X_{jlm} = \mu_j + \eta_{jm} + \gamma_{jlm} + \epsilon_{jlm}$$

   • $\mu_j$: Regionaal gemiddeld niveau.
   • $\eta_{jm}$: Een random effect met een gemiddelde van nul dat de regionale signaalcomponent representeert, wat inter-regionale afhankelijkheid induceert via de correlatiematrix $R = \{\rho_{jj'}\}$.
   • $\gamma_{jlm}$: Een lokaal random effect met een gemiddelde van nul dat intra-regionale ruimtelijke afhankelijkheid vastlegt, ongecorreleerd met $\eta_{jm}$.
   • $\epsilon_{jlm}$: Ongecorreleerde meetruis.
     Het model opereert in het wavelet-domein om schaal-specifieke temporele afhankelijkheid te isoleren en hoogfrequente ruis te verminderen. De covariantiestructuur gaat uit van scheidbaarheid tussen de ruimtelijke en wavelet-domeinen, waarbij gebruik wordt gemaakt van Matérn- en Gaussische kernels.

2. Tweestaps Schattingsprocedure:
   Om de computationele onhandelbaarheid van volledige maximum likelihood-schatting (MLE) voor hoog-dimensionale spatiotemporele data aan te pakken (waarbij $N$ de 20.000 voxels kan overschrijden), is een tweestaps benadering ontwikkeld:

   • Fase 1: Restricted Maximum Likelihood (ReML) wordt onafhankelijk toegepast op de data van elke regio om intra-regionale covariantieparameters (ruimtelijke correlatie, ruisvariantie) te schatten. Deze stap behandelt de regionale random effects als vaste effecten om de lokale structuur te isoleren.
   • Fase 2: De inter-regionale correlatieparameters ($\rho_{jj'}$) worden geschat door te conditioneren op de Fase 1-schattingen. Voor paren van regio's maken de auteurs gebruik van de Vecchia's likelihood benadering (Vecchia, 1988; Guinness, 2018). Deze methode benadert de gezamenlijke likelihood als een product van conditionele distributies gebaseerd op een ijle set van dichtstbijzijnde buren, waardoor de computationele complexiteit wordt teruggebracht van $O(N^3)$ naar een hanteerbaar niveau terwijl de schattingsnauwkeurigheid behouden blijft.

3. Inferentie:
   De asymptotische eigenschappen van de estimators worden vastgesteld onder een regime van een uitbreidend tijdsvenster (divergerend aantal wavelet-coëfficiënten). Benaderde betrouwbaarheidsintervallen worden geconstrueerd met behulp van de Fisher-informatie matrix afgeleid van de Vecchia-benadering.

Belangrijkste Bijdragen

• Nieuwe Estimator: Het artikel introduceert een op maximum likelihood gebaseerde estimator voor inter-regionale connectiviteit afgeleid van een mixed-effects model, waarmee het verder gaat dan de moment-gebaseerde CA-estimator.
• Expliciete Modellering van Variabiliteit: Het framework modelleert expliciet en scheidt bronnen van variabiliteit: inter-regionale correlatie, intra-regionale ruimtelijke afhankelijkheid en meetruis. Dit maakt een onbevooroordeelde schatting van connectiviteitsparameters mogelijk die niet worden verward met artefacten van regionale signaalaggregatie.
• Computationele Schaalbaarheid: Dit werk presenteert de eerste toepassing van Vecchia's likelihood benadering op functionele connectiviteit met voxel-niveau fMRI-data, wat de analyse van massieve datasets mogelijk maakt die voorheen computationeel onhaalbaar waren voor volledige likelihood-methoden.
• Verbeterde Reproduceerbaarheid: De methode is ontworpen om consistentere individuele hersennetwerken te produceren, waarmee het probleem van de test-retest betrouwbaarheid in fMRI wordt aangepakt.

Resultaten

• Simulatiestudies: Uitgebreide simulaties die de voorgestelde ReML- en Vecchia-estimators vergelijken met CA, de Average of Covariances (ACE) en EBLUE estimators, tonen een superieure prestatie aan. De voorgestelde estimators vertonen aanzienlijk lagere bias en gemiddelde kwadratische fout, met name in scenario's met een hoge intra-regionale ruimtelijke correlatie of lage signaal-ruisverhoudingen. Cruciaal is dat de voorgestelde methode een nominale dekking bereikt voor 95% betrouwbaarheidsintervallen, terwijl CA-gebaseerde intervallen systematisch falen om de ware parameter te dekken.
• Toepassing op de Human Connectome Project (HCP): De methode werd toegepast op test-retest data van 42 jonge volwassenen. Connectiviteitsnetwerken werden geconstrueerd voor zowel de voorgestelde methode als de CA-estimator. Gebruikmakend van de Concordance Correlation Coefficient (CCC) om de overeenstemming tussen herhaalde scans te meten (Exam 0 vs. Exam 1), leverde de voorgestelde mixed-model benadering consequent een hogere concordantie op dan CA. Dit geeft aan dat de voorgestelde methode meer betrouwbare en reproduceerbare individuele hersennetwerken produceert.

Betekenis
Het artikel voert aan dat het wijdverbreide gebruik van de CA-estimator een systematische bias introduceert in functionele connectiviteitsstudies, wat potentieel de validiteit van groepsvergelijkingen en individuele karakterisering beïnvloedt. Door een statistisch gefundeerde, mixed-model benadering te bieden die rekening houdt met voxel-niveau afhankelijkheden en ruis, bieden de auteurs een methode die nauwkeurigere en betrouwbaardere schattingen van de hersennetwerkarchitectuur oplevert. De aangetoonde verbetering in test-retest reproduceerbaarheid suggereert dat deze benadering kan leiden tot een robuuster wetenschappelijk begrip van individuele variabiliteit in connectiviteitsnetwerken en hun associaties met gezondheidsuitkomsten. De auteurs merken op dat hoewel de methode is ontworig voor fMRI, het onderliggende model ook toepasbaar kan zijn in andere contexten met ruimtelijk gegroepeerde tijdreeksen.