Quantile Function-Based Models for Neuroimaging Classification Using Wasserstein Regression

TL;DR

This paper introduces a novel quantile function-based Wasserstein regression method for neuroimaging classification, significantly improving diagnostic accuracy for mTBI by utilizing full distributional information of brain signals.

Contribution

It develops a Wasserstein-Fréchet regression framework that preserves complete distributional characteristics of neuroimaging data for improved classification of mTBI.

Findings

Achieved up to 98% test accuracy on unseen data.

Preserved full distributional information improves diagnostic performance.

Provides a statistically principled approach for neuroimaging classification.

Abstract

We propose a novel quantile function-based approach for neuroimaging classification using Wasserstein-Fr\'echet regression, specifically applied to the detection of mild traumatic brain injury (mTBI) based on the MEG and MRI data. Conventional neuroimaging classification methods for mTBI detection typically extract summary statistics from brain signals across the different epochs, which may result in the loss of important distributional information, such as variance, skewness, kurtosis, etc. Our approach treats complete probability density functions of epoch space results as functional response variables within a Wasserstein-Fr\'echet regression framework, thereby preserving the full distributional characteristics of epoch results from $L_{1}$ minimum norm solutions. The global Wasserstein-Fr\'echet regression model incorporating covariates (age and gender) allows us to directly compare the distributional patterns between healthy control subjects and mTBI patients. The classification procedure computes Wasserstein distances between estimated quantile functions from control and patient groups, respectively. These distances are then used as the basis for diagnostic decisions. This framework offers a statistically principled approach to improving diagnostic accuracy in mTBI detection. In practical applications, the test accuracy on unseen data from Innovision IP's dataset achieves up to 98\%.

Tables

Table 1. Table 1: Prediction accuracy comparison for mild traumatic brain injury (mTBI) detection. Baselines 1 and 2 correspond to the linear and SVM classifiers based on average spectral magnitudes, respectively. We only report the accuracies for all frequency bands.

References	Accuracy	Based-on	Method
Huang et al. (2014)	0.83	Resting-state MEG, Age, Gender	Statistical Model (VESTAL)
Vakorin et al. (2016)	0.88	Resting-state MEG, Connectivity	Machine Learning
Vergara et al. (2017)	0.841	Resting-state MEG, Connectivity	Support Vector Machine (SVM)
		and Fractional Anisotropy
McNerney et al. (2019)	0.91	EEG+Self-reported Symptoms	Machine Learning
Thanjavur et al. (2021)	0.926	EEG	Recurrent Neural Network (RNN)
Itälinna et al. (2023)	0.79	Resting-state MEG, Source Magnitude	Support Vector Machine (SVM)
Baseline 1	0.9388	Resting-state MEG, Source Magnitude	OLS Classifiers
Baseline 2	0.9337	Resting-state MEG, Source Magnitude	SVM Classifiers
Our method	0.9802	Resting-state MEG, Source Magnitude	Wasserstein-Fréchet Regression

Equations

$$d_{\mathrm{W}}^2(f,g):=\underset{M\in {\mathcal{M}}_{f,g}}{\inf }{\int }_{\mathbb{R}}{\left(M(u)-u\right)}^{2}f(u)\mathrm{d}u.$$

$$d_{\mathrm{W}}^2(f,g)={\int }_{\mathbb{R}}{\left(M_{f,g}^{\text{opt}}(u)-u\right)}^{2}f(u)\mathrm{d}u={\int }_0^1{\left(F^{-1}(t)-G^{-1}(t)\right)}^2\mathrm{d}t,$$

$$f_{\oplus }^{\ast }:=\underset{f\in \mathcal{D}}{\mathrm{argmin}}E\left(d_{\mathrm{W}}^2(\mathfrak{F},f)\right),{\mathrm{Var}}_{\oplus }(\mathfrak{F}):=E\left(d_{\mathrm{W}}^2(\mathfrak{F},f_{\oplus }^{\ast })\right).$$

$$f_{\oplus }(\mathbf{z}):=\underset{f\in \mathcal{D}}{\mathrm{argmin}}E\left[d_{\mathrm{W}}^2(\mathfrak{F},f)|Z=\mathbf{z}\right],\mathbf{z}\in {\mathfrak{S}}_Z.$$

$$\begin{array}{rl}{\mathrm{Var}}_{\oplus }(\mathfrak{F}|Z=\mathbf{z})& =E\left[d_{\mathrm{W}}^2(\mathfrak{F},f_{\oplus }(\mathbf{z}))|Z=\mathbf{z}\right]={\int }_{\mathbb{R}}E\left[(T_{\mathbf{z}}(u)-u{)}^2|Z=\mathbf{z}\right]f_{\oplus }(\mathbf{z},u)\mathrm{d}u\\ & ={\int }_{\mathbb{R}}\mathrm{Var}(T_{\mathbf{z}}(u)|Z=\mathbf{z})f_{\oplus }(\mathbf{z},u)\mathrm{d}u.\end{array}$$

$$f_{\oplus }(\mathbf{z})=\underset{f\in \mathcal{D}}{\mathrm{argmin}}E\left[s(Z,\mathbf{z})d_{\mathrm{W}}^2(\mathfrak{F},f)\right],$$

$${\hat{Q}}_{\oplus }(\mathbf{z})=\underset{Q\in \mathfrak{Q}}{\mathrm{argmin}}\sum _{i=1}^n{s}_{in}(\mathbf{z}){\Vert Q-Q_i\Vert }_{L^2}^2.$$

$$d_{1}\leftarrow\Big{(}\!\int_{0}^{1}\!\big{[}\hat{Q}(\mathbf{z}_{0})-\hat{Q}_{\oplus}^{\mathrm{ctl}}(\mathbf{z}_{0})\big{]}^{2}\mathrm{d}t\Big{)}^{1/2},\quad d_{2}\leftarrow\Big{(}\!\int_{0}^{1}\!\big{[}\hat{Q}(\mathbf{z}_{0})-\hat{Q}_{\oplus}^{\mathrm{mTBI}}(\mathbf{z}_{0})\big{]}^{2}\mathrm{d}t\Big{)}^{1/2}.$$

$${\hat{y}}_0\leftarrow \{\begin{array}{ll}\mathrm{control},& \text{if }{d}_1\le k{d}_2,\\ \mathrm{mTBI},& \text{otherwise}.\end{array}$$

$$\hat{f}(x;h)=\frac{1}{nh}\sum _{i=1}^{n}K\left(\frac{x-x_i}{h}\right),$$

$$\frac{\partial }{\partial h}f(x;h)=Lf(x;h),x\in \mathcal{X},h>0,$$

$${}_{\ast }t={\left(\frac{1}{2n\sqrt{\pi }{\Vert f^{\prime \prime }\Vert }^2}\right)}^{2/5}.$$

$${\Vert f^{(j)}\Vert }^2=\frac{(-1{)}^j}{n^2}\sum _{k=1}^n\sum _{m=1}^n{\varphi }^{(2j)}(x_k,x_m;2t_j),$$

$${\Vert f^{(j)}\Vert }^2=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }{\left|\mathcal{F}[{\hat{f}}^{(j)}](\omega ,t)\right|}^2\mathrm{d}\omega =\frac{1}{2\pi }{\int }_{-\infty }^{\infty }{\omega }^{2j}{e}^{-{\omega }^{2}t}{\left|{\hat{f}}_{\text{emp}}(\omega )\right|}^2\mathrm{d}\omega .$$

$${\Vert f^{(j)}\Vert }^2\approx 2{\pi }^{2j}\sum _k{I}_k^j\ast a_k\ast e^{-I_k{\pi }^2{t}_j}.$$

$${}_{\ast }{t}_j={\left(\frac{1+1/2^{j+1/2}}{3}\frac{1\times 3\times 5\times \cdots \times (2j-1)}{n\sqrt{\pi /2}{\Vert f^{(j+1)}\Vert }^2}\right)}^{2/(3+2j)}.$$

$${}_{\ast }{t}_{cdf}={\left(\frac{1}{\sqrt{\pi }n{\Vert f^{\prime }\Vert }^2}\right)}^{2/3}.$$

$$\underset{{\mathbf{q}}^{\ast }}{\min }\sum _{j=1}^m{\left[t_j\left(q_j-q_j^{\ast }\right)\right]}^2\text{subject to}{q}_{j+1}^{\ast }-q_j^{\ast }\ge 0,j=1,\dots ,m-1,$$

$$f(x)=0.3\frac{1}{\sqrt{2\pi }}{e}^{-\frac{1}{2}(x-6{)}^2}+0.6\frac{1}{0.7\sqrt{2\pi }}{e}^{-\frac{1}{2}{\left(\frac{x-9.5}{0.7}\right)}^2}+0.1\frac{1}{0.5\sqrt{2\pi }}{e}^{-\frac{1}{2}{\left(\frac{x-12}{0.5}\right)}^2}.$$

$$\mu ={\beta }_0+{\beta }_1\cdot \text{age}+{\beta }_2\cdot \text{gender}$$

$${\bar{y}}_i={\alpha }_0+{\alpha }_1\cdot {\text{age}}_i+{\alpha }_2\cdot {\text{gender}}_i+{ϵ}_i,i=1,2,\dots ,m,$$

Full text

Keywords— Wasserstein distance, quantile function, mild traumatic brain injury, regression

1 Introduction

Mild traumatic brain injury (mTBI) affects millions of people annually and can lead to persistent cognitive and affective symptoms despite often unremarkable findings on conventional structural neuroimaging (Lewine et al., 1999, 2007; Huang et al., 2014; Allen et al., 2021). This diagnostic gap has driven interest in functional neurophysiology and network connectivity approaches that are sensitive to subtle post-injury alterations (Dunkley et al., 2015; Vakorin et al., 2016; Antonakakis et al., 2017; Vergara et al., 2017; Wang et al., 2017; Allen et al., 2021). Magnetoencephalography (MEG) and electroencephalography (EEG) offer millisecond-level insights into brain activity. These techniques have uncovered consistent patterns of abnormality in mTBI patients. Research shows increased slow-wave activity in these patients compared to healthy controls. Changes in resting-state network connectivity have also been documented (Huang et al., 2014; Vakorin et al., 2016; Vergara et al., 2017). These researches motivate quantitative biomarkers that can generalise across subjects and sites (or systems) while remaining transparent and clinically interpretable (Gross et al., 2013; Niso et al., 2019; Boon et al., 2021; Ferrante et al., 2022).

A substantial literature applies machine learning (ML) to neurophysiology and multimodal data for mTBI detection (Vakorin et al., 2016; Itälinna et al., 2023; Vergara et al., 2017; McNerney et al., 2019; Thanjavur et al., 2021; Vivaldi et al., 2021; Caiola et al., 2023; Wall et al., 2022). Examples include normative modelling and support vector machines (SVMs) on MEG source features (Itälinna et al., 2023), recurrent neural networks on raw resting-state EEG (Thanjavur et al., 2021), deep audio classifiers for voice-based screening (Wall et al., 2022), and combined frontopolar EEG with self-reported symptoms (McNerney et al., 2019). Beyond electrophysiology, resting-state functional network connectivity and diffusion MRI (fractional anisotropy) have achieved encouraging performance with classical ML (Vergara et al., 2017). Reported accuracies in these studies typically range between about 0.79 and 0.93 (Itälinna et al., 2023; Vergara et al., 2017; McNerney et al., 2019; Thanjavur et al., 2021; Wall et al., 2022), also see Table˜1. However, these advances in black-box AI approaches encounter significant obstacles when it comes to clinical implementation and medico-legal acceptance (Rudin, 2019; Leslie, 2020; World Health Organization, 2021). In particular, under current UK practice, automated, non-transparent models are unlikely to be accepted as standalone evidence in court; expert opinion requires methods that are explainable, auditable, and grounded in established statistical principles (Law Commission of England and Wales, 2011). At the same time, purely statistical models that focus on interpretability often sacrifice predictive accuracy, say 83% in Huang et al. (2014) and 79% in Itälinna et al. (2023). This challenge highlights the importance of developing statistical methods that balance interpretability with strong predictive performance. (Breiman, 2001; Barredo Arrieta et al., 2020; Bollaerts et al., 2006).

We bridge this gap by proposing an interpretable framework based on Wasserstein–Fréchet regression (Petersen et al., 2021) of quantile functions derived from source-space power of MEG signals. Rather than summarising the epoch magnitudes with low-dimensional moments (Li et al., 2025), we treat the entire empirical distribution of source magnitudes within a region and frequency band as the response variable. Quantile functions preserve full distributional information and admit a natural geometry via the Wasserstein metric; conditional Fréchet means under this metric can be expressed as weighted averages of quantile functions, yielding estimators that are easy to compute and explain (Petersen et al., 2021; Petersen and Müller, 2016; Petersen et al., 2022; Ramsay and Silverman, 2006; Wang et al., 2016). In our pipeline, epoch-wise $L_{1}$-constrained source magnitude imaging (VESTAL) produces robust regional magnitude distributions (Huang et al., 2014). Diffusion-based kernel density estimation (Pelz et al., 2023) with fixed-point bandwidth selection provides stable density and quantile estimates. Global Wasserstein regression then links these distributional responses to covariates (age, gender), producing group- and covariate-conditioned “prototype” distributions. For an individual, classification reduces to comparing Wasserstein distances between the subject’s quantile function and covariate-matched group prototypes, yielding a transparent decision rule across quantiles and regions.

This design has several implications for accuracy and translation. First, modelling full distributions mitigates information loss from ad hoc summary statistics and captures heteroscedastic and skewed changes commonly observed post-injury (Huang et al., 2014; Allen et al., 2021; Vakorin et al., 2016). Moreover, our nonparametric, distribution-free density estimation approach makes no assumptions on the distribution of the data, allowing the method to adapt flexibly to arbitrary distribution that may arise in neuroimaging applications. Second, the Wasserstein geometry inherently captures mass transport between distributions, which improves sensitivity to changes in slow and gamma frequency bands previously linked to mTBI (Huang et al., 2014; Vakorin et al., 2016). Third, the decision process remains fully transparent: practitioners can examine distances to covariate-matched group means and quantile-level contributions, providing clear advantages for clinical reporting and medico-legal review.

The remainder of the paper is organised as follows. Section 2 formulates the Wasserstein–Fréchet regression on quantile functions and details density/quantile estimation via diffusion-based kernel density estimator (KDE). Section 3 presents simulations comparing the proposed estimator to alternatives and evaluates classification performance relative to standard linear baselines. Section 4 applies the full MRI+MEG pipeline to Cam-CAN controls (Cam-CAN et al., 2014) and a mTBI cohort, including multiband analyses (delta, slow gamma, fast gamma) and cross-validated evaluation. Section 5 discusses limitations, generalisability, and the translational potential of interpretable, distributional MEG biomarkers for mTBI.

2 Methodology

In many applications, probability density functions arise as functional response objects in regression with Euclidean covariates. Petersen et al. (2021) develop a regression framework for density-valued responses under the Wasserstein metric with vector predictors. Their global model requires no local smoothing or other tuning parameters. We summarise the elements relevant to our setting and then detail the Wasserstein–Fréchet regression used here.

Let $f_{i}$ denote the probability density function for subject $i=1,\dots,n$. Each subject $i$ provides a sample $\mathbf{x}_{i}=(x_{i1},\dots,x_{in_{i}})^{\top}$ drawn independently from $f_{i}$. Sample sizes $n_{i}$ may vary between subjects. We also record covariates $\mathbf{z}_{i}=(z_{i1},\dots,z_{ip})^{\top}$ such as age and gender. Each subject has a binary clinical status label $y_{i}\in\{0,1\}$. The complete training data consists of $\{(y_{i},\mathbf{x}_{i},\mathbf{z}_{i})\}_{i=1}^{n}$. For a new subject, we observe $(\mathbf{x}_{0},\mathbf{z}_{0})$ and aim to predict the corresponding $y_{0}$.

Let $\mathcal{D}$ denote the set of univariate probability density functions $f$ on $\mathbb{R}$ that are absolutely continuous and have finite second moment, i.e., $\int_{\mathbb{R}}u^{2}f(u)\,\mathrm{d}u<\infty$. For $f,g\in\mathcal{D}$, let $\mathcal{M}_{f,g}$ be the set of measurable maps $M:\mathbb{R}\to\mathbb{R}$ such that, if $U\sim f$, then $M(U)\sim g$. The squared Wasserstein distance between $f$ and $g$ is defined by

[TABLE]

In one dimension, we know the infimum has a specific solution. The optimal transport map is $M_{f,g}^{\textrm{opt}}=G^{-1}\circ F$. Here, $F$ and $G$ represent the cumulative distribution functions (CDFs) of $f$ and $g$. This leads to a closed-form expression:

[TABLE]

where the second equality follows from the change of variables $t=F(u)$Here, $F^{-1}$ and $G^{-1}$ are the quantile functions for the two distributions. This shows that in one dimension, the Wasserstein distance becomes the $L^{2}$-distance between quantile functions.

2.1 Random densities

A random density $\mathfrak{F}$ is defined as a random variable that takes values almost surely in $\mathcal{D}$ (Petersen et al., 2021). The associated random CDF and quantile function are denoted as $F$ and $Q=F^{-1}$ respectively (Parzen, 1979). For clarity, we use specific notation conventions throughout this work. We use $u,v\in\mathbb{R}$ as arguments for densities and CDFs. In contrast, we use $s,t\in[0,1]$ as arguments for quantile and quantile-density functions. The Wasserstein-Fréchet mean and variance of $\mathfrak{F}$ are defined as follows:

[TABLE]

In regression applications, we often need to model how $\mathfrak{F}$ changes with predictor variables $Z\in\mathbb{R}^{p}$. The pair $(Z,\mathfrak{F})$ follows a joint distribution $\mathcal{G}$ on $\mathbb{R}^{p}\times\mathcal{D}$. The quantities in (2) represent the overall Fréchet mean and variance of $\mathfrak{F}$ without conditioning on $Z$. Let $\mathfrak{S}_{Z}$ be the support of $Z$’s marginal distribution. We want to find the Fréchet regression function. This maps each covariate value to its corresponding conditional Fréchet mean:

[TABLE]

Next, we define the conditional Fréchet variance. Let $F_{\oplus}(\mathbf{z})$ be the CDF that corresponds to $f_{\oplus}(\mathbf{z})$. We denote its quantile function as $Q_{\oplus}(\mathbf{z})$. To clarify our notation: $f_{\oplus}(\mathbf{z},u)$ represents the value of the conditional mean density $f_{\oplus}(\mathbf{z})$ evaluated at point $u\in\mathbb{R}$. Similarly, $F_{\oplus}(\mathbf{z},u)$ and $Q_{\oplus}(\mathbf{z},t)$ denote the corresponding CDF and quantile function values at $u$ and $t\in[0,1]$, respectively.

For a pair $(Z,\mathfrak{F})$, we define an optimal transport map $T_{\mathbf{z}}:\mathbb{R}\to\mathbb{R}$ for each $\mathbf{z}$. This map transports from $f_{\oplus}(\mathbf{z})$ to the random density $\mathfrak{F}$. We express it as $T_{\mathbf{z}}(u)\coloneqq Q\!\left(F_{\oplus}(\mathbf{z},u)\right)$, where $Q$ represents the random quantile function of $\mathfrak{F}$. From equation (1), we know that $E\!\left\{Q(t)\mid Z=\mathbf{z}\right\}=Q_{\oplus}(\mathbf{z},t)$. This relationship tells us that $E\!\left\{T_{\mathbf{z}}(u)\mid Z=\mathbf{z}\right\}=u$ for all $u$ where $f_{\oplus}(\mathbf{z},u)>0$. We can now express the conditional Fréchet variance as

[TABLE]

The approach above requires that both marginal and conditional Wasserstein mean densities exist and have unique solutions. This property does not hold automatically in all cases. However, the sufficient conditions (A1) through (A3) established by Petersen et al. (2021) ensure these requirements are met. We work under these assumptions and move directly to the regression framework.

2.2 Global Wasserstein-Fréchet regression

To examine the effects of covariates $Z$ on conditional Wasserstein means, we use a global regression approach for the conditional means $f_{\oplus}(\mathbf{z})$ from (3). Following Petersen and Müller (2019), this Fréchet regression model uses a weighted Fréchet mean structure:

[TABLE]

where $s(Z,\mathbf{z})=1+(Z-\mu)^{\top}\Sigma^{-1}(\mathbf{z}-\mu),\quad\mu=E(Z),\ \Sigma=\mathrm{Var}(Z)$ is the weight function, and $\Sigma$ is positive definite. This model extends linear regression to work with density-valued responses. Instead of using $Y$ and the usual Euclidean space $(\mathbb{R},|\cdot|)$, we substitute $\mathfrak{F}$ and the Wasserstein space $(\mathcal{D},d_{\mathrm{W}})$. Model (5) gives us a way to specify the conditional Wasserstein mean of $\mathfrak{F}$. In what follows, we concentrate on this mean specification.

2.3 Estimation

To estimate the regression function $f_{\oplus}(\mathbf{z})$, we use an empirical weighted least squares criterion in Wasserstein space as in (5). We calculate the sample mean $\bar{Z}=n^{-1}\sum_{i=1}^{n}Z_{i}$ and the sample covariance matrix $\hat{\Sigma}=n^{-1}\sum_{i=1}^{n}(Z_{i}-\bar{Z})(Z_{i}-\bar{Z})^{\top}$. We then define empirical weights as $s_{in}(\mathbf{z})=1+(Z_{i}-\bar{Z})^{\top}\hat{\Sigma}^{-1}(\mathbf{z}-\bar{Z})$. Let $\mathfrak{Q}$ denote the set of quantile functions in $L^{2}[0,1]$. Using ${\lVert\cdot\rVert}_{L^{2}}$ as the standard Hilbert norm on $L^{2}[0,1]$, we estimate the conditional mean quantile function $Q_{\oplus}(\mathbf{z})$ by solving:

[TABLE]

where $Q_{i}$ is the subject-specific quantile function.

The above introduction of Wasserstein-Fréchet regression are summarised based on Petersen et al. (2021). Implementation details are provided in Algorithm 1 of Petersen et al. (2021), the R package WRI for Wasserstein-Fréchet regression is available online here. In finite samples, $\hat{Q}_{\oplus}(\mathbf{z})$ need not be strictly increasing and thus may fail to induce a density. However, the Lemma 2 in Petersen et al. (2021) guarantees that, with high probability for large samples, $\hat{Q}_{\oplus}(\mathbf{z})$ is strictly increasing.

Having obtained the conditional mean quantile estimator $\hat{Q}_{\oplus}(\mathbf{z})$ from the weighted Fréchet criterion in (6), we fit two separate Wasserstein–Fréchet regression models. We train one model using healthy controls and another using mTBI patients. We restrict the optimization in (6) to subjects within each group $g\in\{\mathrm{ctl},\mathrm{mTBI}\}$. This produces covariate-matched prototype quantiles $\hat{Q}_{\oplus}^{g}(\mathbf{z})$ that capture the distributional characteristics at covariate value $\mathbf{z}$. For a new subject with covariates $\mathbf{z}_{0}$ and estimated quantile function $\hat{Q}(\mathbf{z}_{0})$, we classify by comparing Wasserstein distances. We compute distances between $\hat{Q}(\mathbf{z}_{0})$ and both prototypes: $\hat{Q}_{\oplus}^{\mathrm{ctl}}(\mathbf{z}_{0})$ and $\hat{Q}_{\oplus}^{\mathrm{mTBI}}(\mathbf{z}_{0})$. Algorithm 1 provides the complete procedure and decision rule. We select the classification threshold through cross-validation on the training set.

The theoretical framework of Petersen et al. (2021) assumes that subject-level densities and quantile functions are directly observable. However, real neuroimaging studies must estimate these functions from finite, noisy data samples. This creates several practical challenges. We must ensure smoothness and proper boundary behaviour. Most importantly, we need to guarantee that estimated quantile functions remain monotonic. We tackle these problems using a two-step approach in Section˜2.4. First, we estimate densities through diffusion-based kernel density estimation with fixed-point bandwidth selection. Second, we enforce quantile monotonicity using a lightweight quadratic programming procedure. These steps produce valid CDF and quantile function pairs that work well with Wasserstein-Fréchet regression and classification methods.

2.4 Quantile Function Estimation

In this section, we estimate quantile functions using diffusion-based kernel density estimation (Botev et al., 2010). For clarity, we temporarily omit the subject index $i$ in the observations $\mathbf{x}_{i}=(x_{i1},x_{i2},\ldots,x_{n_{i}})^{\top}$, and let $\mathbf{x}=(x_{1},x_{2},\ldots,x_{n})^{\top}$ when calculating the quantile estimation. The kernel density estimator (KDE) is a non-parametric method for estimating the probability density function of $\mathbf{x}$. The KDE of the true density $f$ is defined as

[TABLE]

where $K$ is the kernel function, $h>0$ is the bandwidth, and $n$ is the number of observations. The kernel $K$ is a symmetric probability density that integrates to 1; common choices of kernel $K$ include Gaussian, Epanechnikov, and uniform kernels. However, standard KDE has several problems. First, data-driven plug-in bandwidth selectors often assume normality, which is unrealistic in practice. Second, the Gaussian kernel lacks local adaptivity (Terrell and Scott, 1992). Third, most KDEs suffer boundary bias when the random variable has bounded support (Marron and Ruppert, 1994; Park et al., 2003).

The diffusion-based KDE addresses these issues. The basic strategy relies on approximating (7) via linear diffusion process smoothing (Botev et al., 2010):

[TABLE]

where $\mathcal{X}$ is the domain of $x$. The second-order linear differential operator $L$ is defined as $\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}x}(a(x)\frac{\mathrm{d}}{\mathrm{d}x}(\frac{\cdot}{p(x)}))$, where $a(x)$ and $p(x)$ are arbitrary positive functions on $\mathcal{X}$. Equation (8) describes an Itô diffusion in the pseudo-time $h$, and its solution yields the evolving density of $x$. In this paper we take $a(x)=p(x)\propto 1$, in which case the solution of (8) recovers the Gaussian-kernel KDE in (7). Next, we describe how to choose the bandwidth $h$ in the diffusion-based KDE. Let $t=h^{2}$. Using a first-order Taylor expansion, the asymptotic mean integrated squared error (AMISE) is $t^{2}{\left\lVert f^{\prime\prime}\right\rVert}^{2}/4+1/(2n\sqrt{\pi t})$ (Botev et al., 2010; Wand and Jones, 1995). The asymptotically optimal bandwidth is

[TABLE]

However, $f^{\prime\prime}$ is unknown in practice. Within diffusion-based KDE, the $L^{2}$ norm of the $j$-th derivative of $f$ can be expressed as

[TABLE]

where $\phi^{(2j)}(x_{k},x_{m};2t_{j})$ is the $2j$-th derivative of the Gaussian kernel $K(x_{k},x_{m};2t_{j})$ with respect to $x_{k}$ and $x_{m}$, and $t_{j}$ is the bandwidth for the $j$-th derivative (Botev et al., 2010, eq. (26)). Computing (10) requires $O(n^{2})$ operations, making it too slow for large $n$. To reduce complexity, Botev et al. (2010) exploit the discrete Fourier transform (DFT) to compute the $L^{2}$-norm of the $j$-th derivative. For the Gaussian kernel, the Fourier transform is also Gaussian, i.e., $\mathcal{F}[\phi](\omega,t)=e^{-\omega^{2}t/2}.$ The Fourier transform of $\hat{f}(x,t)$ is therefore $\mathcal{F}[\hat{f}](\omega,t)=\mathcal{F}[\phi](\omega,t)\cdot\hat{f}_{\text{emp}}(\omega)$ where $\hat{f}_{\text{emp}}(\omega)=\frac{1}{n}\sum_{i=1}^{n}e^{-i\omega x_{i}}$ is the empirical characteristic function. The power spectrum is ${\left\lvert\mathcal{F}[\hat{f}](\omega,t)\right\rvert}^{2}=e^{-\omega^{2}t}{\left\lvert\hat{f}_{\text{emp}}(\omega)\right\rvert}^{2}.$ Parseval’s theorem states that the $L^{2}$-norm of a function equals the $L^{2}$-norm of its Fourier transform (up to a $2\pi$ factor). Using $\mathcal{F}[f^{(j)}](\omega,t)=\omega^{j}\mathcal{F}[\hat{f}](\omega,t)$, we obtain

[TABLE]

In practice, the integral in (11) is approximated using a discrete cosine transform (DCT):

• •

Let $v_{k}$ be the Type-II DCT coefficients of the binned data;

• •

The squared DCT coefficients $a_{k}\coloneqq(v_{k}/2)^{2}$ approximate ${\left\lvert\hat{f}_{\text{emp}}(\omega)\right\rvert}^{2}$ at discrete frequencies $\omega_{k}=\frac{k\pi}{R},k=1,2,\ldots$, where $R$ is the data range;

• •

Let $I_{k}=k^{2}$ denote squared frequency terms (Laplacian eigenvalues in the DCT basis);

• •

Gaussian smoothing (i.e., $e^{-\omega^{2}t}$) in frequency space is approximated by $e^{-I_{k}\pi^{2}t}$.

Then, the $L^{2}$-norm of the $j$-th derivative of $f$ with bandwidth $t_{j}$ can be approximated as

[TABLE]

In addition, Proposition 2 in Botev et al. (2010) shows that the optimal bandwidth ${}_{*}t_{j}$ for the $j$-th derivative is

[TABLE]

Equations˜9, 12 and 13 motivate an $l$-stage algorithm to find the optimal bandwidth; see Algorithm 2.

We set the initial order $l=7$ following the recommendation in Botev et al. (2010).

Once the optimal bandwidth ${}_{*}t$ is obtained, the density $f$ is computed by smoothing $v_{k}$ in the frequency domain. Let $v_{k}^{*}\coloneqq v_{k}*\exp(-_{*}t(\pi k)^{2}/2)$. To construct the smoothed density functions, we apply the inverse DCT to the coefficients $v_{k}^{*},k=1,2,\ldots$. Furthermore, we need to normalise the smoothed density to integrate to 1. To estimate the cumulative distribution function (CDF), we reuse the DCT representation. With the assumption of a Gaussian kernel, the optimal bandwidth for smoothing the distribution function (Bowman et al., 1998) is

[TABLE]

By substituting $j=1$ and $t_{j}=\,_{*}t$ into (12), we directly obtain $\left\|f^{\prime}\right\|^{2}$. Hence, ${}_{*}t_{cdf}$ can be computed via (14). The CDF estimation procedure is given in Algorithm 3.

Finally, we compute the quantile function $Q(u)=F^{-1}(u)$ by using one-dimensional linear interpolation (e.g., the function “interp1d” in the Python package “scipy”) based on binned $\mathbf{x}$ and the estimated CDF in Algorithm˜3. To enforce quantile monotonicity, we directly refit the quantiles via the following optimization problem:

[TABLE]

where $\mathbf{q}^{*}=(q_{1}^{*},q_{2}^{*},\ldots,q_{m}^{*})$, $t_{j}$ is the $j$-th quantile weight, $q^{*}_{j}$ is the $j$-th refitted quantile of the estimated distribution function, and $q_{j}$ is the $j$-th interpolated quantile from Algorithm˜3. This problem can be solved using quadratic programming or Lagrange multipliers. Such monotone refitting is widely used to ensure valid quantile functions (Takeuchi et al., 2006; Bollaerts et al., 2006; Petersen et al., 2021). The code of Algorithms˜1, 2 and 3 is available in the Python package NeuroWAR.

3 Simulation Study

3.1 Kernel Density Estimation

In this section, we conduct a simulation study to evaluate the performance of the proposed method. First, we benchmark the implementation against a known probability density in the context of kernel density estimation (KDE). Pelz et al. (2023) proposed another diffusion-based KDE (diffKDE). The primary difference concerns the selection of the optimal bandwidth: diffKDE uses a numerical approximation, whereas our method employs the fixed-point algorithm of Botev et al. (2010).

Here, we use the mixture of Gaussian distributions (Pelz et al., 2023, eq. (39)) as the underlying density function:

[TABLE]

This density is a three-component Gaussian mixture with means $(6,9.5,12)$ and standard deviations $(1,0.7,0.5)$. The component weights are 0.3, 0.6, and 0.1, respectively. While this simulation employs a Gaussian mixture model for illustration, we emphasize that the subsequent Wasserstein-Fréchet regression framework makes no distributional assumptions. The method is entirely distribution-free and does not require knowledge of the underlying data. We use $g_{0}$ to denote the density estimate obtained by the diffusion-based KDE of Pelz et al. (2023), and $g_{1}$ to denote the density estimate from our procedure. We compare performance using the total variation distance. Let $p$ and $q$ be any density functions; the total variation distance is defined as $\delta(p,q)=\frac{1}{2}\int_{\mathbb{R}}\left|p(x)-q(x)\right|\mathrm{d}x$. Our algorithm is implemented in the Python package NeuroWAR, while the diffKDE implementation is available in the Python package diffKDE.

We consider sample sizes $n=50,100,200,400$. The densities $f$, $g_{0}$, $g_{1}$ and the corresponding distances $\delta(f,g_{0})$, $\delta(f,g_{1})$ are shown in Figure˜1. The results indicate that our method achieves uniformly smaller total variation distance than diffKDE across all sample sizes. As expected, the distance decreases with larger $n$, reflecting improved recovery of the underlying distribution.

3.2 Wasserstein-Fréchet Regression

This section presents simulation experiments to test our Wasserstein-Fréchet regression approach. We examine how well it performs compared to standard linear regression methods. First, we generate data according to the following model:

[TABLE]

where $\beta_{0}\sim N(0,0.5)$, $\beta_{1}\sim N(\nu_{1},\sigma_{1})$, and $\beta_{2}\sim N(\nu_{2},\sigma_{2})$ are mutually independent coefficients; ‘age’ denotes the subject’s age and ‘gender’ denotes the sex indicator (1 for male, 0 for female). $\nu_{1},\nu_{2},\sigma_{1},\sigma_{2}$ are the means and standard deviations of the normal distributions for $\beta_{1}$ and $\beta_{2}$. Equation˜15 specifies a linear relation between the response $\mu$ and covariates ‘age’ and ‘gender’. Given ‘age’ and ‘gender’, we randomly sample the tuple $(\beta_{0}^{(1)},\beta_{1}^{(1)},\beta_{2}^{(1)})$ as subject-specific coefficients. The conditional response is $\beta_{0}^{(1)}+\beta_{1}^{(1)}\cdot\text{age}+\beta_{2}^{(1)}\cdot\text{gender}$, with variance $0.25+\sigma_{1}^{2}\cdot\text{age}^{2}+\sigma_{2}^{2}\cdot~\text{gender}^{2}$. For brevity, write $\mu_{1|\cdot}=\beta_{0}^{(1)}+\beta_{1}^{(1)}\cdot\text{age}+\beta_{2}^{(1)}\cdot\text{gender}$ and $\sigma_{1|\cdot}^{2}=0.25+\sigma_{1}^{2}\cdot\text{age}^{2}+\sigma_{2}^{2}\cdot~\text{gender}^{2}$ for the conditional mean and variance given ‘age’ and ‘gender’. Moreover, $E[\mu]=\mu^{*}=\nu_{1}\cdot\text{age}+\nu_{2}\cdot\text{gender}$, so the density $N(\mu_{1|\cdot},\sigma_{1|\cdot})$ can be viewed as a random mean-shift of $N(\mu^{*},\sigma_{1|\cdot})$. We then sample observations from $N(\mu_{1|\cdot},\sigma_{1|\cdot})$ conditional on each subject’s ‘age’ and ‘gender’.

For the control group, we set $\nu_{1}=0.1,\nu_{2}=2,\sigma_{1}=0.1,\sigma_{2}=0.5$. For the mild traumatic brain injury (mTBI) group, we set $\nu_{2}=1,\sigma_{2}=0.5$, while $(\nu_{1},\sigma_{1})$ vary over the Cartesian product $\{(\nu_{1},\sigma_{1})\mid\nu_{1}\in A,\sigma_{1}\in B\}$, with $A=\{0.1,0.3,0.5,0.7\}$ and $B=\{0.1,0.2,0.3,0.4,0.5,0.6\}$. For each group, we simulate 2000 subjects with ‘age’ uniformly sampled from $[18,90]$ and ‘gender’ sampled from $\{0,1\}$. For each subject, we randomly draw one tuple $(\beta_{0}^{(1)},\beta_{1}^{(1)},\beta_{2}^{(1)})$ and then generate 1000 observations from $N(\mu_{1|\cdot},\sigma_{1|\cdot})$.

We apply the Wasserstein–Fréchet regression method to our simulated data following the approach outlined in Section 2. To establish a performance baseline, we additionally implement a standard linear regression approach. The linear baseline operates according to this specification:

We denote $y_{ij}$ as the observation for subject $i$ and sample $j$, where $i=1,2,\ldots,m$ and $j=1,2,\ldots,1000$. The subject-level summary is calculated as $\bar{y}_{i}=\sum_{j=1}^{1000}y_{ij}$, which represents the total response for each individual given their ‘age’ and ‘gender’. The linear regression model is

[TABLE]

where $\epsilon_{i}$ represents zero-mean Gaussian noise with variance $\sigma^{2}$. The coefficients $\alpha_{0},\alpha_{1},\alpha_{2}$ are estimated by ordinary least squares (OLS). Let $\hat{\alpha}_{0}^{1},\hat{\alpha}_{1}^{1},\hat{\alpha}_{2}^{1}$ and $\hat{\alpha}_{0}^{2},\hat{\alpha}_{1}^{2},\hat{\alpha}_{2}^{2}$ be the fitted coefficients for the control and mTBI groups, respectively. For an unseen subject with ‘age’ and ‘gender’, the corresponding group-wise predictions are $\hat{\bar{y}}_{1|\cdot}$ and $\hat{\bar{y}}_{2|\cdot}$. Given the subject’s average $\bar{y}_{0|\cdot}$, we classify this subject as control if ${\left\lvert\hat{\bar{y}}_{1|\cdot}-\bar{y}_{0|\cdot}\right\rvert}\leq{\left\lvert\hat{\bar{y}}_{2|\cdot}-\bar{y}_{0|\cdot}\right\rvert}$, and as mTBI otherwise. We evaluate linear regression and Wasserstein-Fréchet regression using classification accuracy, which is defined as $ACC=(TP+TN)/(TP+TN+FP+FN)$, where $TP$, $TN$, $FP$, and $FN$ denote true positives, true negatives, false positives, and false negatives. Data splitting follows a 70/30 training-testing ratio. The implementation code of Wasserstein-Fréchet regression and linear regression classifiers is available in the Python package NeuroWAR. The accuracies for Wasserstein-Fréchet regression and linear regression are shown in Figure˜2. As seen in Figure˜2, we can see that the Wasserstein–Fréchet regression method outperforms the traditional linear regression method in terms of accuracy. The accuracy of the Wasserstein–Fréchet regression method is higher than that of the traditional linear regression method, indicating that the Wasserstein–Fréchet regression method is more effective in capturing the underlying distribution of the data.

4 Real Data Analysis

In this section, we apply the proposed Wasserstein-Fréchet regression method to real-world neuroimaging data. The control dataset is obtained from the Cambridge Centre for Ageing and Neuroscience (Cam-CAN) data repository (Cam-CAN et al., 2014), containing approximately 700 subjects with T1-weighted structural MRI scans and 9-minute resting-state, eyes-closed magnetoencephalography (MEG) recordings. Subjects are aged between 18 and 87 years. The dataset is publicly available at Cam-CAN webpage. The mTBI dataset is obtained from Innovision IP Ltd. This dataset contains 143 subjects with T1-weighted structural MRI scans and two sessions of 5-minute resting-state, eyes-closed MEG recordings. Subjects are aged between 16 and 70 years. The dataset is not publicly available due to confidentiality agreements with patients. Both datasets were recorded using the MEGIN TRIUX system, a 306-channel MEG system with 102 magnetometers and 204 gradiometers. The goal of this analysis is to predict mTBI status based on demographic information (age, gender) and MRI/MEG data.

A promising MEG biomarker for mTBI detection is low-frequency power in resting-state recordings (Allen et al., 2021). The prominent finding is that low-frequency power in resting-state MEG is significantly higher in mTBI patients than in healthy controls (Lewine et al., 1999, 2007; Huang et al., 2014; Dunkley et al., 2015). Recent research (Dunkley et al., 2015; Huang et al., 2014; Wang et al., 2017; Marsh et al., 2025) has also demonstrated statistically significant spectral power increases in the gamma band. To investigate the joint effects of delta and gamma bands, we focus on the following frequency bands of interest: delta (0.5–4 Hz) (Hejazi et al., 2024), slow gamma (30–60 Hz), and fast gamma (60–100 Hz) (Bieri et al., 2014).

Next, we present a reproducible MRI+MEG analysis pipeline covering steps from raw data preprocessing to source magnitude imaging (Gramfort, 2013; Gross et al., 2013; Niso et al., 2019; Ferrante et al., 2022).

Anatomical reconstruction. We begin by constructing detailed brain models using FreeSurfer 8.0.0 (Fischl, 2012), an automated software suite for cortical surface reconstruction and volumetric segmentation. The T1-weighted structural MRI scans undergo the standard FreeSurfer pipeline ‘recon-all’. Following anatomical reconstruction, we construct Boundary Element Method (BEM) head models to obtain realistic forward models for MEG source analysis. BEM head model construction is performed using MNE-Python (Gramfort, 2013).

Data preprocessing. The raw MEG data undergo comprehensive preprocessing using MNE-Python (Gramfort, 2013). The preprocessing pipeline comprises several sequential steps: (1) Signal-Space Separation (SSS) is applied to eliminate environmental noise and artifacts from MEG recordings (Taulu and Simola, 2006); (2) Temporal cropping is performed to extract relevant time intervals, specifically from stimulus onset to termination when stimulus channels are present; (3) Power line artifact detection and removal is conducted to eliminate environmental interference manifesting as persistent oscillations at AC power line frequencies; (4) Temporal segmentation of continuous recordings into 2-second epochs is performed, excluding previously identified bad channels from analysis; (5) Computation of the noise covariance matrix is performed using approximately 5 minutes of empty room recordings to characterise environmental noise baseline for subsequent source modelling. After preprocessing, the control group (Cam-CAN) contains 588 healthy subjects, and the mTBI group contains 131 subjects.

Source magnitude imaging. Source magnitude imaging is performed using vector-based spatial-temporal analysis (VESTAL, Huang et al., 2014). As the VESTAL algorithm is based on the $L_{1}$-minimum-norm solution, it is robust to environmental noise and artifacts in the data. MEG source magnitude imaging is carried out on preprocessed MEG data using the BEM head model. For each 2-second epoch, we apply band-pass filtering for the following frequency bands: delta (0.5–4 Hz), slow gamma (30–60 Hz), and fast gamma (60–100 Hz). For each frequency band, we obtain source magnitude images using the VESTAL algorithm (Huang et al., 2014). The brain volume is divided into a grid of 5120 nodes in the BEM model, meaning that the source magnitude image is a 5120-dimensional vector for each 2-second epoch. Rather than analyzing individual nodes, we focus on regions of interest (ROIs) defined by the Desikan-Killiany atlas (Desikan et al., 2006), which yields 68 regions across left and right brain hemispheres. For each region, we compute the average source magnitude across all nodes within the region. Finally, we concatenate the average source magnitude vectors across all epochs to form the final feature representation for each subject.

Wasserstein-Fréchet regression. First, training and test datasets are generated as follows: For the Cam-CAN dataset, we randomly select 131 subjects as the test dataset, with the remainder forming the training dataset. For the Innovision dataset, we randomly select 65 subjects as the test dataset, with the remaining 66 subjects forming the training dataset.

For each subject in the training datasets and each frequency band, we estimate quantile functions of magnitude using our diffusion-based kernel density estimation across all epochs for each brain area. We discretise the quantile function on the support $[0,1]$ with equal-grid spacing (1/1024), resulting in a finite number of quantile levels. Wasserstein-Fréchet regression from Section 2 is then performed on these quantile functions to model the relationship between MEG features and demographic information (age, gender).

In the second step, we generate predictions for each brain region and subject in our training data. We examine three frequency bands: delta, slow gamma, and fast gamma. The true labels are stored in vector $Y_{0}$, where $Y_{0,i}=1$ indicates mTBI and $Y_{0,i}=0$ indicates control status for subject $i$. Our prediction results create a binary matrix $Z_{0}$ of size $(457+66,68\times 3)$.We consider using random forest to build two models: $Y_{0}\sim Z_{0}[:,1:68]$ and $Y_{0}\sim Z_{0}$, where $Z_{0}[:,1:68]$ represents predicted status in the delta band.

Finally, we use the Wasserstein-Fréchet regression model to obtain a binary matrix $Z_{1}$ with dimensions $(196,204)$ for the test dataset. Let $Y_{1}$ be the binary vector of underlying labels of test subjects. We use the trained random forest models to predict the status of each subject in the test dataset. Define $\hat{Y}_{1,\text{delta}}$ and $\hat{Y}_{1,\text{all}}$ as the predicted labels for delta-only and all-band analyses, respectively. We measure performance using standard classification accuracy $ACC$. Results show that delta-band classification achieves 0.8893 accuracy. When incorporating all frequency bands, performance improves substantially to 0.9802.

To provide a more comprehensive comparison with the linear and SVM classifiers mentioned in the literature, we trained and evaluated both classifiers on the same dataset used for our Wasserstein-Fréchet framework. First, since diffusion-based KDE can handle varying numbers of epochs across frequency bands and brain regions, while linear and SVM require balanced input data, we modified our feature extraction approach. Instead of using the full KDE representation, we computed the average spectral magnitude across all epochs for each region and frequency band as input features for the SVM and linear classifier.

Second, for the linear classifier baseline, we fitted the linear regression model from (16) (described in Section˜3.2) to the balanced preprocessed data. This approach generates binary predictions for individual brain regions and frequency bands. Following the same ensemble strategy as our Wasserstein-Fréchet approach, we employed random forest to aggregate individual predictions into final subject-level classifications. In essence, this baseline substitutes ordinary linear regression for the Wasserstein-Fréchet regression component while maintaining the same overall architecture.

Third, for the SVM classifier baseline, we constructed feature vectors using age, gender, and the averaged spectral magnitudes from each frequency band and brain region. For the delta-only analysis, this resulted in 70 features (68 brain regions plus age and gender). When incorporating all three frequency bands, the feature space expanded to 206 dimensions (68 regions $\times$ 3 bands + 2 demographic variables). In our study, before implementation of SVM classification, we applied a standard scaler, i.e., standardised the features to have zero mean and unit variance, to the training and test data. We used the radial basis function (RBF) kernel for SVM classification, as it is effective in handling non-linear relationships in high-dimensional spaces (Vapnik, 1995; Chang and Lin, 2011). It is easy to implement the SVM classification in Python using libraries such as scikit-learn. Finally, when analyzing only the delta band, the linear and SVM classifiers achieved test accuracies of 0.8776 and 0.9031, respectively. For the analysis incorporating all three frequency bands, the linear and SVM classifiers obtained test accuracies of 0.9388 and 0.9337, respectively. Both accuracies are lower than 0.9802 achieved by our Wasserstein-Fréchet regression approach. This is not surprising because the linear and SVM classifiers use only the first moment of distributions instead of the full density. The linear and SVM classifiers discards higher-order moments and distributional shape information that may be important for classification.

From the Table˜1, we can see that our method achieves the highest prediction accuracy of 0.9802, which is significantly higher than the previous methods. This indicates that the Wasserstein-Fréchet regression method is effective in capturing the underlying distribution of the source magnitude across epochs and can be used as a powerful tool for mTBI detection.

5 Conclusion

This study introduced a quantile function-based representation of epoch-wise MEG source magnitude distributions and a global Wasserstein-Fréchet regression framework to link these functional responses to demographic covariates. Methodologically, the work contributes (i) a diffusion-based kernel density procedure with fixed-point bandwidth selection and fast spectral approximation to estimate subject- and region-specific quantile functions; (ii) a covariate-adjusted Wasserstein regression that yields conditional mean distributions without local smoothing or tuning; and (iii) a two-stage classification strategy that aggregates region-level predictions across three frequency bands. The pipeline relies on anatomically realistic forward modelling using BEM, and utilises $L_{1}$-based source magnitude imaging (VESTAL) to stabilise estimates under realistic noise conditions.

Simulation results demonstrate that the proposed density estimation reduces total variation error relative to a recent diffusion-based KDE alternative across sample sizes. Wasserstein-Fréchet regression outperforms standard linear regression when the response is a fundamental quantile function. On real data combining Cam-CAN controls and mTBI cases, the approach attains high predictive performance, with accuracy of 0.9802 when integrating delta, slow gamma, and fast gamma bands, exceeding delta-only performance. Our result also align with reports of low-frequency and gamma-band abnormalities in mTBI. These gains reflect the joint benefits of (a) preserving full distributional information via quantile functions, (b) modelling covariate effects at the distribution level, (c) multiband aggregation that mitigates frequency-specific leakage and enhances robustness.

Limitations include the need for broader external validation beyond a single cohort, assessment of generalisability across sites, scanners, and preprocessing choices (e.g., SSS parameters, epoching, ROI definitions). Future work will pursue multicenter validation, joint modelling of power and connectivity. Overall, the results support Wasserstein-Fréchet regression on quantile representations as a principled and effective framework for MEG biomarkers in mTBI detection, with clear translational potential given its accuracy, modularity, and compatibility with standard neuroimaging.

Data and Code Availability

The control data comes from the Cambridge Centre for Ageing and Neuroscience (Cam-CAN) dataset (Cam-CAN et al., 2014). The code for the simulation studies are available in the Python package NeuroWAR.

Competing interests

The authors from the University of Kent, J.L., J.Z.. have no competing interests. G.G. is an employee of Innovision IP Ltd, which provides commercial reports on individuals who may have had a head injury.

Funding

This research received financial support jointly from UKRI through Innovate UK and Innovision IP Ltd. (Grant No: 10053865).

Acknowledgments

This work was supported by the High-Performance Computing Cluster at the University of Kent.