DIVE: A spatiotemporal progression model of brain pathology in neurodegenerative disorders

TL;DR

DIVE is a novel image-based model that reconstructs detailed brain pathology progression patterns from short-term data, enabling better disease understanding and patient stratification.

Contribution

It introduces a vertexwise disease progression model that clusters cortical measurements and estimates individual disease stages, providing new insights into neurodegenerative disorders.

Findings

Successfully recovers ground truth clusters in simulated data.

Identifies consistent atrophy patterns in Alzheimer's across datasets.

Distinguishes pathology patterns between different neurodegenerative diseases.

Abstract

Here we present DIVE: Data-driven Inference of Vertexwise Evolution. DIVE is an image-based disease progression model with single-vertex resolution, designed to reconstruct long-term patterns of brain pathology from short-term longitudinal data sets. DIVE clusters vertex-wise biomarker measurements on the cortical surface that have similar temporal dynamics across a patient population, and concurrently estimates an average trajectory of vertex measurements in each cluster. DIVE uniquely outputs a parcellation of the cortex into areas with common progression patterns, leading to a new signature for individual diseases. DIVE further estimates the disease stage and progression speed for every visit of every subject, potentially enhancing stratification for clinical trials or management. On simulated data, DIVE can recover ground truth clusters and their underlying trajectory, provided the average trajectories are sufficiently different between clusters. We demonstrate DIVE on data from two cohorts: the Alzheimer's Disease Neuroimaging Initiative (ADNI) and the Dementia Research Centre (DRC), UK, containing patients with Posterior Cortical Atrophy (PCA) as well as typical Alzheimer's disease (tAD). DIVE finds similar spatial patterns of atrophy for tAD subjects in the two independent datasets (ADNI and DRC), and further reveals distinct patterns of pathology in different diseases (tAD vs PCA) and for distinct types of biomarker data: cortical thickness from Magnetic Resonance Imaging (MRI) vs amyloid load from Positron Emission Tomography (PET). Finally, DIVE can be used to estimate a fine-grained spatial distribution of pathology in the brain using any kind of voxelwise or vertexwise measures including Jacobian compression maps, fractional anisotropy (FA) maps from diffusion imaging or other PET measures. DIVE source code is available online: https://github.com/mrazvan22/dive

Figures

Tables

Table 1. Table 1 : Comparison of our model with two more simplistic models on the ADNI MRI dataset. For each of the three models, we show the correlation of the disease progression scores (DPS) with respect to several cognitive tests: CDRSOB, ADAS13, MMSE and RAVLT. The correlation numbers represent the mean correlation across the 10 cross-validation folds.

Model	CDRSOB ($\rho$)	ADAS13 ($\rho$)	MMSE ($\rho$)	RAVLT ($\rho$)	Prediction (RMSE)
DIVE	0.37 +/- 0.09	0.37 +/- 0.10	0.36 +/- 0.11	0.32 +/- 0.12	1.021 +/- 0.008
ROI-based model	0.36 +/- 0.10	0.35 +/- 0.11	0.34 +/- 0.13	0.30 +/- 0.13	1.019 +/- 0.010
No-staging model	*0.09 +/- 0.06	*0.03 +/- 0.09	*0.05 +/- 0.06	*0.02 +/- 0.06	*1.062 +/- 0.024

Equations

$$M^{(u)}=\underset{M}{\mathrm{argmax}}\sum _{z1,\dots ,z_L}^{K}p(Z=(z_1,...,z_L)|V,M^{(u-1)})\left[logp(V,Z|M)\right]+logp(M)$$

$$p(Z|V,M)=\frac{1}{C}p(V,Z|M)=\prod _l^L\left[\prod _{(i,j)\in I}N(V_l^{ij}|f({\alpha }_i{t}_{ij}+{\beta }_i|{\theta }_{Z_l}),{\sigma }_{Z_l})\prod _{l_2\in N_l}\Psi (Z_l,Z_{l_2})\right]$$

$$p(Z|V,M)\approx \prod _l^L{\mathbb{E}}_{Z_{N_l}^{(u-1)}|V_l,M}\left[p(Z_l|V_l,M,Z_{N_l}^{(u-1)})\right]$$

$$p(Z_l|V_l,M)\approx \frac{1}{C}\sum _{Z_{N_l}^{(u-1)}}p(V_l|Z_l,M)p(Z_l|Z_{N_l}^{(u-1)})p(Z_{N_l}^{(u-1)}|V_l,M)$$

$$p(Z_l|V_l,M)\approx \frac{1}{C}p(V_l|Z_l,M)\sum _{z_{N_l(1)},..,z_{N_l(|N_l|)}}\prod _{m\in \{1,...,N_l\}}p(Z_l|Z_{N_l(m)}^{(u-1)}=z_{N_l(m)})p(Z_{N_l(m)}^{(u-1)}=z_{N_l(m)}|V_l,M)$$

$$p(Z_l|V_l,M)=p(V_l|Z_l,M)\prod _{l_2\in N_l}\sum _{z_{l_2}}p(Z_l|Z_{l_2}^{(u-1)}=z_{l_2})p(Z_{l_2}^{(u-1)}=z_{l_2}|V_l,M)$$

$$p(Z_l|V_l,M)=p(V_l|Z_l,M)\prod _{l_2\in N_l}\sum _{k_2}p(Z_l|Z_{l_2}^{(u-1)}=k_2)p(Z_{l_2}^{(u-1)}=k_2|V_l,M)$$

$$z_{lk}^{(u)}\propto \left[\prod _{i,j\in I}N(V_l^{ij}|f({\alpha }_i{t}_{ij}+{\beta }_i;{\theta }_k),{\sigma }_k)\right]\left[\prod _{l_2\in N_l}\sum _{k_2=1}^K{z}_{l_2{k}_2}^{(u-1)}\Psi (Z_l=k,Z_{l_2}=k_2)\right]$$

$$log{z}_{lk}^{(u)}\propto \left[\sum _{i,j\in I}-\frac{log(2\pi {\sigma }_k^2)}{2}-\frac{1}{2{\sigma }_k^2}(V_l^{ij}-f({\alpha }_i{t}_{ij}+{\beta }_i;{\theta }_k){)}^2\right]+\left[\sum _{l_2\in N_l}log\sum _{k_2=1}^K{z}_{l_2{k}_2}^{(u-1)}({\delta }_{k_{2}k}exp(\lambda )+(1-{\delta }_{k_{2}k})exp(-{\lambda }^2))\right]$$

$$D_{lk}=-\frac{log(2\pi {\sigma }_k^2)|I|}{2}-\sum _{i,j\in I}\frac{1}{2{\sigma }_k^2}(V_l^{ij}-f({\alpha }_i{t}_{ij}+{\beta }_i;{\theta }_k){)}^2$$

$$log{z}_{lk}^{(u)}\propto D_{lk}+\left[\sum _{l_2\in N_l}log\sum _{k_2=1}^K{z}_{l_2{k}_2}^{(u-1)}({\delta }_{k_{2}k}(exp(\lambda )-exp(-{\lambda }^2))+exp(-{\lambda }^2))\right]$$

$$\begin{array}{r}log{z}_{lk}^{(u)}\propto D_{lk}+\left[\sum _{l_2\in N_l}log\left[exp(-{\lambda }^2)+z_{l_{2}k}^{(u-1)}(exp(\lambda )-exp(-{\lambda }^2))\right]\right]\end{array}$$

$$z_{lk}=\frac{e^{x(k)}}{e^{x(1)}+e^{x(2)}+\cdots +e^{x(K)}}=\frac{1}{e^{x(1)-x(k)}+e^{x(2)-x(k)}+\cdots +e^{x(K)-x(k)}}$$

$$\theta =\underset{\theta }{\mathrm{argmax}}\sum _{z1,\dots ,z_L}^{K}p(Z=(z_1,...,z_L)|V,M^{(u-1)})\left[logp(V,Z|M)\right]+logp(\theta )$$

$${\theta }_k=\underset{{\theta }_k}{\mathrm{argmax}}\sum _{z1,\dots ,z_L}^{K}p(Z=(z_1,...,z_L)|V,M^{(u-1)})\left[logp(V,Z|M)\right]+logp({\theta }_k)$$

$${\theta }_k=\underset{{\theta }_k}{\mathrm{argmax}}\sum _{z1,\dots ,z_L}^{K}p(Z=(z_1,...,z_L)|V,M^{(u-1)})log\left[\prod _{l=1}^L\prod _{(i,j)\in I}N(V_l^{ij}|f({\alpha }_i{t}_{ij}+{\beta }_i;{\theta }_{z_l}),{\sigma }_{z_l})\right]+logp({\theta }_k)$$

$${\theta }_k=\underset{{\theta }_k}{\mathrm{argmax}}\sum _{z1,\dots ,z_L}^{K}p(Z=(z_1,...,z_L)|V,M^{(u-1)})\sum _{l=1}^L\sum _{(i,j)\in I}logN(V_l^{ij}|f({\alpha }_i{t}_{ij}+{\beta }_i;{\theta }_{z_l}),{\sigma }_{z_l})+logp({\theta }_k)$$

$${\theta }_k=\underset{{\theta }_k}{\mathrm{argmax}}\sum _{l=1}^L\sum _{z_l=1}^{K}p(Z_l=z_l|V,M^{(u-1)})\sum _{(i,j)\in I}logN(V_l^{ij}|f({\alpha }_i{t}_{ij}+{\beta }_i;{\theta }_{z_l}),{\sigma }_{z_l})+logp({\theta }_k)$$

$${\theta }_k=\underset{{\theta }_k}{\mathrm{argmax}}\sum _{l=1}^{L}p(Z_l=k|V,M^{(u-1)})\sum _{(i,j)\in I}logN(V_l^{ij}|f({\alpha }_i{t}_{ij}+{\beta }_i;{\theta }_k),{\sigma }_k)+logp({\theta }_k)$$

$${\theta }_k=\underset{{\theta }_k}{\mathrm{argmax}}\sum _{l=1}^{L}p(Z_l=k|V,M^{(u-1)})\sum _{(i,j)\in I}\left[log(2\pi {\sigma }_k{)}^{-1/2}-\frac{1}{2{\sigma }_k^2}(V_l^{ij}-f({\alpha }_i{t}_{ij}+{\beta }_i;{\theta }_k){)}^2\right]+logp({\theta }_k)$$

$${\theta }_k=\underset{{\theta }_k}{\mathrm{argmax}}\sum _{l=1}^{L}p(Z_l=k|V,M^{(u-1)})\sum _{(i,j)\in I}\left[-\frac{1}{2{\sigma }_k^2}(V_l^{ij}-f({\alpha }_i{t}_{ij}+{\beta }_i;{\theta }_k){)}^2\right]+logp({\theta }_k)$$

$${\sigma }_k=\underset{{\sigma }_k}{\mathrm{argmax}}\sum _{l=1}^{L}p(Z_l=k|V,M^{(u-1)})\sum _{(i,j)\in I}logN(V_l^{ij}|f({\alpha }_i{t}_{ij}+{\beta }_i;{\theta }_k),{\sigma }_k)$$

$${\sigma }_k=\underset{{\sigma }_k}{\mathrm{argmax}}\sum _{l=1}^{L}p(Z_l=k|V,M^{(u-1)})\sum _{(i,j)\in I}\left[log(2\pi {\sigma }_k^2{)}^{-1/2}-\frac{1}{2{\sigma }_k^2}(V_l^{ij}-f({\alpha }_i{t}_{ij}+{\beta }_i;{\theta }_k){)}^2\right]$$

$$\frac{\delta l({\sigma }_k|.)}{\delta {\sigma }_k}=\sum _{l=1}^{L}p(Z_l=k|V,M^{(u-1)})\sum _{(i,j)\in I}\frac{\delta }{\delta {\sigma }_k}\left[log(2\pi {\sigma }_k^2{)}^{-1/2}-\frac{1}{2{\sigma }_k^2}(V_l^{ij}-f({\alpha }_i{t}_{ij}+{\beta }_i;{\theta }_k){)}^2\right]$$

$$\frac{\delta l({\sigma }_k|.)}{\delta {\sigma }_k}=\sum _{l=1}^{L}p(Z_l=k|V,M^{(u-1)})\sum _{(i,j)\in I}\left[-\frac{\delta }{\delta {\sigma }_k}\frac{log{\sigma }_k^2}{2}-\frac{\delta }{\delta {\sigma }_k}\frac{1}{2{\sigma }_k^2}(V_l^{ij}-f({\alpha }_i{t}_{ij}+{\beta }_i;{\theta }_k){)}^2\right]$$

$$\frac{\delta l({\sigma }_k|.)}{\delta {\sigma }_k}=\sum _{l=1}^{L}p(Z_l=k|V,M^{(u-1)})\sum _{(i,j)\in I}\left[-\frac{1}{{\sigma }_k}-\frac{-2}{2{\sigma }_k^3}(V_l^{ij}-f({\alpha }_i{t}_{ij}+{\beta }_i;{\theta }_k){)}^2\right]$$

$$\frac{\delta l({\sigma }_k|.)}{\delta {\sigma }_k}=\sum _{l=1}^{L}p(Z_l=k|V,M^{(u-1)})\sum _{(i,j)\in I}\left[-\frac{{\sigma }_k^2}{{\sigma }_k^3}+\frac{1}{{\sigma }_k^3}(V_l^{ij}-f({\alpha }_i{t}_{ij}+{\beta }_i;{\theta }_k){)}^2\right]$$

$$\frac{\delta l({\sigma }_k|.)}{\delta {\sigma }_k}=\sum _{l=1}^{L}p(Z_l=k|V,M^{(u-1)})\sum _{(i,j)\in I}\left[-{\sigma }_k^2+(V_l^{ij}-f({\alpha }_i{t}_{ij}+{\beta }_i;{\theta }_k){)}^2\right]=0$$

$${\sigma }_k^2=\frac{1}{|I|}\sum _{l=1}^{L}p(Z_l=k|V,M^{(u-1)})\sum _{(i,j)\in I}(V_l^{ij}-f({\alpha }_i{t}_{ij}+{\beta }_i;{\theta }_k){)}^2$$

$${\alpha }_i,{\beta }_i=\underset{{\alpha }_i,{\beta }_i}{\mathrm{argmax}}\sum _{l=1}^L\sum _{k=1}^{K}p(Z_l=k|V,M^{(u-1)})\sum _{(i^{\prime },j)\in I}logN(V_l^{i^{\prime }j}|f({\alpha }_{i^{\prime }}{t}_{i^{\prime }j}+{\beta }_{i^{\prime }};{\theta }_k),{\sigma }_k)+logp({\alpha }_i,{\beta }_i)$$

Full text

See pages - of final_manuscript.pdf

Supplementary Material

1 Simulations - Error in Estimated Trajectories and DPS

2 Comparison between DIVE and other models

2.1 Motivation

We were also interested to compare the performance of DIVE with other disease progression models. In particular, we were interested to test whether:

• •

Modelling dynamic clusters on the brain surface improves subject staging and biomarker prediction

• •

Modelling subject-specific stages with a linear transformation (the $\alpha_{i}$ and $\beta_{i}$ terms) improves biomarker prediction

2.2 Experiment design

We compared the performance of our model to two simplified models:

• •

ROI-based model: groups vertices according to an a-priori defined ROI atlas. This model is equivalent to the model by Jedynak et al., Neuroimage, 2012 and is a special case of our model, where the latent variables $z_{lk}$ are fixed instead of being marginalised as in equation 6.

• •

No-staging model: This is a model that doesn’t perform any time-shift of patients along the disease progression timeline. It fixes $\alpha_{i}=1$, $\beta_{i}=0$ for every subject, which means that the disease progression score of every subject is age.

We performed this comparison using 10-fold cross-validation. For each subject in the test set, we computed their DPS score and correlated all the DPS values with the same four cognitive tests used previously. We also tested how well the models can predict the future vertex-wise measurements as follows: for every subject i in the test set, we used their first two scans to estimate $\alpha_{i}=1$, $\beta_{i}=0$ and then used the rest of the scans to compute the prediction error. For one vertex location on the cortical surface, the prediction error was computed as the root mean squared error (RMSE) between its predicted measure and the actual measure. This was then averaged across all subjects and visits.

2.3 Results

Table 1 shows the results of the model comparison, on ADNI MRI dataset. Each row represents one model tested, while each column represents a different performance measure: correlations with four different cognitive tests and accuracy in the prediction of future vertexwise measurements. In each entry, we give the mean and standard deviation of the correlation coefficients or RMSE across the 10 cross-validation folds.

3 Validation of subject parameters against APOE

4 Derivation of the Generalised EM algorithm

We seek to calculate $M^{(u)}=\operatorname*{arg\,max}_{M}\mathbb{E}_{Z|V,M^{(u-1)}}\left[log\ p(V,Z|M)\right]+log\ p(M)$ where $M^{(u)}=(\alpha^{(u)},\beta^{(u)},\theta^{(u)},\sigma^{(u)},\lambda^{(u)})$ are the set of model parameters at iteration $u$ of the EM algorithm, and $Z=(Z_{1},\dots,Z_{L})$ is the set of discrete latent variables, where $Z_{l}$ represents the cluster that voxel $l$ was assigned to, so $Z_{l}\in{1,\dots,K}$. While $Z_{l}$ (with capital letter) is a random variable, we will also use the notation $z_{l}$ (small letter) to represent the value that the random variable $Z_{l}$ was instantiated to. Finally, $p(M^{(u)})$ is a prior on these parameters that is chosen by the user. Expanding the expected value, we get:

[TABLE]

The E-step involves computing $p\left(Z=(z_{1},...,z_{L})|V,M^{(u-1)}\right)$, while the M-step comprises of solving the above equation.

4.1 E-step

In this step we need to estimate $p(Z|V,M^{(u-1)})$. For notational simplificy we will drop the $(u-1)$ superscript from $M$

[TABLE]

where $N_{l}$ is the set of neighbours of vertex $l$. However, this doesn’t directly factorise over the vertices $l$ due to the MRF terms $\Psi(Z_{l},Z_{l_{2}})$. It is however necessary to find a form that factorises over the vertices, otherwise we won’t be able to represent in memory the joint distribution over all $Z$ variables. If we make the approximation $p(Z|V,M)\approx\prod_{l}^{L}p(V_{l}|Z_{l},M)$ then we loose out all the MRF terms and the model won’t account for spatial correlation. We instead do a first-degree approximation by conditioning on the values of $Z_{N_{l}}^{(u-1)}$, the labels of nearby vertices from the previous iteration. The approximation is thus:

[TABLE]

This form allows us to factorise over all the vertices to get $p(Z_{l}|V_{l},M)$:

[TABLE]

where $C$ is a normalistion constant that can be dropped. We can now further factorise $p(Z_{l}|Z_{N_{l}}^{(u-1)})\approx\prod_{m\in\{1,...,N_{l}\}}p(Z_{l}|M,Z_{N_{l}(m)}^{(u-1)}=z_{N_{l}(m)})$ and apply a similar factorisation to the prior $p(Z_{N_{l}}^{(u-1)}|V_{l},M)$, resulting in:

[TABLE]

Factorising the summation over $z_{N_{l}}$’s we get:

[TABLE]

Replacing $z_{l_{2}}$ with $k_{2}$ we get:

[TABLE]

We shall also denote $z_{lk}=p(Z_{l}|V_{l},M)$. Further simplifications result in:

[TABLE]

[TABLE]

We further define the data-fit term $D_{lk}$ as follows:

[TABLE]

This results in:

[TABLE]

Finally, we simplify the sum over $k_{2}$ to get the update equation for $z_{lk}$:

[TABLE]

In practice, we cannot naively compute the exponential term $z_{lk}=exp(log(z_{lk}))$ due to precision loss. However, we go around this by recomputing the exponentiation and normalisation of $z_{lk}$ simultaneously. Denoting $x(k)=log\ z_{lk}$, for $k\in[1\dots K]$, we get:

[TABLE]

4.2 M-step

The M-step itself does not have a closed-form analytical solution. We choose to solve it by successive refinements of the cluster trajectory parameters and the subject time shifts.

4.2.1 Optimising trajectory parameters

**Trajectory shape - $\theta$

**

Taking equation 1 and fixing the subject time-shifts $\alpha$, $\beta$ and measurement noise $\sigma$, we can find its maximum with respect to $\theta$ only. More precisely, we want:

[TABLE]

We observe that for each cluster the individual $\theta_{k}$’s are conditionally independent, i.e. $\theta_{k}\perp\!\!\!\perp\theta_{m}|\{Z,\alpha,\beta,\sigma$} $\forall k,m$. We also assume that the prior factorizes for each $\theta_{k}$: $\log p(\theta)=\prod_{k}^{K}\log p(\theta_{k})$. This allows us to optimise each $\theta_{k}$ independently:

[TABLE]

Replacing the full data log-likelihood, we get:

[TABLE]

Note that we didn’t include the MRF clique terms, since they are not a function of $\theta_{k}$. We propagate the logarithm inside the products:

[TABLE]

We next assume that $Z_{l}$, the hidden cluster assignment for vertex $l$, is conditionally independent of the other vertex assignments $Z_{m}$, $\forall m\neq l$ (See E-step approximation from Eq. 3). This independence assumption induces the following factorization: ${p(Z=(z_{1},...,z_{L})|V,M^{(u-1)})=\prod_{l}^{L}p(Z_{l}=z_{l}|V,M^{(u-1)})}$. Propagating this product inside the sum over the vertices, we get:

[TABLE]

The terms which don’t contain $\theta_{k}$ dissapear:

[TABLE]

We further expand the gaussian noise model:

[TABLE]

Constants dissapear due to the $\operatorname*{arg\,max}$ and we get the final update equation for $\theta_{k}$:

[TABLE]

**Measurement noise - $\sigma$

**

We first assume a uniform prior on the $\sigma$ parameters to simplify derivations. Using a similar approach as with $\theta$, after propagating the product inside the logarithm and removing the terms which don’t contain $\sigma_{k}$, we get:

[TABLE]

Note that, just as for $\theta$ above, the MRF clique terms were not included because they are not a function of $\sigma_{k}$. Expanding the noise model we get:

[TABLE]

The maximum of a function $l(\sigma_{k})$ can be computed by taking the derivative of the function $l$ and setting it to zero. This is under the assumption that $l$ is differentiable, which it is but we won’t prove it here. This gives:

[TABLE]

Propagating the differential operator further inside the sums we get:

[TABLE]

We next perform several small manipulations to reach a more suitable form of the derivative and then set it to be equal to zero:

[TABLE]

[TABLE]

[TABLE]

Finally, we solve for $\sigma_{k}$ and get its update equation:

[TABLE]

4.2.2 Estimating subject time shifts - $\alpha$, $\beta$

For estimating $\alpha$, $\beta$, we adopt a similar strategy as in the case of $\theta$, up to Eq. 18. This gives us the following problem:

[TABLE]

The terms $\alpha_{i^{\prime}},\beta_{i^{\prime}}$ for other subjects $i^{\prime}\neq i$ dissappear:

[TABLE]

Expanding the gaussian noise model we get:

[TABLE]

After removing constant terms we end up with the final update equation for $\alpha_{i}$, $\beta_{i}$:

[TABLE]

4.2.3 Estimating MRF clique term - $\lambda$

We optimise $\lambda$ using the following formula:

[TABLE]

Note that $p(Z|V,M^{(u-1)},\lambda,Z^{(u-1)})$ is a function of $\lambda$, so for each lambda we estimate $z_{lk}$ through approximate inference. We do this because otherwise the optimisation of $\lambda$ will only take into account the clique terms and completely exclude the data terms. We further simplify the objective function for lambda as follows:

[TABLE]

We take the logarithm:

[TABLE]

Let us denote $z_{lk}=p(Z_{l}=k|V,M^{(u-1)},\lambda,Z^{(u-1)})$. Assuming independence between the latent variables $Z_{l}$ we get:

[TABLE]

However, we now want to make $z_{lk}$ a function of $\lambda$ as previously mentioned, so $z_{lk}=\zeta_{lk}(\lambda)$, for some function $\zeta_{lk}$. More precisely, using the E-step update from Eq. 12 we define for each vertex $l$ and cluster $k$ a function $\zeta_{lk}(\lambda)$ as follows:

[TABLE]

where $D_{lk}$ is as defined in Eq 10. We then replace $z_{lk}$ with $\zeta_{lk}(\lambda)$ and introduce the chosen MRF clique model to get:

[TABLE]

We separate the cliques that have matching clusters to the ones that don’t:

[TABLE]

We also factorise the clique terms:

[TABLE]

Finally, we simplify to get the objective function for $\lambda$.

[TABLE]

For implementation speed-up, data-fit terms $D_{lk}$ can be pre-computed.

5 Fast DIVE Implementation - Proof of Equivalence

Fitting DIVE can be computationally prohibitive, especially given that the number of vertices/voxels can be very high, e.g. more than 160,000 in our datasets. We derived a fast implementtion of DIVE, which is based on the idea that for each subject we compute a weighted mean of the vertices within a particular cluster, and then compare that mean with the corresponding trajectory value. This is in contrast with comparing the value at each vertex with the corresponding trajectory of its cluster. In the next few sections, we will present the mathematical formulation of the fast implementation for parameters [$\theta$, $\alpha$, $\beta$]. Parameter $\sigma$ already has a closed-form update, while parameter $\lambda$ has a more complex update procedure for which this fast implementation doesn’t work. For each parameter, we will also provide proofs of equivalence.

5.1 Trajectory parameters - $\theta$

5.1.1 Fast implementation

The fast implementation for $\theta$ implies that, instead of optimising Eq. 29 we optimise the following problem:

[TABLE]

where $<V^{ij}>_{\hat{Z}_{k}}$ is the mean value of the vertices belonging to cluster $k$. Mathematically, we define $\hat{Z}_{k}=[z_{1k}\gamma_{k},\ z_{2k}\gamma_{k},\ \dots,\ z_{Lk}\gamma_{k}]$ where $\gamma_{k}=(\sum_{l=1}^{L}z_{lk})^{-1}$ is the normalisation constant. Moreover, we have that ${<V^{ij}>_{\hat{Z}_{k}}=\sum_{l=1}^{L}z_{lk}\gamma_{k}V^{ij}}$. We take the derivative of the likelihood function $l_{fast}$ of the fast implementation (Eq. 36) with respect to $\theta_{k}$ and perform several simplifications:

[TABLE]

[TABLE]

using the fact that $\sum_{l=1}^{L}\gamma_{k}z_{lk}=1$ we get:

[TABLE]

[TABLE]

By setting the derivative to zero, the optimal $\theta$ is thus a solution of the following equation:

[TABLE]

5.1.2 Slow implementation

We will prove that if theta is a solution of the slow implementation, it is also a solution of Eq. 41, which will prove that the fast implementation is equivalent. The slow implementation is finding $\theta$ from the following equation:

[TABLE]

Taking the derivative of the function above ($l_{slow}$) with respect to $\theta_{k}$ we get:

[TABLE]

After swapping terms around and using distributivity we get:

[TABLE]

This is the same optimisation problem as in Eq. 41, which proves that the two formulations are equivalent with respect to $\theta$.

5.2 Noise parameter - $\sigma$

The noise parameter $\sigma$ can actually be computed in a closed-form solution for the original slow model implementation, so there is no benefit in implementing the fast update for $\sigma$. Moreover, the $\sigma$ in the fast implementation computed the standard deviation in the mean value of the vertices within a certain cluster, and not the deviation withing the actual value of the vertices.

5.3 Subjects-specific time shifts - $\alpha$, $\beta$

5.3.1 Fast implementation

The equivalent fast formulation for the subject-specific time shifts is similar to the one for the trajectory parameters. It should be noted however that we need to weight the sums corresponding to each cluster by $\gamma_{k}^{-1}$. This gives the following equation for the fast formulation:

[TABLE]

In order to prove that this is equivalent to the slow version, we need to take the derivative of the likelihood function ($l_{fast}$) from the above equation with respect to $\alpha_{i}$, $\beta_{i}$ and set it to zero:

[TABLE]

We expand the average across the vertices and slide the derivative operator inside the sums:

[TABLE]

Since $\sum_{l=1}^{L}\gamma_{k}z_{lk}=1$ we get:

[TABLE]

Removing the factor 2 and sliding $\gamma_{k}$:

[TABLE]

Further sliding $\sum_{l=1}^{L}z_{lk}$ to the left we get the final optimisation problem:

[TABLE]

5.3.2 Slow implementation

In a similar way to the trajectory parameters, we want to prove that solving the problem from Eq. 50 (fast implementation) is the same as solving the original slow implementation problem, which is defined as:

[TABLE]

Taking the derivative of the function above with respect to $\alpha_{i},\beta_{i}$, we get:

[TABLE]

This is the same problem as the fast implementation one from Eq. 50, thus the fast model is equivalent to the slow model with respect to $\alpha$, $\beta$.