Block-Term Decomposition Approach to Blind Multi-trial Functional Ultrasound Unmixing

Resumen Técnico: Un enfoque de Descomposición de Bloques de Términos para el desmezclado ciego de ultrasonido funcional multiensayos

Planteamiento del Problema
El ultrasonido funcional (fUS) es una modalidad de neuroimagen de alta resolución que mide los cambios en el volumen sanguíneo cerebral para inferir la actividad neuronal mediante el acoplamiento neurovascular (NVC). Sin embargo, las señales de fUS proporcionan únicamente observaciones indirectas, lo que requiere la separación ciega de las fuentes neuronales latentes y la deconvolución de sus funciones de respuesta hemodinámica (HRF). Los métodos de desmezclado ciego existentes suelen enfrentar desafíos al estimar conjuntamente los mapas espaciales, las señales de activación temporal y las funciones de respuesta hemodinámica (HRF), particularmente cuando se trata de datos de múltiples ensayos donde el promedio de los ensayos puede simplificar excesivamente las respuestas fisiológicas. Los enfoques previos, como los métodos de dos pasos, corren el riesgo de propagación de errores desde la estimación de la firma temporal hacia la etapa de deconvolución. Además, aunque algunos modelos tienen en cuenta la variabilidad de los ensayos, a menudo carecen de restricciones fisiológicas sobre la forma de la HRF o no logran imponer la estructura de bajo rango inherente a los mapas espaciales en la neuroimagen.

Metodología
Los autores proponen un modelo de descomposición de bloques de términos convolutiva (convBTD) basado en datos, diseñado para mediciones de fUS de múltiples ensayos. La innovación central radica en modelar la NVC como un sistema lineal e invariante en el tiempo (LTI) donde la señal medida es la convolución de una señal de activación neuronal de constante por tramos con una HRF dependiente del ensayo y de la fuente.

Características clave del modelado:

1. Formulación tensorial: Los datos de múltiples ensayos se organizan en un tensor de cuarto orden, reestructurado en un tensor de tercer orden $\mathbf{X} \in \mathbb{R}^{N_s \times N_t \times K}$ (donde $N_s$ es el número de píxeles, $N_t$ es el número de muestras temporales y $K$ es el número de ensayos).
2. Estructura espacial de bajo rango: A diferencia del convCPD estándar que vectoriza los mapas espaciales, este modelo impone una estructura rank-$(L_r, L_r, 1)$ en los mapas espaciales ($\mathbf{U}_r$). Esto permite que el mapa espacial de cada fuente se represente como una suma de un pequeño número de productos de Kronecker, consistente con el aprendizaje de diccionarios probabilísticos de atlas multisujeto y estudios previos de fUS.
3. Restricciones fisiológicas:
   • Mapas espaciales: Restringidos a ser no negativos y de bajo rango.
   • Señales de activación: Se asumen constantes por tramos, aplicándose una regularización de Variación Total (TV).
   • HRFs: Modeladas utilizando una forma de función gamma fisiológicamente plausible (Ecuación 5), lo que permite variaciones dependientes del ensayo y de la fuente en lugar de una HRF canónica fija.
4. Marco de optimización: El problema de estimación conjunta se plantea como una tarea de optimización restringida que minimiza un término de fidelidad de datos y un término de regularización TV. La solución emplea un algoritmo de Descenso de Coordenadas por Bloques (BCD) donde las variables ($\mathbf{U}, \mathbf{Z}, \mathbf{H}$) se actualizan alternadamente utilizando Descenso de Gradiente Proyectado (PGD).
   • Las actualizaciones para $\mathbf{U}$ implican proyecciones alternas sobre el ortante no negativo y el conjunto de matrices de bajo rango (vía SVD truncada).
   • Las actualizaciones para $\mathbf{H}$ implican la proyección sobre la variedad de la HRF mediante el ajuste de curvas de Mínimos Cuadrados No Lineales (NLS).
   • Las actualizaciones para $\mathbf{Z}$ utilizan un operador proximal para manejar el término TV no suave.

Contribuciones Clave
El artículo identifica dos contribuciones primarias:

1. Extensión de trabajos previos: Los autores extienden el modelo BTD de rank-$(L_r, L_r, 1)$ de [10] a un escenario más realista que involucra HRFs dependientes del ensayo y/o de la fuente. Crucialmente, fusionan los pasos previamente separados de separación de fuentes y deconvolución en un único marco de optimización conjunta para mitigar la propagación de errores.
2. Mejoras estructurales y fisiológicas sobre convCPD: En comparación con los modelos de CPD invariante por desplazamiento (shiftCPD) y CPD convolutiva (convCPD) (por ejemplo, [25]), este trabajo promueve explícitamente la estructura de bajo rango de los mapas espaciales en lugar de vectorizarlos. Además, adjunta una interpretación fisiológica específica a los filtros dependientes del ensayo al restringirlos a una forma de HRF de función gamma, en lugar de tratarlos como respuestas de impulso genéricas y dispersas.

Resultados de la Simulación
Los autores evaluaron su método utilizando datos sintéticos de fUS de múltiples ensayos con $R=2$ fuentes, $K=4$ ensayos y variaciones de la relación señal-ruido (SNR) de -5 dB a 20 dB.

• Mapas espaciales: El método demostró una alta precisión en la recuperación de los mapas espaciales en todos los niveles de SNR, con errores relativos medianos alrededor de 0.1 en -5 dB, disminuyendo a medida que la SNR aumentaba.
• Señales de activación: Los perfiles de activación temporal se estimaron de manera fiable, aunque con errores medianos ligeramente superiores (~0.2) en comparación con los mapas espaciales.
• Estimación de la HRF: Los autores señalan que la estimación de la HRF sigue siendo el aspecto más desafiante del problema. Si bien los errores medianos fueron comparables a los de otros componentes (~0.2), la variabilidad (percentil 75 y valores atípicos) fue significativamente mayor. Esto se atribuye a la ambigüedad inherente de la deconvolución, donde pequeños cambios en un factor (por ejemplo, el tiempo de activación) pueden compensarse con cambios opuestos en el otro (la forma de la HRF) manteniendo un buen ajuste de los datos.
• Robustez: El método mostró robustez a la inicialización; el rendimiento utilizando la inicialización que minimizó el error de reconstrucción de datos (un escenario práctico sin verdad de campo) fue solo marginalmente peor que el uso de la inicialización óptima de tipo oráculo.

Significancia y Reivindicaciones
El artículo posiciona este trabajo como un primer paso en un proyecto de investigación destinado a desarrollar una metodología completa basada en tensores para la descomposición ciega de fUS. Los autores afirman modestamente que el marco propuesto demuestra con éxito:

• La viabilidad de estimar conjuntamente mapas espaciales, señales de activación y HRFs dependientes de los ensayos de manera dirigida por datos.
• La capacidad de recuperar mapas espaciales y perfiles de activación con una precisión razonable a través de varios niveles de ruido.
• La confirmación de que, mientras que la recuperación espacial y temporal es robusta, la estimación de la HRF sigue siendo el cuello de botella debido a la ambigüedad de la convolución.

Los autores declaran explícitamente que el enfoque actual asume que el número de fuentes ($R$) y sus rangos espaciales ($L_r$) se conocen a priori, una limitación que pretenden abordar en trabajos futuros mediante diagnósticos de consistencia central o inferencia bayesiana. También señalan que, si bien el estudio actual se centra en datos de un solo sujeto, el marco podría extenderse a entornos de múltiples sujetos con la correspondiente registro anatómico.