Monotonic Gaussian Process for Spatio-Temporal Disease Progression Modeling in Brain Imaging Data

TL;DR

This paper presents a probabilistic model combining spatio-temporal matrix factorization and monotonic Gaussian Processes to analyze high-dimensional brain imaging data, effectively disentangling disease progression patterns.

Contribution

It introduces a novel spatio-temporal Gaussian Process model with monotonic constraints for brain disease progression analysis, incorporating anatomically plausible priors and multi-scale sparse spatial coding.

Findings

Successfully disentangles disease trajectories in synthetic data

Reveals key brain regions involved in neurodegeneration

Identifies disease-specific temporal progression patterns

Abstract

We introduce a probabilistic generative model for disentangling spatio-temporal disease trajectories from series of high-dimensional brain images. The model is based on spatio-temporal matrix factorization, where inference on the sources is constrained by anatomically plausible statistical priors. To model realistic trajectories, the temporal sources are defined as monotonic and time-reparametrized Gaussian Processes. To account for the non-stationarity of brain images, we model the spatial sources as sparse codes convolved at multiple scales. The method was tested on synthetic data favourably comparing with standard blind source separation approaches. The application on large-scale imaging data from a clinical study allows to disentangle differential temporal progression patterns mapping brain regions key to neurodegeneration, while revealing a disease-specific time scale associated to the clinical diagnosis.

Tables

Table 1. Table 1 : MSE and SSIM between respectively the ground truth temporal and spatial sources with respect to the ones estimated by the different methods.

	Temporal (MSE)	Spatial (SSIM)
MGPA	$\left(\mathrm{𝟖}\pm \mathrm{𝟒}\right){.10}^{-\mathrm{𝟓}}$	$\mathrm{𝟗𝟖}\%\pm \mathrm{𝟏}$
ICA	$\left(6\pm 3\right){.10}^{-4}$	$97\%\pm 2$
NMF	$\left(3\pm 2\right){.10}^{-2}$	$40\%\pm 17$
PCA	$0.44\pm {10}^{-3}$	$15\%\pm 1$

Table 2. Table 2 : MSE and SSIM between respectively the ground truth temporal and spatial sources with respect to the ones estimated by MGPA. R 2 coefficient of the linear regression between the original time-line and the estimated time-shift.

	Temporal (MSE)	Spatial (SSIM)	R2
MGPA	$\left(\mathrm{𝟐}\pm 0.8\right){.10}^{-\mathrm{𝟐}}$	$\mathrm{𝟗𝟓}\%\pm \mathrm{𝟒}$	$0.98\pm 0.005$

Table 3. Table 3 : Baseline socio-demographic and clinical information for study cohort. Average values and standard deviation in parenthesis. NL: normal individuals, NL converter: normal subjects who converted to MCI or to AD, MCI: mild cognitive impairment, MCI converter: MCI subjects who converted to AD, AD: Alzheimer’s patients. ADAS13: Alzheimer’s Disease Assessment Scale-cognitive subscale, 13 items. FAQ: Functional Assessment Questionnaire. FDG: (18)F-fluorodeoxyglucose Positron Emission Tomography (PET) imaging. AV45: (18)F-florbetapir Amyloid PET imaging.

Group	NL	\pbox20cm      NL
converter	MCI	\pbox20cm     MCI
converter	AD
N	103	34	164	129	114
Age	73 (6)	78 (5)	73 (7)	73 (7)	74 (8)
Education (yrs)	16.3 (3)	16 (3)	15.7 (3)	16 (3)	15.6 (3)
ADAS13	9.1 (4.4)	11.4 (4.3)	14.6 (5.5)	20.4 (6.5)	31.6 (8.5)
FAQ	0.3 (0.7)	0.2 (0.6)	1.9 (2.8)	5.0 (4.6)	13.5 (6.9)
Entorhinal (cm3)	3.8 (0.5)	3.5 (0.5)	3.6 (0.6)	3.2 (0.7)	2.8 (0.6)
Hippocampus (cm3)	7.4 (0.9)	6.9 (0.7)	6.9 (0.9)	6.4 (0.9)	5.9 (0.8)
Ventricles (cm3)	31 (16)	42 (21)	39 (23)	40 (19)	48 (23)
Whole brain (cm3)	1033 (104)	1019 (91)	1058 (103)	1037 (102)	1005 (115)
FDG	1.3 (0.1)	1.3 (0.1)	1.2 (0.1)	1.1 (0.1)	1.0 (0.1)
AV45	1.3 (0.2)	1.3 (0.1)	1.3 (0.2)	1.4 (0.2)	1.5 (0.2)

Table 4. Table 1 : Learning rates (LR) of the different parameters of the model.

	$𝜽$	$𝑴$	$𝑷$	$𝒁$	$\sigma ,\nu$	$𝜹$
lr	${10}^{-2}$	${10}^{-3}$	${10}^{-1}$	${10}^{-1}$	${10}^{-2}$	${10}^{-4}$

Table 5. Table 1 : MSE and SSIM between respectively the ground truth temporal and spatial sources with respect to the ones estimated by the different standard methods.

	Temporal (MSE)	Spatial (SSIM)
ICA	$0.24\pm 0.08$	$54\%\pm 2$
NMF	$0.25\pm 0.03$	$22\%\pm 14$
PCA	$0.66\pm 0.05$	$9\%\pm 3$

Equations

$${\mathbf{f}}_p(\tau )=\tau +{\delta }_p.$$

$$\begin{array}{rl}& {\mathbf{Y}}_{p:}^m({\mathbf{f}}_p(\tau ),{\theta }_m,{\psi }_m)={\mathbf{S}}_{p:}^m({\mathbf{f}}_p(\tau ),{\theta }_m){\mathbf{A}}^m({\psi }_m)+{\mathbf{Z}}_{p:}^m+{\mathcal{E}}_{p:}^m,\\ & {\mathbf{V}}_{p,c}({\mathbf{f}}_p(\tau ),{\theta }_c)={\mathbf{U}}_{p,c}({\mathbf{f}}_p(\tau ),{\theta }_c)+{\mathcal{H}}_{p,c}.\end{array}$$

$$\displaystyle\begin{split}\displaystyle p(\boldsymbol{Y},\boldsymbol{V}|\boldsymbol{A},\boldsymbol{S},\boldsymbol{Z},\boldsymbol{U},\boldsymbol{\delta},\sigma,\nu)&=\Big{(}\prod_{m}p(\boldsymbol{Y}^{m}|\boldsymbol{A}^{m},\boldsymbol{S}^{m},\boldsymbol{Z}^{m},\boldsymbol{\delta},\sigma_{m})\Big{)}\Big{(}\prod_{c}p(\boldsymbol{V}_{:c}|\boldsymbol{U}_{:c},\boldsymbol{\delta},\nu_{c})\Big{)}.\end{split}$$

$$\displaystyle\begin{split}p(\boldsymbol{Y}^{m}|\boldsymbol{A}^{m},\boldsymbol{S}^{m},\boldsymbol{Z}^{m},\boldsymbol{\delta},\sigma_{m})=\Big{(}\prod_{p}&\frac{1}{(2\pi\sigma_{m}^{2})^{\frac{F_{m}}{2}}}\exp(-\frac{1}{2\sigma_{m}^{2}}||\boldsymbol{Y}_{p:}^{m}(\boldsymbol{f}_{p}(\tau),\theta_{m},\psi_{m})\\ &-\boldsymbol{S}_{p:}^{m}(\boldsymbol{f}_{p}(\tau),\theta_{m})\boldsymbol{A}^{m}(\psi_{m})-\boldsymbol{Z}_{p:}^{m}||^{2})\Big{)}.\end{split}$$

$$p(\mathcal{C}|{\mathbf{S}}^{\prime },\gamma )=\prod _{p,n}(1+\exp (-\gamma {\mathbf{S}}_{p,n}^{\prime }(\mathbf{t})){)}^{-1},$$

$${\mathbf{S}}_{:n}^{\prime }(t)=\frac{d\varphi (\mathbf{t}({\mathbf{\omega }}^n{)}^T)}{d\mathbf{t}}{\mathbf{w}}^n.$$

$$q_1(\mathbf{B})=\prod _{n=1}^N\mathcal{N}({\mathbf{M}}_{n:},diag({\alpha }_{n,1}{\mathbf{M}}_{n,1}^2...{\alpha }_{n,F}{\mathbf{M}}_{n,F}^2)).$$

$$\begin{array}{rl}\log (p(\mathbf{Y},\mathbf{V},\mathcal{C}|\mathbf{Z},\mathbf{\delta },\sigma ,\nu ,\gamma ))& =\sum _m\log (p({\mathbf{Y}}^m,{\mathcal{C}}^m|{\mathbf{Z}}^m,\mathbf{\delta },{\sigma }_m,{\gamma }_m))+\sum _c\log (p({\mathbf{V}}_{:c},{\mathcal{C}}^c|\mathbf{\delta },{\nu }_c,{\gamma }_c)).\end{array}$$

$$\begin{array}{rl}\log (p({\mathbf{Y}}^m,{\mathcal{C}}^m|{\mathbf{Z}}^m,\mathbf{\delta },{\sigma }_m,{\gamma }_m))⩾& E_{q_1,q_2,q_3}[\log (p({\mathbf{Y}}^m|{\mathbf{B}}^m,{\mathbf{\Omega }}^m,{\mathbf{W}}^m,{\mathbf{Z}}^m,\mathbf{\delta },{\sigma }_m))]\\ & +E_{q_2,q_3}[\log (p({\mathcal{C}}^m|{\mathbf{\Omega }}^m,{\mathbf{W}}^m,\mathbf{\delta },{\gamma }_m))]\\ & -\mathcal{D}[q_1({\mathbf{B}}^m)||p({\mathbf{B}}^m)]-\mathcal{D}[q_2({\mathbf{\Omega }}^m)||p({\mathbf{\Omega }}^m)]-\mathcal{D}[q_3({\mathbf{W}}^m)||p({\mathbf{W}}^m)],\\ \log (p({\mathbf{V}}_{c:},{\mathcal{C}}^c|\mathbf{\delta },{\nu }_c,{\gamma }_c))⩾& E_{q_2,q_3}[\log (p({\mathbf{V}}_{c:}|{\mathbf{\Omega }}^c,{\mathbf{W}}^c,\mathbf{\delta },{\sigma }_c))]\\ & +E_{q_2,q_3}[\log (p({\mathcal{C}}^c|{\mathbf{\Omega }}^c,{\mathbf{W}}^c,\mathbf{\delta },{\gamma }_c))]\\ & -\mathcal{D}[q_2({\mathbf{\Omega }}^c)||p({\mathbf{\Omega }}^c)]-\mathcal{D}[q_3({\mathbf{W}}^c)||p({\mathbf{W}}^c)]\end{array}$$

$$\begin{array}{rl}\log (p(\mathbf{Y},\mathbf{V},\mathcal{C}|\mathbf{Z},\mathbf{\delta },\sigma ,\nu ,\gamma ))& ⩾\sum _m{\mathcal{L}}_m+\sum _c{\mathcal{L}}_c.\end{array}$$

$$\begin{array}{rl}{q}_2({\mathbf{\Omega }}^m)& =\prod _{n=1}^{N_m}{q}_2({\mathbf{\omega }}^n{)}^m=\prod _{n=1}^{N_m}\prod _{j=1}^{N_{rf}}\mathcal{N}({\mathbf{R}}_{n,j},{\mathbf{Q}}_{n,j}^2{)}^m,\\ q_3({\mathbf{W}}^m)& =\prod _{n=1}^{N_m}{q}_3({\mathbf{w}}^n{)}^m=\prod _{n=1}^{N_m}\prod _{j=1}^{N_{rf}}\mathcal{N}({\mathbf{T}}_{n,j},{\mathbf{V}}_{n,j}^2{)}^m,\end{array}$$

$$\begin{array}{rl}\mathbf{\theta }& =\{{\mathbf{R}}_{n:}^m,{\mathbf{Q}}_{n:}^m,{\mathbf{T}}_{n:}^m,{\mathbf{V}}_{n:}^m,l_n,n\in [1,N_m]{\}}_{m=1}^M\cup \{{\mathbf{R}}_{c:},{\mathbf{Q}}_{c:},{\mathbf{T}}_{c:},{\mathbf{V}}_{c:}^c,l_c{\}}_{c=1}^C,\\ \mathbf{\psi }& =\{{\mathbf{M}}_{n:}^m,{\mathbf{P}}_{n:}^m,n\in [1,N_m]{\}}_{m=1}^M.\end{array}$$

$$\displaystyle\begin{split}\log(p(\boldsymbol{Y}^{m},\mathcal{C}^{m}|\boldsymbol{Z}^{m},\boldsymbol{\delta},\sigma_{m},\gamma_{m}))=\log\Big{[}\int&p(\boldsymbol{Y}^{m}|\boldsymbol{B}^{m},\boldsymbol{S}^{m},\boldsymbol{Z}^{m},\boldsymbol{\delta},\sigma_{m})p(\mathcal{C}|\frac{d\boldsymbol{S}^{m}}{d\boldsymbol{t}},\boldsymbol{\delta},\gamma_{m})p(\boldsymbol{B}^{m})\\ &p(\boldsymbol{S}^{m},\frac{d\boldsymbol{S}^{m}}{d\boldsymbol{t}}|\boldsymbol{\delta},\gamma)d\boldsymbol{B}^{m}d\boldsymbol{S}^{m}\Big{]}\\ &=\log\Big{[}\int p(\boldsymbol{Y}^{m}|\boldsymbol{B}^{m},\boldsymbol{S}^{m},\boldsymbol{Z}^{m},\boldsymbol{\delta},\sigma_{m})p(\mathcal{C}|\frac{d\boldsymbol{S}^{m}}{d\boldsymbol{t}},\boldsymbol{\delta},\gamma_{m})p(\boldsymbol{B}^{m})\\ &p(\frac{d\boldsymbol{S}^{m}}{d\boldsymbol{t}}|\boldsymbol{S}^{m},\boldsymbol{\delta},\gamma)p(\boldsymbol{S}^{m})d\boldsymbol{B}^{m}d\boldsymbol{S}^{m}\Big{]}.\end{split}$$

$$\displaystyle\begin{split}\log(p(\boldsymbol{Y}^{m},\mathcal{C}^{m}|\boldsymbol{Z}^{m},\boldsymbol{\delta},\sigma_{m},\gamma_{m}))=\log\Big{[}\int&p(\boldsymbol{Y}^{m}|\boldsymbol{B}^{m},\boldsymbol{S}^{m},\boldsymbol{Z}^{m},\boldsymbol{\delta},\sigma_{m})p(\mathcal{C}|\frac{d\boldsymbol{S}^{m}}{d\boldsymbol{t}},\boldsymbol{\delta},\gamma_{m})p(\boldsymbol{B}^{m})\\ &p(\boldsymbol{S}^{m})d\boldsymbol{B}^{m}d\boldsymbol{S}^{m}\Big{]}.\end{split}$$

$$\log (p({\mathbf{Y}}^m,{\mathcal{C}}^m|{\mathbf{Z}}^m,\mathbf{\delta },{\sigma }_m,{\gamma }_m))$$

$$\displaystyle\hskip 45.00006ptp(\boldsymbol{S}^{m})d\boldsymbol{B}^{m}d\boldsymbol{S}^{m}\Big{]}$$

$$\displaystyle=\log\Big{[}\int p(\boldsymbol{Y}^{m}|\boldsymbol{B}^{m},\boldsymbol{\Omega}^{m},\boldsymbol{W}^{m},\boldsymbol{Z}^{m},\boldsymbol{\delta},\sigma_{m})p(\mathcal{C}|\boldsymbol{\Omega}^{m},\boldsymbol{W}^{m},\boldsymbol{\delta},\gamma_{m})p(\boldsymbol{B}^{m})$$

$$\displaystyle\hskip 45.00006ptp(\boldsymbol{\Omega}^{m})p(\boldsymbol{W}^{m})d\boldsymbol{B}^{m}d\boldsymbol{\Omega}^{m}d\boldsymbol{W}^{m}\Big{]}$$

$$\displaystyle=\log\Big{[}\int p(\boldsymbol{Y}^{m}|\boldsymbol{B}^{m},\boldsymbol{\Omega}^{m},\boldsymbol{W}^{m},\boldsymbol{Z}^{m},\boldsymbol{\delta},\sigma_{m})p(\mathcal{C}|\boldsymbol{\Omega}^{m},\boldsymbol{W}^{m},\boldsymbol{\delta},\gamma_{m})p(\boldsymbol{B}^{m})$$

$$\displaystyle\hskip 45.00006ptp(\boldsymbol{\Omega}^{m})p(\boldsymbol{W}^{m})\frac{q_{1}(\boldsymbol{B}^{m})q_{2}(\boldsymbol{\Omega}^{m})q_{3}(\boldsymbol{W}^{m})}{q_{1}(\boldsymbol{B}^{m})q_{2}(\boldsymbol{\Omega}^{m})q_{3}(\boldsymbol{W}^{m})}d\boldsymbol{B}^{m}d\boldsymbol{\Omega}^{m}d\boldsymbol{W}^{m}\Big{]}$$

$$\displaystyle=\log\Big{[}\operatorname{\mathbb{E}}_{q_{1},q_{2},q_{3}}\frac{p(\boldsymbol{Y}^{m}|\boldsymbol{B}^{m},\boldsymbol{\Omega}^{m},\boldsymbol{W}^{m},\boldsymbol{Z}^{m},\boldsymbol{\delta},\sigma_{m})p(\mathcal{C}|\boldsymbol{\Omega}^{m},\boldsymbol{W}^{m},\boldsymbol{\delta},\gamma_{m})}{q_{1}(\boldsymbol{B}^{m})q_{2}(\boldsymbol{\Omega}^{m})q_{3}(\boldsymbol{W}^{m})}$$

$$\displaystyle\hskip 70.0001pt\frac{p(\boldsymbol{B}^{m})p(\boldsymbol{\Omega}^{m})p(\boldsymbol{W}^{m})}{q_{1}(\boldsymbol{B}^{m})q_{2}(\boldsymbol{\Omega}^{m})q_{3}(\boldsymbol{W}^{m})}\Big{]}$$

$$\displaystyle\geq\operatorname{\mathbb{E}}_{q_{1},q_{2},q_{3}}\Big{(}\log\Big{[}\frac{p(\boldsymbol{Y}^{m}|\boldsymbol{B}^{m},\boldsymbol{\Omega}^{m},\boldsymbol{W}^{m},\boldsymbol{Z}^{m},\boldsymbol{\delta},\sigma_{m})p(\mathcal{C}|\boldsymbol{\Omega}^{m},\boldsymbol{W}^{m},\boldsymbol{\delta},\gamma_{m})}{q_{1}(\boldsymbol{B}^{m})q_{2}(\boldsymbol{\Omega}^{m})q_{3}(\boldsymbol{W}^{m})}$$

$$\displaystyle\hskip 76.00017pt\frac{p(\boldsymbol{B}^{m})p(\boldsymbol{\Omega}^{m})p(\boldsymbol{W}^{m})}{q_{1}(\boldsymbol{B}^{m})q_{2}(\boldsymbol{\Omega}^{m})q_{3}(\boldsymbol{W}^{m})}\Big{]}\Big{)}$$

$$=E_{q_1,q_2,q_3}[\log (p({\mathbf{Y}}^m|{\mathbf{B}}^m,{\mathbf{\Omega }}^m,{\mathbf{W}}^m,{\mathbf{Z}}^m,\mathbf{\delta },{\sigma }_m))]$$

$$+E_{q_2,q_3}[\log (p({\mathcal{C}}^m|{\mathbf{\Omega }}^m,{\mathbf{W}}^m,\mathbf{\delta },{\gamma }_m))]$$

$$-\mathcal{D}[q_1({\mathbf{B}}^m)||p({\mathbf{B}}^m)]-\mathcal{D}[q_2({\mathbf{\Omega }}^m)||p({\mathbf{\Omega }}^m)]-\mathcal{D}[q_3({\mathbf{W}}^m)||p({\mathbf{W}}^m)].$$

$$\begin{array}{rl}& \mathcal{D}[q_1(\mathbf{B})||p(\mathbf{B})]=\sum _m\mathcal{D}[q_1({\mathbf{B}}^m)||p({\mathbf{B}}^m)],\\ & \mathcal{D}[q_2(\mathbf{\Omega })||p(\mathbf{\Omega })]=\sum _m\mathcal{D}[q_1({\mathbf{\Omega }}^m)||p({\mathbf{\Omega }}^m)]+\sum _c\mathcal{D}[q_1({\mathbf{\Omega }}^c)||p({\mathbf{\Omega }}^c)],\\ & \mathcal{D}[q_3(\mathbf{W})||p(\mathbf{W})]=\sum _m\mathcal{D}[q_3({\mathbf{W}}^m)||p({\mathbf{W}}^m)]+\sum _c\mathcal{D}[q_3({\mathbf{W}}^c)||p({\mathbf{W}}^c)].\end{array}$$

$$-\mathcal{D}[q_1(\mathbf{B})||p(\mathbf{B})]=\sum _{n,f}$$

$$\begin{array}{rl}& \mathcal{D}[q_2(\mathbf{\Omega })|p(\mathbf{\Omega })]=\frac{1}{2}\sum _{n,j}{\mathbf{Q}}_{n,j}^2{l}_n+{\mathbf{R}}_{n,j}^2{l}_n-1-\log ({\mathbf{Q}}_{n,j}^2{l}_n),\\ & \mathcal{D}[q_3(\mathbf{W})|p(\mathbf{W})]=\frac{1}{2}\sum _{n,j}{\mathbf{V}}_{n,j}^2+{\mathbf{T}}_{n,j}^2-1-\log ({\mathbf{V}}_{n,j}^2).\end{array}$$

$$\begin{array}{rl}\mathbf{\theta }& =\{{\mathbf{R}}_{n:}^m,{\mathbf{Q}}_{n:}^m,{\mathbf{T}}_{n:}^m,{\mathbf{V}}_{n:}^m,l_n,n\in [1,N_m]{\}}_{m=1}^M\cup \{{\mathbf{R}}_{c:},{\mathbf{Q}}_{c:},{\mathbf{T}}_{c:},{\mathbf{V}}_{c:},l_c,{\}}_{c=1}^C,\\ \mathbf{\psi }& =\{{\mathbf{M}}_{n:}^m,{\mathbf{P}}_{n:}^m,n\in [1,N_m]{\}}_{m=1}^M.\end{array}$$

Full text

††journal: NeuroImage

1 Introduction

Neurodegenerative disorders such as Alzheimer’s disease (AD) are characterized by morphological and molecular changes of the brain, ultimately leading to cognitive and behavioral decline. Clinicians suggested hypothetical models of the disease evolution, showing how different types of biomarkers interact and lead to the final dementia stage [15]. In the past years, efforts have been made in order to collect large databases of imaging and clinical measures, hoping to obtain more insights about the disease progression through data-driven models describing the trajectory of the disease over time. This kind of models are of critical importance for understanding the pathological progression in large scale data, and would represent a valuable reference for improving the individual diagnosis.

Current clinical trials in AD are based on longitudinal monitoring of biomarkers. Disease progression modelling aims at providing an interpretable way of modelling the evolution of biomarkers according to an estimated history of the pathology, as proposed for example in [10], [12], [16], [23], and [41]. Therefore, disease progression models are promising methods for automatically staging patients, and quantifying their progression with respect to the underlying model of the pathology. These approaches entail a great potential for automatic stratification of individuals based on their estimated stage and progression speed, and for assessment of efficacy of disease modifying drugs. Within this context, we propose a spatio-temporal generative model of disease progression, aimed at disentangling and quantifying the independent dynamics of changes observed in datasets of multi-modal data. With this term we indicate data acquired via different imaging modalities such as Magnetic Resonance Imaging (MRI) or Positron-Emission Tomography (PET), as well as non-imaging data such as clinical scores assessed by physicians. Moreover, we aim at automatically inferring the disease severity of a patient with respect to the estimated trajectory. Defining such a disease progression model raises a number of methodological challenges.

AD spreads over decades with a temporal mismatch between the onset of the disease and the moment where the clinical symptoms appear. Either age of diagnosis, or the chronological age, are therefore not suitable as a temporal reference to describe the disease progression in time. Moreover, as the follow-up of patients doesn’t exceed a few years, the development of a model of long-term pathological changes requires to integrate cross-sectional data from different individuals, in order to consider a longer period of time. In virtue of the lack of a well defined temporal reference, observations from different individuals are characterized by large and unknown variability in the onset and speed of the disease. It is therefore necessary to account for a time-reparameterization function, mapping each individuals’ observations to a common temporal axis associated to the absolute disease trajectory [16, 36]. This would allow to estimate an absolute time-reference related to the natural history of the pathology.

The analysis of MRI and PET data, requires to account for spatio-temporally correlated features (voxels, i.e. volumetric pixels) defined over arrays of more than a million entries. The development of inference schemes jointly considering these correlation properties thus raises scalability issues, especially when accounting for the non-stationarity of the image signal. Furthermore, the brain regions involved in AD exhibit various dynamics in time, and evolve at different speed [39]. From a modeling perspective, accounting for differential trajectories over space and time raises the problem of source identification and separation. This issue has been widely addressed in neuroimaging via Independent Component Analysis (ICA) [8], especially on functional MRI (fMRI) data [7]. Nevertheless, while fMRI time-series are usually defined over a few hundreds of time points acquired per subject, our problem consists in jointly analyzing short-term and cross-sectional data observations with respect to an unknown time-line. This problem cannot be tackled with standard ICA, as time is generally an independent variable on which inference is not required. Moreover, ICA retrieves spatial sources based on the assumption of statistical independence. This assumption does not necessarily lead to clinically interpretable findings. Indeed, dependency across temporal patterns can be still highly relevant to the pathology, for example when modeling temporal delay across similar sources.

The problem of providing a realistic description of the biological processes is critical when analyzing biomedical data, such as medical images. For example, to describe a plausible evolution of AD from normal to pathological stages, smoothness and monotonicity are commonly assumed for the temporal sources. It is also necessary to account for the non-stationarity of changes affecting the brain from global to localized spatio-temporal processes. As a result, spatial sources need to account for different resolutions at which these changes take place. While several multi-scale analysis approaches have been proposed to model spatio-temporal signals [26, 6, 14], extending this type of methods to the high-dimension of medical images is generally not trivial due to scalability issues. Finally, the noisy nature of medical images, along with the large signal variability across observations, requires a modeling framework robust to bias and noise.

In this work, we propose to jointly address these issues within a Bayesian framework for the spatio-temporal analysis of large-scale collections of multi-modal brain data. We show that this framework allows us to naturally encode plausibility constraints through clinically-inspired priors, while accounting for the uncertainty of the temporal profiles and brain structures we wish to estimate. Similarly to the ICA setting, we formulate the problem of trajectory modeling through matrix factorization across temporal and spatial sources. This is done for each modality by inferring their specific spatio-temporal sources. To promote smoothness in time and avoid any unnecessary hypothesis on the temporal trajectories, we rely on non-parametric modeling based on Gaussian Process (GP). We account for a plausible evolution from healthy to pathological stages thanks to a monotonicity constraint applied on the GP. Moreover, individuals’ observations are temporally re-aligned on a common scale via a time-warping function. In case of imaging data, to model the non-stationarity of the spatial signal, the spatial sources are defined as sparse activation maps convolved at different scales. We show that our framework can be efficiently optimized through stochastic variational inference, allowing to exploit automatic differentiation and GPU support to speed up computations.

The paper is organized as follows: Section 2 analyzes related work on spatio-temporal modeling of neurodegeneration, while Section 3 details our method. In Section 4 we present experiments on synthetic data in which we compare our model to standard blind source separation approaches. We finally provide a demonstration of our method on the modeling of imaging data from a large scale clinical study. Prospects for future work and conclusions are drawn in section 5. Derivations that we could not fit in the paper are detailed in the Appendices.

2 Related Work in Neurodegeneration Modeling

To deal with the uncertainty of the time-line of neurodegenerative pathologies, the concept of time-reparameterization of imaging-derived features has been used in several works. The underlying principle consists in estimating an absolute time-scale of disease progression by temporally re-aligning data from different subjects. For instance, in [42] the time-evolution was approximated as a sequence of events which need to be re-ordered for each patient. This approach thus considers the evolution of neurodegenerative diseases as a collection of transitions between discrete stages. This hypothesis is however limiting, as it doesn’t reflect the continuity of changes affecting the brain along the course of the pathology.

To address this limitation, we rely on a continuous parameterization of the time-axis as in [23, 10]. In particular, individuals’ observations are time-realigned on a common temporal scale via a time-warping function. Using a set of relevant scalar biomarkers, this kind of approach allows to learn a time-scale describing the pathology evolution, and to estimate a data-driven time-line markedly correlated with the decline of cognitive abilities. Similarly, in [4] a disease progression score was estimated using biomarkers from molecular imaging. These methods are however based on the analysis of low-dimensional measures, such as collections of clinical variables. Therefore, they do not allow to scale to the high dimension of multi-modal medical images. Our work tackles this shortcoming thanks to a scalable inference scheme based on stochastic variational inference.

Concerning the spatio-temporal representation of neurodegeneration, a mixed-effect model was proposed by [21] to learn an average spatio-temporal trajectory of brain evolution on cortical thickness data. The fixed-effect describes the average trajectory, while random effects are estimated through individual spatio-temporal warping functions, modeling how each subject differs from the global progression. Still, the extension of this approach to image volumes raises scalability issues. It has also to be noted that, to allow computational tractability, the brain evolution was assumed to be stationary both in space and time, thus limiting the ability of the model to disentangle the multiple dynamics of the brain structures involved in AD.

An attempt to source separation is proposed in [27], through the decomposition of cortical thickness measurements as a mixture of spatio-temporal processes. This is performed by associating to each cortical vertex a temporal progression modeled by a sigmoid function, which may be however too simplistic to describe the progression of AD temporal processes. We propose to overcome this issue by non-parametric modeling of the temporal sources through GPs. Moreover, the model in [27] is lacking of an explicit vertex-wise correlation model, as it only assumes correlation between clustering parameters at the resolution of the mesh graph. For this reason, it may still be sensitive to spatial variation at different scales and noise. We address this problem by modeling the spatial sources through convolution of sparse maps at multiple resolutions, allowing to deal with signal non-stationarity and robustness to noise.

3 Methods

In the following sections a matrix will be denoted by an uppercase letter $\boldsymbol{X}$, its n-th row will be given by $\boldsymbol{X}_{n:}$ and its n-th column by $\boldsymbol{X}_{:n}$. A column vector will be denoted by a lowercase letter $\boldsymbol{x}$. Subscript indices will be used to index the elements of matrices, vectors or sets of scalars. Superscipt indices will allow to index the blocks of block diagonal matrices.

3.1 Individual time-shift

To account for the uncertainty of the time-line of individual measurements, we assume that the observations are defined with respect to an absolute temporal reference $\tau$. This is performed through a time-warping function $t_{p}=\boldsymbol{f}_{p}(\tau)$, that models the individual time-reparameterization. We choose an additive parameterization such that:

[TABLE]

Within this setting the individual time-shift $\delta_{p}$ encodes the temporal position of subject $p$, which in our application can be interpreted as the disease stage of subject $p$ with respect to the long-term disease trajectory. We denote by $\boldsymbol{\delta}=\{\delta_{p}\}_{p=0}^{P}$ the set of time-shift parameters.

3.2 Data modeling

We represent the spatio-temporal data $\boldsymbol{D}$ by a block diagonal matrix in which we differentiate two main blocks $\boldsymbol{Y}$ and $\boldsymbol{V}$ as illustrated in Figure 1. Each sub-block $\boldsymbol{Y}^{m}$ is a matrix containing the data represented by one of the $M$ imaging modalities we wish to consider. These matrices have dimensions $P\times F_{m}$, where $P$ denotes the number of subjects and $F_{m}$ the number of imaging features for modality $m$, which in our case is the number of voxels. The matrix $\boldsymbol{V}$ accounts for non-imaging or scalar data such as clinical scores and has dimensions $P\times C$, where $C$ is the number of scalar features considered. We postulate a generative model and decompose the data as shown in Figure 1.

For each sub-block $\boldsymbol{Y}^{m}$, the data is factorized in a set of $N_{m}$ spatio-temporal sources $\boldsymbol{Y}^{m}=\boldsymbol{S}^{m}\boldsymbol{A}^{m}$. The columns of the matrix $\boldsymbol{S}^{m}$ describe the non-linear temporal evolution of the corresponding spatial maps contained in the rows of $\boldsymbol{A}^{m}$. Therefore, their product represents the voxel-wise linear combination of the spatial maps modulated by the corresponding temporal sources. The subjects share the same set of temporal sources across $\boldsymbol{S}^{1},..,\boldsymbol{S}^{M}$, as these sources describe the temporal evolution of the group-wise images through the regression problem specified in Figure 1. The data in matrix $\boldsymbol{V}$ is modelled by a matrix $\boldsymbol{U}$ whose columns depict the temporal trajectories of the different scalar scores. In the case of imaging data, we also consider a constant term modeling brain areas which don’t exhibit any intensity changes over time. This is done by including constant matrix terms $\boldsymbol{Z}^{m}$ that we need to estimate. We assume for a given modality $m$ that the vectors $\boldsymbol{Z}^{m}_{p:}$ are common to every subjects. Finally, for each modality $m$, scalar score $c$, and subject $p$, we assume Gaussian observational noise $\boldsymbol{\mathcal{E}}^{m}_{p:}\sim\mathcal{N}(\boldsymbol{0},\sigma^{2}_{m}\boldsymbol{I})$, and $\boldsymbol{\mathcal{H}}_{p,c}\sim\mathcal{N}(0,\nu^{2}_{c})$ for respectively imaging and scalar information.

Therefore, if we consider the data from modality $m$ and scalar $c$ of patient $p$ observed at time $\boldsymbol{f}_{p}(\tau)$ we have:

[TABLE]

We denote by $\theta_{m}$ and $\theta_{c}$ the temporal parameters related respectively to the modality m and scalar feature c, while $\psi_{m}$ represents the set of spatial parameters of modality m. We assume conditional independence across modalities and scalar scores given the time-shift information:

[TABLE]

Relying on classical regression formulation, we assume exchangeability across subjects allowing us to derive the data likelihood for a given modality m. According to the generative model we can write:

[TABLE]

Naturally, a similar equation holds for $p(\boldsymbol{V}_{:c}|\boldsymbol{U}_{:c},\boldsymbol{\delta},\nu_{c})$.

Within a Bayesian modeling framework, we wish to maximize the marginal log-likelihood $\log(p(\boldsymbol{Y},\boldsymbol{V}|\boldsymbol{Z},\boldsymbol{\delta},\sigma,\nu))$, to obtain posterior distributions for the spatio-temporal processes. Since the derivation of this quantity in a closed-form is not possible, we tackle this optimization problem through stochastic variational inference. Based on this formulation, in what follows we illustrate our model by detailing the variational approximations imposed on the spatio-temporal sources, along with the priors and constraints we impose to represent the data (Sections 3.3 and 3.4). Finally, we detail the variational lower bound and optimization strategy in Section 3.5.

For ease of notation we will drop the $m$ and $c$ indexes in Sections 3.3 and 3.4. As a result the matrix $\boldsymbol{S}$ will indistinctly refer to either any $\boldsymbol{S}^{m}$ or $\boldsymbol{U}$, while matrix $\boldsymbol{A}$ will refer to any $\boldsymbol{A}^{m}$, and $\boldsymbol{Y}$ to any $\boldsymbol{Y}^{m}$. For a given modality $m$, the number of patients $P$ will be indexed by $p$, the number of sources $N^{m}$ or the number of scalar scores $C$ will be indexed by $n$, and finally $f$ will index the number of imaging features $F^{m}$.

3.3 Spatio-temporal processes

3.3.1 Temporal sources

In order to flexibly account for non-linear temporal patterns, the temporal sources are encoded in a matrix $\boldsymbol{S}$ in which each column $\boldsymbol{S}_{:n}$ is a GP representing the evolution of source $n$ and is independent from the other sources. To allow computational tractability within a variational setting, we rely on the GP approximation proposed in [9], through kernel approximation via random feature expansion [32]. Within this framework, a GP can be approximated as a Bayesian Neural Network with form: $\boldsymbol{S}_{:n}(\boldsymbol{t})=\phi(\boldsymbol{t}(\boldsymbol{\omega}^{n})^{T})\boldsymbol{w}^{n}$. For example, in the case of the Radial Basis Function (RBF) covariance, $\boldsymbol{\omega}^{n}$ is a linear projection in the spectral domain. It is equipped with a Gaussian distributed prior $p(\boldsymbol{\omega}^{n})\sim\mathcal{N}(\boldsymbol{0},l_{n}\boldsymbol{I})$ with a zero-mean and a covariance parameterized by a scalar $l_{n}$, acting as the length-scale parameter of the RBF covariance. The non-linear basis functions activation is defined by setting $\phi(\cdot)=(\cos(\cdot),\sin(\cdot))$, while the regression parameter $\boldsymbol{w}^{n}$ is given with a standard normal prior. The GP inference problem can be conveniently performed by estimating approximated variational distributions for all the $\boldsymbol{\omega}^{n}$ and $\boldsymbol{w}^{n}$ (Section 3.5). We will respectively denote by $\boldsymbol{\Omega}$ and $\boldsymbol{W}$ the block diagonal matrices whose blocks are the $(\boldsymbol{\omega}^{n})^{T}$ and $\boldsymbol{w}^{n}$. Considering the $N$ temporal sources, we can write $p(\boldsymbol{\Omega})=\prod_{n}p(\boldsymbol{\omega}^{n})$ and $p(\boldsymbol{W})=\prod_{n}p(\boldsymbol{w}^{n})$.

We wish also to account for a steady evolution of the temporal processes, hence constraining the temporal sources to monotonicity. This is relevant in the medical case, where one would like to model the steady progression of a disease from normal to pathological stages. In our case, we want to constrain the space of the temporal sources to the set of solutions $\mathcal{C}_{n}=\{\boldsymbol{S}_{:n}(\boldsymbol{t})\mid\boldsymbol{S}_{:n}^{\prime}(\boldsymbol{t})\geq 0\quad\forall\ \boldsymbol{t}\}$. This can be done done consistently within the regression setting of [33], and in particular with the GP random feature expansion framework as shown in [22]. In that work, the constraint is introduced as a second likelihood term on the temporal sources dynamics:

[TABLE]

where $\boldsymbol{S}^{\prime}$ contains every derivatives $\boldsymbol{S}^{\prime}_{:n}$, $\gamma$ controls the magnitude of the monotonicity constraint, and $\mathcal{C}=\bigcap_{n}\mathcal{C}_{n}$. According to [22] this constraint can be specified through the parametric form for the derivative of each $\boldsymbol{S}_{:n}$:

[TABLE]

This setting leads to an efficient scheme for estimating the temporal sources through stochastic variational inference (Section 3.5).

3.3.2 Spatial sources.

According to the model introduced in Section 3.2, each observation $\boldsymbol{Y}_{p:}$ is obtained as the linear combination at a specific time-point between the temporal and spatial sources. In order to deal with the multi-scale nature of the imaging signal, we propose to represent the spatial sources at multiple resolutions. To this end, we encode the spatial sources in a matrix $\boldsymbol{A}$ whose rows $\boldsymbol{A}_{n:}$ represent a specific source at a given scale. The scale is prescribed by a convolution operator $\boldsymbol{\Sigma}^{n}$, which is a applied to a map $\boldsymbol{B}_{n:}$ that we wish to infer. This problem can be specified by defining $\boldsymbol{A}_{n:}=\boldsymbol{B}_{n:}\boldsymbol{\Sigma}^{n}$, where $\boldsymbol{\Sigma}^{n}$ is an $F\times F$ Gaussian kernel matrix imposing a specific spatial resolution. The length-scale parameter $\lambda_{n}$ of the Gaussian kernel is fixed for each source, to force the model to pick details at that specific scale. Due to the high-dimension of the data we are modeling, performing stochastic variational inference in this setting raises scalability issues. For instance, if we assume a Gaussian distribution $\mathcal{N}(\mu_{\boldsymbol{B}_{n:}},diag(\boldsymbol{\Lambda}))$ for $\boldsymbol{B}_{n:}$, the distribution of the spatial signal would be $p(\boldsymbol{A}_{n:})\sim\mathcal{N}(\boldsymbol{\mu}_{\boldsymbol{B}_{n:}}\boldsymbol{\Sigma}^{n},\boldsymbol{\Sigma}^{n}diag(\boldsymbol{\Lambda})(\boldsymbol{\Sigma}^{n})^{T})$. As a result, sampling from $p(\boldsymbol{A}_{n:})$ is not computationally tractable due to the size of the covariance matrix, which prevents the use of standard inference schemes on $\boldsymbol{B}_{n:}$. This can be overcome thanks to the separability of the Gaussian convolution kernel [28, 25], according to which the 3D convolution matrix $\boldsymbol{\Sigma}^{n}$ can be decomposed into the Kronecker product of 1D matrices, $\boldsymbol{\Sigma}^{n}=\boldsymbol{\Sigma}_{x}^{n}\otimes\boldsymbol{\Sigma}_{y}^{n}\otimes\boldsymbol{\Sigma}_{z}^{n}$. This decomposition allows to efficiently perform standard operations such as matrix inversion, or matrix-vector multiplication [35]. Thanks to this choice, we recover tractability for the inference of $\boldsymbol{B}_{n:}$ through sampling, as required by stochastic inference methods [20].

3.4 Sparsity

In order to detect specific brain areas involved in neurodegeneration, we propose to introduce a sparsity constraint on the maps (or codes) $\boldsymbol{B}_{n:}$. Consistently with our variational inference scheme, we induce sparsity via Variational Dropout as proposed in [19]. This approach leverages on an improper log-scale uniform prior $p(|\boldsymbol{B}_{n:}|)\propto\prod_{f}\ 1/|\boldsymbol{B}_{n,f}|$, along with an approximate posterior distribution:

[TABLE]

In this formulation, the dropout parameter $\alpha_{n,f}$ is related to the individual dropout probability $p_{n,f}$ of each weight by $\alpha_{n,f}=p_{n,f}(1-p_{n,f})^{-1}$. When the parameter $\alpha_{n,f}$ exceeds a fixed threshold, the dropout probability $p_{n,f}$ is considered high enough to ignore the corresponding weight $\boldsymbol{M}_{n,f}$ by setting it to zero. However, this framework raises stability issues affecting the inference of the dropout parameters due to large-variance gradients, thus limiting $p_{n,f}$ to values smaller than $0.5$. To tackle this problem, we leverage on the extension of Variational Dropout proposed in [29]. In this setting, the variance parameter is encoded in a new independent variable $\boldsymbol{P}_{n,f}=\alpha_{n,f}\boldsymbol{M}^{2}_{n,f}$, while the posterior distribution is optimized with respect to ($\boldsymbol{M},\boldsymbol{P}$). Therefore, in order to minimize the cost function for large variance $\boldsymbol{P}_{n,f}\to\infty$ ($\alpha_{n,f}\to\infty$ i.e $p_{n,f}\to 1$), the value of the weight’s magnitude must be controlled by setting to zero the corresponding parameter $\boldsymbol{M}_{n,f}$. As a result, by dropping out weights in the code, we sparsify the estimated spatial maps, thus better isolating relevant spatial sub-structures. Spatial correlations in the images are obtained thanks to the convolution operation detailed in Section 3.3.2.

3.5 Variational inference

We detailed in the previous sections the choices of priors and constraints that we apply to the spatio-temporal processes in order to plausibly model the data. To illustrate the overall formulation of the method, we provide in Figure 2 the graphical model over the M modalities in the case of imaging data. Naturally, this graph simplifies when we deal with scalar data as we don’t need to account for any spatial dependence.

To infer the time-shift parameter $\boldsymbol{\delta}$, the sets of parameters $\theta_{m}$, $\theta_{c}$, and $\psi_{m}$, as well as $\boldsymbol{Z}$, $\sigma$ and $\nu$, we need to jointly optimize the data evidence according to priors and constraints:

[TABLE]

We tackle the optimization of Equation (8) via stochastic variational inference. Following [9] and [22] we introduce approximations, $q_{2}(\boldsymbol{\Omega}^{m})$ and $q_{3}(\boldsymbol{W}^{m})$ in addition to $q_{1}(\boldsymbol{B}^{m})$ in order to derive a lower bound $\boldsymbol{\mathcal{L}}_{m}$ for each modality. We recall that the temporal trajectories $\boldsymbol{S}^{m}$ and $\boldsymbol{U}$ are treated similarly as described in Section 3.3.1. We also note that the choice of distributions $q_{1},q_{2}$ and $q_{3}$ is the same across modalities, while their parameters will be inferred independently. This leads to:

[TABLE]

Where $\mathcal{D}$ refers to the Kullback-Leibler (KL) divergence. Combining the lower bounds of the different modalities we obtain:

[TABLE]

A detailed derivation of the lower bound is given in Appendix A.

The approximated distributions $q_{2}(\boldsymbol{\Omega}^{m})$ and $q_{3}(\boldsymbol{W}^{m})$ are factorized across GPs such that:

[TABLE]

where $N_{rf}$ is the number of random features used for the projection in the spectral domain. Using Gaussian priors and approximations we introduced above, we can obtain a closed-form formula for the KL divergence. Moreover, the choice of prior and approximate posterior distribution for the maps of $\boldsymbol{B}^{m}$ leads to an approximation for the divergence $\mathcal{D}[q_{1}(\boldsymbol{B}^{m})||p(\boldsymbol{B}^{m})]$ detailed in [29]. This allows to analytically compute all the KL terms in our cost function. Formulas for the KL divergences are detailed in Appendix B.

Finally, we optimize the individual time-shifts $\boldsymbol{\delta}=\{\delta_{p}\}_{p=0}^{P}$, $\boldsymbol{Z}$, $\sigma=\{\sigma_{m}\}_{m=1}^{M}$, $\nu=\{\nu_{c}\}_{c=1}^{C}$ as well as the overall sets of spatio-temporal parameters $\boldsymbol{\theta}=\{\theta_{m}\}_{m=1}^{M}\cup\{\theta_{c}\}_{c=1}^{C}$ and $\boldsymbol{\psi}=\{\psi_{m}\}_{m=1}^{M}$.

[TABLE]

Following [20] and using the reparameterization trick, we can efficiently sample from the approximated distributions $q_{1},q_{2}$ and $q_{3}$ to compute the two expectation terms from (9) for each modality. We chose to alternate the optimization between the spatio-temporal parameters and the time-shift. We set $\gamma_{m}$ to the minimum value that gives monotonic sources. This was done through multiple tests on data batches with different numbers of imaging features $F_{m}$ and sources $N_{m}$. We empirically found that monotonicity was enforced when the magnitude of $\gamma_{m}$ was in the order of $F_{m}\times N_{m}$. The threshold for the dropout probability above which we set a weight $\boldsymbol{B}^{m}_{n,f}$ to zero was fixed at $95\%$ (i.e $\alpha=19$), while the $\sigma_{m}$ and $\nu_{m}$ were optimized during training along with the spatio-temporal parameters. The model is implemented and trained using the Pytorch library [30]. The complete experimental setting is detailed in Appendix C. We also provide a pseudo-code detailing the optimization procedure in Appendix D. In the following sections we will refer to our method as Monotonic Gaussian Process Analysis (MGPA).

4 Experiments and Results

In this section we first benchmark MGPA on synthetic data to demonstrate its reconstruction and separation properties while comparing it to standard sources separation methods. We finally apply our model on a large set of medical data from a publicly available clinical study, demonstrating the ability of our method to retrieve spatio-temporal processes relevant to AD, along with a time-scale describing the course of the disease.

4.1 Synthetic tests on spatio-temporal trajectory separation

For the synthetic tests we considered the case where the data is associated to a single imaging modality only. We tested MGPA on synthetic data generated as a linear combination of temporal functions and 3D activation maps at prescribed resolutions. The goal was to assess the method’s ability to identify the spatio-temporal sources underlying the data. We benchmarked our method with respect to ICA, Non-Negative Matrix Factorization (NMF), and Principal Component Analysis (PCA), which were applied from the standard implementation provided in the Scikit-Learn library [31].

The benchmark was specified by defining a $10$-folds validation setting, generating the data at each fold as a linear combination of temporal sources $\widetilde{\boldsymbol{S}}(\boldsymbol{t})=[\widetilde{\boldsymbol{S}}_{:0}(\boldsymbol{t}),\widetilde{\boldsymbol{S}}_{:1}(\boldsymbol{t})]$, and spatial maps $\widetilde{\boldsymbol{A}}=[\widetilde{\boldsymbol{A}}_{0:},\widetilde{\boldsymbol{A}}_{1:}]$. The data was defined as $\boldsymbol{Y}_{p:}=\widetilde{\boldsymbol{S}}_{p:}(\boldsymbol{t}_{p})\widetilde{\boldsymbol{A}}+\boldsymbol{\mathcal{E}}_{p:}$ over $50$ time points $\boldsymbol{t}_{p}$, where $\boldsymbol{t}_{p}$ was uniformly distributed in the range $[0,0.7]$, and $\boldsymbol{\mathcal{E}}_{p:}\sim\mathcal{N}(\boldsymbol{0},\sigma^{2}\boldsymbol{I})$. The temporal sources were specified as sigmoid functions $\widetilde{\boldsymbol{S}}_{p,i}(\boldsymbol{t}_{p})=1/(1+\exp(-\boldsymbol{t}_{p}+\alpha_{i}))$, while the spatial structures had dimensions $(30\times 30\times 30)$ such that $\widetilde{\boldsymbol{A}}_{i:}=\widetilde{\boldsymbol{B}}_{i:}\widetilde{\boldsymbol{\Sigma}}^{i}$. The $\widetilde{\boldsymbol{\Sigma}}^{i}$ were chosen as Gaussian convolution matrices with respective length-scale of $\lambda=2$ mm and $\lambda=1$ mm. The $\widetilde{\boldsymbol{B}}_{i:}$ were randomly sampled sparse 3D maps.

Variable selection. We applied our method by specifying an over-complete set of six sources with respective spatial length-scale of $\lambda=\{2,2,1,1,0.5,0.5\ mm\}$. Figure 3 shows an example of the sparse maps obtained for a specific fold. The model prunes the signal for most of the maps, while retaining two sparse maps, $\boldsymbol{B_{0:}}$ and $\boldsymbol{B_{4:}}$, whose length-scale are $\lambda=2$ mm and $\lambda=1$ mm, thus correctly estimating the right number of sources and their spatial resolution. As it can be qualitatively observed in Figure 3, we notice that the estimated sparse code convolved with a Gaussian kernel matrix with $\lambda=1$ mm is closer to its ground truth than the one convolved with a length-scale $\lambda=2$ mm. According to our tests, sparse codes associated to high resolution details (low $\lambda$) are indeed more identifiable. On the contrary, the identifiability of images obtained via a convolution operator with larger kernels (large $\lambda$) is lower, since these maps can be equivalently obtained through the convolution of different sparse codes.

Sources separation. We observe in Table 1 that the lowest Mean-Squared Error (MSE) for the temporal sources reconstruction is obtained by MGPA, closely followed by ICA. Similarly, our model and ICA show the highest Structural Similarity (SSIM) score [38], which quantifies the image reconstruction accuracy with respect to the ground truth maps, while accounting for the inter-dependencies between neighbouring pixels. An example of image reconstruction from a sample fold is illustrated in Figure 4. In this standard benchmark, we note that MGPA leads to comparable results with respect to the state of the art. In the following section, we compare the models in the more challenging setting in which the time-line has to be estimated as well.

4.2 Synthetic tests on trajectory separation and time-reparameterization

In this test, we modify the experimental benchmark by introducing a further element of variability associated to the time-axis. The temporal and spatial sources were modelled following the same procedure as in Section 4.1, however the observations were mixed along the temporal axis. To do so we generated longitudinal data as $\boldsymbol{Y}_{p,j,:}=\widetilde{\boldsymbol{S}}_{p:}(\boldsymbol{t})\widetilde{\boldsymbol{A}}+\boldsymbol{\mathcal{E}}_{j:}$, by sampling between $1$ and $10$ images per time-point and randomly re-arranging them along the time-axis (cf. time-shift $t_{p}$ of each observation at initialization in Figures 5 and 6, panel “Time-Shift”). The goal was to assess the sources separation performances of MGPA when the time-line is unknown. The experiment was run on 10 folds and Figures 5 and 6 illustrate the sources estimation for two different folds. We present these two figures to demonstrate how the time-shift inference affects the temporal sources reconstruction. Since the model is agnostic of a time-scale, we note that the time-shift may have a different range than the original time-axis. However, its relative ordering should be consistent with the original time points. We fitted a linear regression model over the 10 folds between the original time and the estimated time-shift parameter, and obtained an average R2 coefficient of $0.98$ with a standard deviation of $0.005$ (cf. Table 2).

This is illustrated for two different folds in the Time-Shift panel of Figures 5 and 6, where we observe a strong linear correlation with the original time-line, meaning that the algorithm correctly re-ordered the data with respect to the original time-axis. However, we notice in Table 2 that the MSE of the temporal sources significantly increased, due to the additional difficulty brought by the time-shift estimation. Indeed, in order to reconstruct the temporal signal we need to perfectly re-align hundreds of observations. This is the case in Figure 5 (optimal reconstruction result), where the time-shift is highly correlated with the original time-line, allowing to distinguish every single observation and reconstruct the original temporal profiles. Whereas in Figure 6 (sub-optimal reconstruction result), the estimated time-shift doesn’t exhibit a perfect fit, and generally underestimates the time-reparameterization for the later and earlier time points. This is related to the challenging setting of reconstructing the time-line identified by the original temporal sources. Indeed, we observe that $\boldsymbol{S}_{:0}$ reaches a plateau for early time points, while $\boldsymbol{S}_{:1}$ is flat for later ones. This behaviour increases the difficulty of differentiating time points with low signal differences. As a result, it impacts the time-shift optimization and adds variability to the time-shift estimation performances, thus deteriorating the reconstruction of the temporal sources over the 10 folds compared to the previous benchmark. The spatial sources estimation remains comparable to the one without time-shift both quantitatively, with an average SSIM of 95$\%$, and qualitatively, as shown in Figures 5 and 6. Within this setting, ICA, NMF and PCA poorly perform as they can’t reconstruct the time-line. Results obtained using these three methods are provided in Appendix E.

4.3 Application to spatio-temporal brain progression modeling

4.3.1 Data processing

Data used in the preparation of this article were obtained from the Alzheimer’s Disease Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu). The ADNI was launched in 2003 as a public-private partnership, led by Principal Investigator Michael W. Weiner, MD. For up-to-date information, see www.adni-info.org.

We selected a cohort of 544 amyloid positive subjects of the ADNI database composed of 103 controls (NL), 164 Mild Cognitive Impairment (MCI), 114 AD patients, 34 healthy individuals converted to MCI or to AD (NL converter) and 129 MCI converted to AD (MCI converter). The term amyloid positive refers to subjects whose amyloid level in the cerebrospinal fluid (CSF) is below the nominal cutoff of $192$ pg/ml. Conversion to MCI or AD was determined using the last follow-up available information. We provide in Table 3 socio-demographic and clinical information across the different groups.

MRI, FDG-PET and AV45-PET of each individual were processed in order to obtain respectively, volumes of gray matter density, glucose uptake, and amyloid load in a standard anatomical space.

MRI processing protocol. Baseline MRI images were analyzed according to the SPM12 processing pipeline [2]. Each image was initially segmented into grey, white matter and CSF probabilistic maps. Grey matter images were used for the following analysis, normalized to a group-wise reference space via DARTEL [1], and modulated using the Jacobian determinant of the subject-to-template transformation. The subsequent modeling was carried out on the normalised images at the original spatial resolution.

PET processing protocol. Individuals' baseline PET images were initially affinely aligned to the corresponding MRI. After scaling the intensities to the cerebellum, the images were normalized to the grey matter template obtained with DARTEL and smoothed with a FWHM parameter of 4.55.

Images have dimension $102\times 130\times 107$ before vectorization, leading to $1,418,820$ spatial features per patient. These spatial features represent for each voxel their gray matter concentration in the case of MRI images, their glucose metabolism for FDG-PET images, or their amyloid concentration for AV45-PET images. To exploit the ability of our model to automatically adapt to different spatial scales, we chose to keep the MRI images at their native resolution for the analysis, and thus do not perform additional smoohting to equalize to the PET FWHM. In addition to the imaging data of each patient, we also integrate the ADAS13 score assessed by clinicians. High values of this score indicate a decline of cognitive abilities. We consider three matrices $\boldsymbol{Y}^{MRI}$, $\boldsymbol{Y}^{FDG}$, and $\boldsymbol{Y}^{AV45}$ of dimension $(543\times 1,418,820)$ containing the images of all the subjects, and a matrix $\boldsymbol{V}$ of dimension $(543\times 1)$ containing their ADAS13 score. From now on we will refer to the data as the block diagonal matrix containing the four matrices $\boldsymbol{Y}^{MRI},\boldsymbol{Y}^{FDG}$, $\boldsymbol{Y}^{AV45}$, and $\boldsymbol{V}$ as described in Section 3.2. We note that the analysis is performed by only considering a single scan per imaging modality and ADAS13 score for each patient. Therefore, the temporal evolution has to be inferred solely through the analysis of relative differences between the brain morphologies, glucose metabolisms, amyloid concentrations and cognitive abilities across individuals.

4.3.2 Model specification

We aim at showing how MGPA applied on the data extracted from the ADNI cohort is able to temporally re-align patients in order to describe AD progression in a plausible way, while detecting relevant spatio-temporal processes at stake in AD. The model estimates AD progression by relying on MR, FDG-PET, AV45-PET scans and ADAS13 score of each patient. The temporal sources $\boldsymbol{S}^{MRI}$ and $\boldsymbol{S}^{FDG}$ associated respectively to the loss of gray matter, and to the decrease of glucose uptake, are enforced to be monotonically decreasing. On the contrary, the temporal sources $\boldsymbol{S}^{AV45}$ and $\boldsymbol{U}_{:ADAS13}$, modeling respectively the evolution of amyloid concentration, and ADAS13 score, are enforced to be monotonically increasing. Since we don’t consider any information about the disease stage of each individual before applying our method, all the observations are initialized at the same time reference $\tau=0$. Therefore, as for the tests in Section 4.2, the time-shift reparameterization describes a relative re-ordering of the subjects not related to a specific time-unit. To decompose the imaging data we apply our model by specifying an over-complete basis of six sources with $\lambda=\{8,8,4,4,2,2\ mm\}$, to cover both different scales and the associated variety of temporal evolution. Due to the high-dimension of the data matrix, the computations were parallelized over six GPUs, and the model required eighteen hours to complete the training. Details on the model convergence during training are provided in Appendix F.

4.3.3 Estimated spatio-temporal brain dynamics

In Figure 7 we show the spatio-temporal processes retained by the model for each imaging modality. Interestingly, the model adapts to the spatial resolution of MRI and PET images. Indeed, we notice that the model accounts for the high-resolution of MRI images by retaining a source associated to the lowest length-scale ($\lambda=2$ mm). Concerning PET data, we observe that the induced sparsity discards the highest resolution codes ($\lambda=2$ mm) for both FDG and AV45, highlighting the ability of the model to adapt to the coarser resolution of the PET signal.

In the case of MRI data, two sources were retained at two different resolutions ($\lambda=4$ mm and $\lambda=2$ mm). Source $\boldsymbol{S}^{MRI}_{4}$ describes gray matter loss encompassing a large extent of the brain with a focus on cortical areas (see $\boldsymbol{A}^{MRI}_{4}$). We note that this map also targets subcortical areas such as the hippocampi, which are key regions of AD. Source $\boldsymbol{S}^{MRI}_{2}$ ($\lambda=4$ mm) indicates a mild decrease of gray matter which accelerates in the latest stages of the disease, and targets the temporal poles (see $\boldsymbol{A}^{MRI}_{2}$). It is interesting to notice that this differential pattern of gray matter loss also affects the parahippocampal region, whose atrophy is known to be prominent in AD [11]. These results underline the complex evolution of brain atrophy, and the ability of the model to disentangle spatio-temporal processes mapping different regions involved in the pathology [3, 13]. Concerning the spatio-temporal processes extracted from the FDG-PET data, we see on Figure 7 that the model retained two sources at the coarsest resolutions ($\lambda=8$ mm). Source $\boldsymbol{S}^{FDG}_{1}$ indicates a pattern of hypometabolism that tends to plateau and which involves most of the brain regions, thus describing a global effect of the pathology on the glucose uptake. Source $\boldsymbol{S}^{FDG}_{0}$ describes a linear pattern of hypometabolism targeting areas such as the precuneus and the parietal lobe, which are known to be strongly affected during the evolution of the disease [5]. Finally, the model extracted two spatio-temporal sources from the AV45-PET data at two different resolutions ($\lambda=8$ mm and $\lambda=4$ mm). We observe that source $\boldsymbol{S}^{AV45}_{2}$ highlights an increase of amyloid deposition mapping a large extent of the brain, such as the parietal and frontal lobes as well as temporal areas, thus concurring with clinical evidence [34]. Similarly to the FDG-PET processes, we have a source $\boldsymbol{S}^{AV45}_{0}$ exhibiting a differential pattern of amyloid deposition targeting mostly frontal, temporal, occipital areas and precuneus.

The estimated spatio-temporal processes can be combined to obtain an estimated evolution $\boldsymbol{S}^{m}\boldsymbol{A}^{m}$ of the brain along the time-shift axis for each modality. In Figure 8, we show the ratio $|\boldsymbol{S}_{p:}^{m}\boldsymbol{A}^{m}-\boldsymbol{S}_{0:}^{m}\boldsymbol{A}^{m}|/\boldsymbol{S}_{0:}^{m}\boldsymbol{A}^{m}$ between the image predicted at four time-points $t_{p}$ and the image predicted at $t_{0}$ for the three imaging modalities. This allows us to visualize the trajectory of a brain going from a healthy to a pathological state in terms of atrophy, glucose metabolism and amyloid load according to our model.

Finally, we also applied ICA, NMF and PCA on the ADNI data, showing that the associated results are characterized by poor interpretability and high variability. The complete experimental setting and results are detailed in Appendix G.

4.3.4 Model Consistency

To verify the plausibility of the fitted model, we compare in Figure 9 the concentration predicted by the model and the raw concentration measures in different brain areas for the three imaging modalities. We observe a decrease of gray matter and glucose metabolism as we progress along the estimated time-line, allowing to relate large time-shift values to lower gray matter density and glucose uptake. Moreover, we notice the agreement between the predictions made by the model (in blue) and the raw concentration measures (in red). In the case of AV45 data there is only a mild increase of amyloid load according to the model, probably due to the fact that the subjects selected in the cohort are already amyloid positive. As a result, they already show a high baseline amyloid level concentration, close to plateau levels.

In Figure 10, we show the estimated GP $\boldsymbol{U}_{:ADAS13}$. We observe that the model is able to plausibly describe the evolution of this cognitive score, while demonstrating a larger variability than in the case of imaging modalities.

4.3.5 Plausibility with respect to clinical evidence

We assessed the clinical relevance of the estimated time-shift by relating it to independent medical information which were not included in the model during training. To this end, we compared the estimated time-shift to ADAS11, MMSE and FAQ scores. High values of ADAS11 and FAQ or low values of MMSE indicate a decline of performances. We show in Figure 11 that the estimated time-shift correlates with a decrease of cognitive and functional abilities. In particular, a cubic model slightly better describes the relationship between ADAS11 and the time-shift (according to BIC and AIC), with a significance for the cubic coefficient of $p=0.04$. Concerning MMSE and FAQ, quadratic and linear models were almost equivalent; the significance of the linear coefficients was $p<0.01$, while the quadratic coefficient was never significant. Pearson correlation coefficients for ADAS11, FAQ and MMSE were respectively of $0.49$, $0.41$, and $-0.45$, with corresponding p-values $p<0.01$.

The box-plot of Figure 12 shows the time-shift distribution across clinical groups. We observe an increase of the estimated time-shift when going from healthy to pathological stages. The high uncertainty associated to the MCI group is due to the broad definition of this clinical category, which includes subjects not necessarily affected by dementia. We note that MCI subjects subsequently converted to AD (MCI converter) exhibit higher time-shift than the clinically stable MCI group, highlighting the ability of the model to differentiate between conversion status. A similar distinction can be noticed between NL and NL converter groups. We found significant differences between median time-shift for NL-NL converter, MCI-MCI converter and MCI converter-AD (comparisons $p<0.01$, Figure 12). It is also important to recall that this result is obtained from the analysis of a single scan per imaging modality and ADAS13 score for each patient.

5 Discussion

We presented a generative approach to spatio-temporal disease progression modeling based on matrix factorization across temporal and spatial sources. The proposed application on a large set of medical images shows the ability of the model to disentangle relevant spatio-temporal processes at stake in AD, along with an estimated time-scale related to the disease evolution.

The model was compared to standard methods such as ICA, NMF and PCA since they perform blind source separation similarly to our method. This allowed us to demonstrate the advantages of building more complex approaches such as MGPA for the problem we tackle in this work. Concerning the comparison with the state of the art in disease progression modelling, to the best of our knowledge the two closest approaches are [27] and [21]. However, these two methods are specifically designed for modelling data defined on brain surfaces. On the contrary, our method aims at progression modeling using full 3D volumetric information. The data dimension we tackle is thus an order of magnitude greater than the one of [27] and [21], preventing these methods to scale to the spatial geometry of our data.

There are several avenues of improvement for the proposed approach. We found that the optimization is highly sensitive to the initialization of the spatial sources. This is typical of such complex non-convex problems, and requires further investigations to better control the algorithm convergence. More generally, the problem of source separation tackled in this work is intrinsically ill-posed, as the given data can be explained by several solutions. This was illustrated for example in our tests on synthetic data (Section 4.2), where the identification of the sources was more challenging in the case of coarse resolution codes and of flat temporal sources. We note however that this issue is general, and intrinsic to the problem of disease progression modeling.

Indeed, identifiability ultimately remains a critical issue when training the model. Concerning the spatio-temporal parameters, their number is extremely high due to the fact that we scale our method to 3D volumetric images. Estimating a single spatial source from a single modality requires to estimate the mean and variance of its sparse code, i.e $1,418,820\times 2=2,837,640$ parameters. In practice, hypotheses are explicitly introduced to reduce the number of effective parameters. For instance, the convolution of the spatial maps using Gaussian kernels allows to enforce smoothness, and thus reduces the number of effective degrees of freedom via spatial correlation across the related parameters. This is equivalent to the regularization applied to image registration problems, in which the number of parameters is of the same order of magnitude than in our setting. Moreover, our sparsity constraint allows to sensibly reduce the number of parameters at test time. Indeed, after training, the sparse codes of the MRI sources have $2,213,359$ non-zero elements instead of $17,025,840$, which amounts in $87\%$ reduction in the number of parameters. In the case of the FDG-PET and AV45-PET sparse codes, the number of non-zero elements at test time is respectively of $9,023,695$ and $1,362,067$, which is equivalent to a reduction in the number of parameters of $53\%$ and $92\%$. Nonetheless, this high number of parameters still remains a factor of potential convergence issues during the parameters estimation procedure. We present graphs in Appendix F showing the evolution of the different terms composing the cost function during training. These figures show convergence profiles typical of those obtained with stochastic variational inference schemes, such as with Variational Autoencoders or Bayesian Neural Networks. Moreover, the stability of the solution has been ensured through multiple runs of the model. Finally, as mentioned in Section 3.4, the Variational Dropout framework leads to stability issues affecting inference, which are mostly due to the use of an improper prior. This problem may motivate the identification of alternative ways to induce sparsity on the spatial maps.

In this work, we modeled the time-shift of each subject as a translation with respect to a common temporal reference. However, since pathological trajectories are different across individuals, it would be valuable to account for individual speed of progressions by introducing a scaling effect, as it has been proposed for example in [21, 36]. This was not in the scope of the current study, as we focused on the analysis of cross-sectional data, thus having only one data point per subject. Therefore, one of the main extensions of this model will be the integration of longitudinal data for each individual, which will allow a more specific time-reparameterization.

Our noise model for the reconstruction problem of Equation LABEL:eq:data_decomposition is homoscedastic and i.i.d. Gaussian with zero mean. For this reason, data variability for the entire image is encoded by the variance parameter of the Gaussian noise. Similarly as in standard regression problems, this modelling choice has been motivated to promote simplicity of the model and computational efficiency. However, around $40\%$ of the values in the brain images do not provide relevant information as they represent zero and constant background areas. For this reason, during training, the model can perfectly fit this background and increases its confidence on the overall regression solution, thus lowering the value of the noise variance $\sigma_{m}$ (cf Figure 9). This is in contrast to what we observe with the ADAS13 data (cf Figure 10), where the problem corresponds to standard univariate regression. A potential way to fix this issue could be to train the model only on non-zero image areas, or by implementing an heteroscedastic noise model. However, this latter solution may further increase the number of model parameters.

The modeling results are also sensitive to the specification of the spatio-temporal processes priors. In our case, the monotonicity constraint imposed to the GPs may be too restrictive to completely capture the complexity of the progression of neurodegeneration. From a clinical point of view, the model could also benefit from the integration of data measuring the concentration of Tau protein via PET imaging, in order to quantify key neurobiological processes associated to AD [17].

In order to guarantee that all the subjects belong to the same pathological trajectory due to AD, the model has only been applied to a cohort of amyloid positive subjects. However, this choice restricts the dynamics of evolution that we could estimate. Indeed, only considering these subjects narrows down the time-line of the pathology, as we study patients at potentially advanced disease stages. Therefore, it would be interesting in a future work to apply the model on a cohort including amyloid negative subjects, to model the brain dynamics over the whole disease natural history. This extension would require to define a proper methodology for disentangling sub-trajectories associated, for example with normal ageing and different pathological subtypes [24, 37, 40]. Moreover, we know that many patients diagnosed with AD can be associated to mixed pathologies such as vascular disease or Lewy bodies. Therefore, a potential clinical application of our method could be to investigate if the spatio-temporal dynamics estimated by MGPA are able to disentangle the contribution of each comorbidity.

Assessment of clinical plausibility of MGPA on the ADNI must be corroborated by further validation on independent datasets. Therefore, in a future work, we wish to validate the model on different cohorts to demonstrate its generalization properties. The validation step for each subject would be done by estimating the time-point minimizing the cost between the images of each tested individual, and the image progression model previously estimated on ADNI. The estimated time-shift would provide a measure of the pathological stage of the individual with respect to the modelled trajectory, and could be then compared with the clinical diagnosis of the subject, allowing to test the reliability of our model. This additional validation step could ultimately allow to use the model as a diagnostic instrument of AD. This validation would require an important effort in terms of data harmonisation across multiple cohorts, as well as in terms of clinical interpretation. For this reason, this work will be part of a subsequent publication.

We planned to release the source-code along with instructions in order for the model to be used by a large audience. It will be available as a complementary tool on the platform http://gpprogressionmodel.inria.fr/, which already offers a simple front-end to Gaussian Process Progression model.

6 Acknowledgements

This work has been supported by the French government, through the UCAJEDI Investments in the Future project managed by the National Research Agency (ref.n ANR-15-IDEX-01), the grant AAP Santé 06 2017-260 DGA-DSH, and by the Inria Sophia Antipolis - Méditerranée, ”NEF” computation cluster.

Data collection and sharing for this project was funded by the Alzheimer’s Disease Neuroimaging Initiative (ADNI) and DOD ADNI. ADNI is funded by the National Institute on Aging, the National Institute of Biomedical Imaging and Bioengineering, and through generous contributions from the following: AbbVie, Alzheimer’s Association; Alzheimer’s Drug Discovery Foundation; Araclon Biotech; BioClinica, Inc.; Biogen; Bristol-Myers Squibb Company;CereSpir, Inc.;Cogstate;Eisai Inc.; Elan Pharmaceuticals, Inc.; Eli Lilly and Company; EuroImmun; F. Hoffmann-La Roche Ltd and its affiliated company Genentech, Inc.; Fujirebio; GE Healthcare; IXICO Ltd.; Janssen Alzheimer Immunotherapy Research & Development, LLC.; Johnson & Johnson Pharmaceutical Research & Development LLC.;Lumosity;Lundbeck;Merck & Co., Inc.; Meso Scale Diagnostics, LLC.;NeuroRx Research; Neurotrack Technologies;Novartis Pharmaceuticals Corporation; Pfizer Inc.; Piramal Imaging;Servier; Takeda Pharmaceutical Company; and Transition Therapeutics.The Canadian Institutes of Health Research is providing funds to support ADNI clinical sites in Canada. Private sector contributions are facilitated by the Foundation for the National Institutes of Health (www.fnih.org). The grantee organization is the Northern California Institute for Research and Education, and the study is coordinated by the Alzheimer’s Therapeutic Research Institute at the University of Southern California. ADNI data are disseminated by the Laboratory for Neuro Imaging at the University of Southern California.

Appendix A.

In this Appendix, we detail the complete derivation of the lower bound.

[TABLE]

By observing that $\frac{d\boldsymbol{S}^{m}}{d\boldsymbol{t}}$ is completely identified by $\boldsymbol{S}^{m}$, the equation can be written as:

[TABLE]

Similarly this derivation can be applied to $\log(p(\boldsymbol{V}_{:c},\mathcal{C}^{c}|\boldsymbol{\delta},\nu_{c},\gamma_{c}))$.

[TABLE]

This derivation gives us the lower bound $\boldsymbol{\mathcal{L}}_{m}$ of a given modality m. The same technique can be used to derive a lower bound for $\log(p(\boldsymbol{V}_{c:},\mathcal{C}^{c}|\boldsymbol{\delta},\nu_{c},\gamma_{c}))$, and by summation over $m$ and $c$ we obtain the lower bound of Equation 10 for $\log(p(\boldsymbol{Y},\boldsymbol{V},\mathcal{C}|\boldsymbol{Z},\boldsymbol{\delta},\sigma,\nu,\gamma))$.

Appendix B.

In this section we provide formulas for computing the three KL terms of the lower bound. The total KL divergences are:

[TABLE]

For ease of notation we will drop the $m$ and $c$ indices and will give formulas for a single modality. In [29], authors provide an approximation of the KL for the maps $\boldsymbol{B}$:

[TABLE]

where $h$ is the sigmoid function and $k_{1}=0.63576,\ k_{2}=1.87320,\ k_{3}=1.48695$.

In the case of $\boldsymbol{\Omega}$ and $\boldsymbol{W}$, we’ve seen that they have Gaussian priors and approximations which are detailed in Sections 3.3.1 and 3.5. As a result we can obtain closed-form formulas for their KL, leading to:

[TABLE]

By summation over the different modalities we finally obtain the total KL divergences.

Appendix C.

We provide in this Appendix details for the experiments on real data.

The number of random features for the GP estimation was set to $10$, as it was enough to recover the temporal sources in the synthetic experiments.

• 2.

The $\gamma$ parameter controlling monotonicity was set to $\gamma_{m}=10^{7}$ for each imaging modality ($F_{m}=1,418,820$ imaging features and $N_{m}$ = 6 sources) and $\gamma_{c}=1$ for ADAS13 ($C_{c}=1$ scalar feature).

• 3.

The lower bound was optimized using the ADAM optimizer [18].

• 4.

We used an alternate optimization scheme between the spatio-temporal parameters and the time-shift of [2000, 1000] iterations repeated 20 times, followed by 30000 iterations in which we only optimized the spatio-temporal parameters.

• 5.

The expectation terms in the lower bound were approximated using only one Monte-Carlo sample as proposed in [20].

• 6.

The table below gives the learning rates (LR) of all the parameters of the model.

Appendix D.

In this Appendix, we first provide a pseudo-code for sampling from a normal distribution using the reparameterization trick (see Algorithm 1). The second pseudo-code (Algorithm 2) details the steps to compute the lower bound $\boldsymbol{\mathcal{L}}_{m}$ for a given imaging modality m. We recall that we want to optimize the following sets of parameters (see Section 3.5): $\boldsymbol{\delta}=\{\delta_{p}\}_{p=0}^{P}$, $\boldsymbol{Z}$, $\sigma=\{\sigma_{m}\}_{m=1}^{M}$, $\nu=\{\nu_{c}\}_{c=1}^{C}$, $\boldsymbol{\theta}=\{\theta_{m}\}_{m=1}^{M}\cup\{\theta_{c}\}_{c=1}^{C}$, and $\boldsymbol{\psi}=\{\psi_{m}\}_{m=1}^{M}$. Where $P$ is the number of subjects, $M$ the number of imaging modalities, $C$ the number of scalar features, and $N_{m}$ the number of spatio-temporal sources for a given modality m.

[TABLE]

Similarly to Algorithm 2, we can derive a function LOSS_SCALAR when dealing with scalar scores by removing the computations on the spatial sources. Finally the last pseudo-code (Algorithm 3) details the model optimization. For sake of clarity we denote by $\boldsymbol{\Pi}$, the set of all the spatio-temporal parameters of the model.

Appendix E.

In this Appendix, we show results obtained with standard methods (ICA, NMF, PCA) when applied within the experimental setting of Section 4.2. We recall that for these experiments observations were randomly aligned along the time-axis. The goal was to assess the ability of the different methods to reconstruct the spatio-temporal sources underlying the data when the time-axis is unknown. Results obtained in Table 1 show a substantial decrease of performances for the MSE and SSIM compared to MGPA (cf Table 2 in Section 4.2). Indeed, these methods do not consider time as a variable on which inference is required, thus preventing them from reconstructing correctly the temporal sources. Figure 1 shows an example of reconstruction when using ICA. We observe that even though the spatial reconstruction remains acceptable, the estimated temporal sources are not interpretable as ICA reconstructs the data using the time-axis on which observations have been mixed.

Appendix F.

We provide in this Appendix details on the model convergence when applied on the ADNI data. The training was divided in three iterations of 30000 epochs each. During the two first iterations the spatio-temporal parameters and the time-shift are trained alternatively following a scheme of [2000,1000] epochs ten times. The third iteration only optimizes the spatio-temporal parameters. In Figure 1, we show the evolution of the total loss and the different terms composing it during training. The term reconstruction cost stands for $\sum_{m}\operatorname{\mathbb{E}}_{q_{1},q_{2},q_{3}}[\log(p(\boldsymbol{Y}^{m}|\boldsymbol{B}^{m},\boldsymbol{\Omega}^{m},\boldsymbol{W}^{m},\boldsymbol{Z}^{m},\boldsymbol{\delta},\sigma_{m}))]$, monotonicity cost for $\sum_{m}\operatorname{\mathbb{E}}_{q_{2},q_{3}}[\log(p(\mathcal{C}^{m}|\boldsymbol{\Omega}^{m},\boldsymbol{W}^{m},\boldsymbol{\delta},\gamma_{m}))]$ and KL for $\sum_{m}\mathcal{D}[q_{1}(\boldsymbol{B}^{m})||p(\boldsymbol{B}^{m})]+\mathcal{D}[q_{2}(\boldsymbol{\Omega}^{m})||p(\boldsymbol{\Omega}^{m})]+\mathcal{D}[q_{3}(\boldsymbol{W}^{m})||p(\boldsymbol{W}^{m})]$. We observe that through the first two iterations the reconstruction and monotonicity costs decrease, and become stable during the last iteration. Differently, the KL cost increases during the first iteration as the model is driven by the reconstruction and monotonicity constraints. The KL term decreases during the second iteration, thus regularizing the model, before becoming stable during the third iteration. We also note that the graphs in Figure 1 show convergence profiles typical of those obtained with stochastic variational inference schemes, such as with Variational Autoencoders or Bayesian Neural Networks.

Appendix G.

In this Appendix, we provide the results obtained when applying ICA, NMF and PCA on the ADNI data of Section 4.3.1. We used the three imaging modalities for each subject and concatenated these images in a $(544\times 4256460)$ matrix. Our goal was to compare the spatio-temporal processes extracted using these standard methods with the ones from MGPA. We recall that in the case of MGPA the model automatically re-aligns the observations following monotonic assumptions for each biomarker, while these standard methods don’t perform any inference on the time variable. Therefore, we created three experimental settings in which we changed the observations’ alignment. In the first one, subjects were aligned by their chronological age (Figures 1, 2 and 3), in the second one by ADAS13 (Figures 4, 5 and 6) and in the last one time was randomly initialized like in the experiments of Section 4.3.3 (Figures 7, 8 and 9). We extracted six spatio-temporal sources for each method and each time-alignment, like in 4.3.2.

We observe that the temporal profiles are generally noisy and hard to interpret due to the lack of constraints on the temporal evolution. This motivates the need of smooth and monotonic constraints as in MGPA. Moreover, due to the concatenation of all the modalities they all share the same temporal patterns. This is an important difference with the modality-specific modelling of MGPA. Finally, we note that the spatial patterns associated with each method are very similar, independently from the time-initialization, while the temporal sources substantially differ. This is also true when time is randomly initialized. These observations point to the challenge of giving a clinical interpretation of the results obtained with these approaches, and therefore to the need of plausible spatio-temporal constraints as provided in MGPA.

Subjects aligned by age.

Subjects aligned by ADAS13.

Subjects randomly aligned.