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

Technische Zusammenfassung: Ein Mixed-Model-Ansatz zur Schätzung der regionalen funktionellen Konnektivität aus BOLD-Signalen auf Voxel-Ebene

Problemstellung
Die funktionelle Konnektivität im Ruhezustand (Resting-state Functional Connectivity, FC) quantifiziert die Synchronität zwischen den Aktivitätsmustern anatomisch getrennter Hirnregionen. In der funktionellen Magnetresonanztomographie (fMRT) werden Daten als blutsauerstoffabhängige (BOLD-) Signale über zahlreiche räumliche Voxel hinweg erfasst. Die gängige Methode zur Schätzung der inter-regionalen Konnektivität ist der „Correlation of Averages“ (CA)-Schätzer, der die Signale auf Voxel-Ebene innerhalb einer Region (typischerweise durch räumliche Mittelung) aggregiert und die Pearson-Korrelation zwischen diesen regionalen Mittelwerten berechnet.

Obwohl der CA-Schätzer rechentechnisch einfach ist, leidet er unter nachweisbarer Verzerrung (Bias) und einem Mangel an statistisch fundierter Begründung. Wie Achard et al. (2023) zeigen, konvergiert der CA-Schätzer gegen einen Wert, der nicht nur von der wahren inter-regionalen Korrelation abhängt, sondern auch von Störparametern wie der intra-regionalen räumlichen Korrelation, dem Signal-Rausch-Verhältnis und der Anzahl der Voxel. Folglich können CA-Schätzungen stark verzerrt sein, insbesondere wenn die räumlichen Korrelationen hoch oder die Signalstärke gering ist. Darüber hinaus konzentrierten sich bestehende Modellierungsansätze auf Voxel-Ebene (z. B. Kang et al., 2012; Castruccio et al., 2018) primär auf aufgabenbasierte Aktivierung oder Gruppeninferenz und versäumten es oft, zwischen intra-regionaler und inter-regionaler Variabilität klar zu unterscheiden oder eine robuste Schätzung individueller Netzwerke bereitzustellen.

Methodik
Die Autoren schlagen ein rigoroses statistisches Framework vor, das auf linearen gemischten Modellen (Linear Mixed-Effects Models) basiert, um die inter-regionale funktionelle Konnektivität direkt aus BOLD-Signalen auf Voxel-Ebene zu schätzen.

1. Spatiotemporales Mixed-Effects-Modell:
   Das beobachtete BOLD-Signal am Voxel $v_{jl}$ in Region $R_j$ und Wavelet-Koeffizient $m$ wird modelliert als:
   $$X_{jlm} = \mu_j + \eta_{jm} + \gamma_{jlm} + \epsilon_{jlm}$$

• $\mu_j$: Regionales Mittelwertniveau.
• $\eta_{jm}$: Ein Mittelwert-Null-Random-Effect, der die regionale Signalkomponente repräsentiert und durch die Korrelationsmatrix $R = \{\rho_{jj'}\}$ inter-regionale Abhängigkeiten induziert.
• $\gamma_{jlm}$: Ein lokaler Nullwert-Random-Effect, der die intra-regionale räumliche Abhängigkeit erfasst und unkorreliert mit $\eta_{jm}$ ist.
• $\epsilon_{jlm}$: Unkorreliertes Messrauschen.
  Das Modell operiert im Wavelet-Bereich, um skalen-spezifische zeitliche Abhängigkeiten zu isolieren und hochfrequentes Rauschen zu reduzieren. Die Kovarianzstruktur nimmt eine Separabilität zwischen den räumlichen und Wavelet-Domänen an und verwendet Matérn- sowie Gauß-Kernel.

2. Zweistufiges Schätzverfahren:
   Um die rechnerische Unhandlichkeit der vollständigen Maximum-Likelihood-Schätzung (MLE) für hochdimensionale spatiotemporale Daten zu adressieren (wobei $N$ die Anzahl der Voxel übersteigen kann, z. B. > 20.000), wurde ein zweistufiger Ansatz entwickelt:

• Stufe 1: Die eingeschränkte Maximum-Likelihood-Methode (Restricted Maximum Likelihood, ReML) wird unabhängig auf die Daten jeder Region angewendet, um die intra-regionalen Kovarianzparameter (räumliche Korrelation, Rauschvarianz) zu schätzen. Dieser Schritt behandelt die regionalen Random-Effects als feste Effekte, um die lokale Struktur zu isolieren.
• Stufe 2: Die inter-regionalen Korrelationsparameter ($\rho_{jj'}$) werden durch Konditionierung auf die Schätzungen aus Stufe 1 geschätzt. Für Paare von Regionen verwenden die Autoren eine Vecchia-Likelihood-Approximation (Vecchia, 1988; Guinness, 2018). Diese Methode approximiert die gemeinsame Likelihood als Produkt bedingter Verteilungen basierend auf einer spärlichen Menge nächstgelegener Nachbarn, wodurch die rechnerische Komplexität von $O(N^3)$ auf ein handhabbares Niveau reduziert wird, während die Schätzgenauigkeit erhalten bleibt.

3. Inferenz:
   Die asymptotischen Eigenschaften der Schätzer werden unter einem Regime expandierender Zeitfenster (divergierende Anzahl von Wavelet-Koeffizienten) etabliert. Approximative Konfidenzintervalle werden mithend der aus der Vecchia-Approximation abgeleiteten Fisher-Informationsmatrix konstruiert.

Zentrale Beiträge

• Neuartiger Schätzer: Die Arbeit führt einen auf Maximum Likelihood basierenden Schätzer für die inter-regionale Konnektivität ein, der aus einem Mixed-Effects-Modell abgeleitet ist, und bewegt sich damit über den momentenbasierten CA-Schätzer hinaus.
• Explizite Modellierung von Variabilität: Das Framework trennt und modelliert explizit die Quellen der Variabilität: inter-regionale Korrelation, intra-regionale räumliche Abhängigkeit und Messrauschen. Dies ermöglicht eine unverzerrte Schätzung der Konnektivitätsparameter, die nicht durch Artefakte der regionalen Signalaggregation konfundiert sind.
• Rechentechnische Skalierbarkeit: Diese Arbeit stellt die erste Anwendung der Vecchia-Likelihood-Approximation auf die Schätzung der funktionellen Konnektivität mit fMRT-Daten auf Voxel-Ebene dar, was die Analyse massiver Datensätze ermöglicht, die zuvor für vollständige Likelihood-Methoden rechnerisch prohibitiv waren.
• Verbesserte Reproduzierbarkeit: Die Methode ist darauf ausgelegt, konsistentere individuelle Gehirnnetzwerke zu erzeugen, und adressiert damit die Herausforderung der Test-Retest-Reliabilität in der fMRT.

Ergebnisse

• Simulationsstudien: Umfangreiche Simulationen, die den vorgeschlagenen ReML- und Vecchia-Schätzer mit CA, dem „Average of Covariances“ (ACE) und EBLUE-Schätzern vergleichen, zeigen eine überlegene Performance. Die vorgeschlagenen Schätzer weisen eine signifikant geringere Verzerrung (Bias) und einen geringeren mittleren quadratischen Fehler auf, insbesondere in Szenarien mit hoher intra-regionaler räumlicher Korrelation oder niedrigem Signal-Rausch-Verhältnis. Entscheidend ist, dass die vorgeschlagene Methode eine nominale Abdeckung von 95 % Konfidenzintervallen erreicht, während CA-basierte Intervalle systematisch scheitern, den wahren Parameter abzudecken.
• Anwendung auf den Human Connectome Project (HCP): Die Methode wurde auf Test-Retest-Daten von 42 jungen Erwachsenen angewendet. Die Konnektivitätsnetzwerke wurden sowohl für die vorgeschlagene Methode als auch für den CA-Schätzer konstruiert. Unter Verwendung des Concordance Correlation Coefficient (CCC) zur Messung der Übereinstimmung zwischen wiederholten Scans (Exam 0 vs. Exam 1) lieferte der vorgeschlagene Mixed-Model-Ansatz konsistent eine höhere Konkordanz als CA. Dies deutde darauf hin, dass die vorgeschlagene Methode zuverlässigere und reproduzierbarere individuelle Gehirnnetzwerke erzeugt.

Bedeutung
Das Paper argumentiert, dass die weit verbreitete Verwendung des CA-Schätzers eine systematische Verzerrung in Studien zur funktionellen Konnektivität einführt, was potenziell die Validität von Gruppenvergleichen und der individuellen Charakterisierung beeinträchtigt. Durch die Bereitstellung eines statistisch fundierten Mixed-Model-Ansatzes, der Voxel-Ebene-Abhängigkeiten und Rauschen berücksichtigt, bietet die Arbeit eine Methode, die genauere und zuverlässigere Schätzungen der Gehirnnetzwerkarchitektur liefert. Die demonstrierte Verbesserung der Test-Retest-Reproduzierbarkeit legt nahe, dass dieser Ansatz zu einem robusteren wissenschaftlichen Verständnis der individuellen Variabilität in Konnektivitätsnetzwerken und deren Assoziationen mit Gesundheitsoutcomes führen kann. Die Autoren merken an, dass die Methode zwar für die fMRT konzipiert ist, das zugrunde liegende Modell jedoch auch für andere Kontexte mit räumlich gruppierten Zeitreihen anwendbar sein könnte.