⚡ electrical engineering

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

This paper proposes a constrained optimization framework based on convolutive block-term tensor decomposition to blind-unmix multi-trial functional ultrasound data, enabling the simultaneous recovery of low-rank spatial maps, neuronal activation signals, and trial-specific hemodynamic response functions.

Original authors: Sofia-Eirini Kotti, Eleftherios Kofidis, Borbála Hunyadi

Published 2026-06-09

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Sofia-Eirini Kotti, Eleftherios Kofidis, Borbála Hunyadi

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine your brain is a bustling city, and Functional Ultrasound (fUS) is a high-tech drone flying overhead, trying to take a video of the city's activity. However, there's a catch: the drone doesn't see the people (neurons) directly. Instead, it sees the traffic flow (blood volume) that changes whenever people start moving. This is called "neurovascular coupling."

The problem is that the traffic video is a messy mix. It shows:

The actual events (people moving because of a specific stimulus).
Random background noise (cars moving for no reason, weather effects, or camera glitches).
A time delay (it takes a few seconds for the traffic to build up after people start moving).

The paper presents a new mathematical "recipe" to untangle this messy video and figure out exactly who was doing what, when, and how fast the traffic responded, without needing to know the answers beforehand (this is called "blind unmixing").

Here is how their approach works, broken down into simple concepts:

1. The "Multi-Trial" Puzzle

Usually, scientists might take one long video and try to guess what happened. But this paper suggests looking at multiple trials (like watching the same movie scene repeated 4 times).

The Analogy: Imagine you are trying to hear a specific instrument in a band, but the sound is muffled. If you listen to the same song played four times, the instrument's melody stays the same, but the background noise changes slightly each time. By comparing all four recordings, you can isolate the instrument much better.
The Innovation: The authors realized that the "traffic delay" (the Hemodynamic Response Function, or HRF) might not be exactly the same every time the brain reacts. Sometimes the traffic builds up faster; sometimes slower. Their model allows this "delay" to change slightly from trial to trial, rather than forcing it to be identical.

2. The "Block-Term" Recipe (The Core Idea)

To solve this, they use a mathematical tool called Block-Term Decomposition (BTD).

The Analogy: Think of the messy video as a giant jigsaw puzzle where the pieces are mixed up.
- Standard puzzles assume every piece is a single, unique shape (Rank-1).
- This paper's approach realizes that some "pieces" (brain regions) are actually made of a small cluster of connected shapes (Low-Rank).
- They treat the data as a 3D block: Space (where in the brain), Time (when it happened), and Trials (which repetition).
- They assume that while the location of the activity is consistent (the same brain area lights up), the timing of the blood flow response might wiggle a bit between trials.

3. The Three Things They Are Trying to Find

The algorithm tries to separate the messy data into three distinct ingredients:

The Map (Where): A low-resolution picture showing which parts of the brain are active. The authors assume these maps are "low-rank," meaning they are smooth and simple shapes, not chaotic noise.
The Signal (When): A timeline showing when the neurons fired. They assume this signal is "piecewise constant," meaning it's like a light switch that is either ON or OFF for a block of time, rather than flickering randomly.
The Delay (How Fast): The shape of the blood flow response (HRF). They use a specific mathematical curve (based on biology) to describe how blood rushes to the brain, ensuring the answer makes physical sense.

4. How They Solve It (The Optimization)

Solving this is like trying to balance a stack of cards while the wind is blowing. You can't just guess; you have to adjust step-by-step.

The Method: They use a "Projected Gradient Descent."
The Analogy: Imagine you are trying to find the lowest point in a foggy valley (the best solution). You take a step down. But, you have rules:
- You can't go into the "negative" zone (blood flow can't be negative).
- You must stay on a specific path (the shape must look like a real brain map).
- The curve must look like a real biological response.
- Every time you take a step, the algorithm checks these rules and "projects" you back onto the valid path if you wandered off. They repeat this until the picture becomes clear.

5. What They Found (The Results)

They tested their recipe using computer simulations (fake data with known answers) at different noise levels.

The Good News: They were very good at finding the Map (where the activity is) and the Signal (when it happened), even when the data was very noisy.
The Challenge: Finding the exact shape of the Delay (HRF) was the hardest part.
- Why? It's like trying to guess the exact shape of a filter by only looking at the filtered sound. If the sound is slightly different, you might guess the filter is slightly different, even if the original sound was the same. This is called "ambiguity."
- However, even though the exact shape was tricky, the model still captured the general trend and relative differences between trials.

Summary

This paper introduces a smarter way to clean up brain ultrasound videos. Instead of assuming the brain's blood response is a rigid, unchanging machine, they treat it as a flexible system that changes slightly every time. By looking at multiple repetitions of an experiment and using a special mathematical "lens" (Block-Term Decomposition), they can separate the brain's activity map, the timing of the event, and the blood flow delay, even when the data is noisy.

The Bottom Line: They successfully built a tool that can "unmix" complex brain data, recovering clear maps and timing signals, though pinning down the exact speed of the blood flow response remains the most difficult part of the puzzle.

Technical Summary: A Block-Term Decomposition Approach to Blind Multi-trial Functional Ultrasound Unmixing

Problem Statement
Functional ultrasound (fUS) is a high-resolution neuroimaging modality that measures cerebral blood volume changes to infer neuronal activity via neurovascular coupling (NVC). However, fUS signals provide only indirect observations, requiring the blind separation of latent neuronal sources and the deconvolution of their hemodynamic responses. Existing blind unmixing methods often face challenges in jointly estimating spatial maps, temporal activation signals, and hemodynamic response functions (HRFs), particularly when dealing with multi-trial data where trial-averaging may oversimplify physiological responses. Previous approaches, such as two-step methods, risk error propagation from temporal signature estimation to the deconvolution stage. Furthermore, while some models account for trial variability, they often lack physiological constraints on the HRF shape or fail to enforce the low-rank structure inherent to spatial maps in neuroimaging.

Methodology
The authors propose a data-driven convolutive Block-Term Decomposition (convBTD) model tailored for multi-trial fUS measurements. The core innovation lies in modeling the NVC as a linear time-invariant (LTI) system where the measured signal is the convolution of a piecewise-constant neuronal activation signal with a trial- and source-dependent HRF.

Key modeling features include:

Tensor Formulation: Multi-trial data is organized into a fourth-order tensor, reshaped into a third-order tensor $\mathbf{X} \in \mathbb{R}^{N_s \times N_t \times K}$ (where $N_s$ is the number of pixels, $N_t$ time samples, and $K$ trials).
Low-Rank Spatial Structure: Unlike standard Convolutional CPD (convCPD) which vectorizes spatial maps, this model enforces a rank- $(L_r, L_r, 1)$ structure on the spatial maps ( $\mathbf{U}_r$ ). This allows each source's spatial map to be represented as a sum of a small number of Kronecker products, consistent with the multisubject dictionary learning probabilistic atlas and previous fUS studies.
Physiological Constraints:
- Spatial Maps: Constrained to be non-negative and low-rank.
- Activation Signals: Assumed to be piecewise-constant, enforced via Total Variation (TV) regularization.
- HRFs: Modeled using a physiologically plausible gamma function shape (Equation 5), allowing for trial- and source-dependent variations rather than a fixed canonical HRF.
Optimization Framework: The joint estimation problem is cast as a constrained optimization task minimizing a data fidelity term and a TV regularization term. The solution employs a Block Coordinate Descent (BCD) algorithm where variables ( $\mathbf{U}, \mathbf{Z}, \mathbf{H}$ $U, Z, H$ ) are updated alternately using Projected Gradient Descent (PGD).
- Updates for $\mathbf{U}$ involve alternating projections onto the non-negative orthant and the set of low-rank matrices (via truncated SVD).
- Updates for $\mathbf{H}$ involve projecting onto the HRF manifold using Non-Linear Least Squares (NLS) curve fitting.
- Updates for $\mathbf{Z}$ utilize a proximal operator to handle the non-smooth TV term.

Key Contributions
The paper identifies two primary contributions:

Extension of Previous Work: The authors extend the rank- $(L_r, L_r, 1)$ BTD model from [10] to a more realistic scenario involving trial- and/or source-dependent HRFs. Crucially, they merge the previously separate steps of source separation and deconvolution into a single, joint optimization framework to mitigate error propagation.
Structural and Physiological Enhancements over convCPD: Compared to the shift-invariant CPD (shiftCPD) and convolutional CPD (convCPD) models (e.g., [25]), this work explicitly promotes the low-rank structure of spatial maps rather than vectorizing them. Additionally, it attaches a specific physiological interpretation to the trial-dependent filters by constraining them to a gamma-function HRF shape, rather than treating them as generic sparse impulse responses.

Simulation Results
The authors evaluated the method using synthetic multi-trial fUS data with $R=2$ sources, $K=4$ trials, and varying Signal-to-Noise Ratios (SNR) from -5 dB to 20 dB.

Spatial Maps: The method demonstrated high accuracy in recovering spatial maps across all SNR levels, with median relative errors around 0.1 at -5 dB and decreasing as SNR increased.
Activation Signals: Temporal activation profiles were reliably estimated, though with slightly higher median errors (~0.2) compared to spatial maps.
HRF Estimation: The authors note that HRF estimation remains the most challenging aspect of the problem. While median errors were comparable to other components (~0.2), the variability (75th percentile and outliers) was significantly higher. This is attributed to the inherent ambiguity in deconvolution, where small changes in one factor (e.g., activation timing) can be compensated by opposite changes in the other (HRF shape) while maintaining a good data fit.
Robustness: The method showed robustness to initialization, with performance using the initialization that minimized data reconstruction error (a practical scenario without ground truth) being only marginally worse than using the oracle-best initialization.

Significance and Claims
The paper positions this work as a first step in a research project aimed at developing a complete tensor-based methodology for blind fUS decomposition. The authors modestly claim that the proposed framework successfully demonstrates:

The feasibility of jointly estimating spatial maps, activation signals, and trial-dependent HRFs in a data-driven manner.
The ability to recover spatial maps and activation profiles with reasonable accuracy across various noise levels.
The confirmation that while spatial and temporal recovery is robust, HRF estimation remains the bottleneck due to convolutional ambiguity.

The authors explicitly state that the current approach assumes the number of sources ( $R$ ) and their spatial ranks ( $L_r$ ) are known a priori, a limitation they intend to address in future work through core consistency diagnostics or Bayesian inference. They also note that while the current study focuses on single-subject data, the framework could be extended to multi-subject settings with appropriate anatomical co-registration.