There is a great need to understand the progression of Alzheimer’s Disease (AD) especially before the clinical symptoms to better target therapeutic interventions [hampel2017]
. During this silent phase, neuroimaging reveals the disease effects on brain structure and function, such as the atrophy of the cortex due to neuronal loss. However, the precise dynamics of the lesions in the brain are not so clear at the group level and even less at the individual level. Personalized models of lesion propagation would enable to relate structural or metabolic alterations to the clinical signs, offering ways to estimate stage of the disease progression in the pre-symptomatic phase. Numerical models have been introduced to describe the temporal and the spatial evolution of these alterations, defining aspatiotemporal trajectory of the disease, i.e. a description of the changes in the brain over time, such as lesion progressions, tissue deformation and atrophy propagation.
Statistical models are well suited to estimate distributions of spatiotemporal patterns of propagation out of series of short-term longitudinal observations [donohue2014, bilgel2016]. However, the absence of time correspondence between patients is a clear obstacle for these types of approaches. Using data series of several individuals requires to re-align the series of observations in a common time-line and to adjust to a standardized pace of progression. Current models either consider a sequential propagation [young2015], without taking into account the continuous dynamics of changes, or develop average scenarios [guerrero2016, iturria2015]. Recently, a generic approach to align patients has been proposed for a set of biomarkers in [schiratti2015nips]: the temporal inter-subject variability results from individual variations of a common time-line granting each patient a unique age at onset and pace of progression. On top of the time-alignment of the observations, there exists a spatial variability of the signal propagation that characterizes a distribution of trajectories.
In order to exhibit a spatial representation of the alterations, we study medical images or image-derived features taking the form of a signal discretized at the vertices of a mesh, for instance the cortical thickness distributed on the mesh of the pial surface or Standardized Uptake Value Ratio (SUVR) distributed on the regular voxel grid of a PET scan. The spatial distribution of the signal is encoded in a distance matrix, giving the physical distance between the graph nodes. A sensible prior to include in the model is to enforce smooth variations of the temporal profile of signal changes across neighbouring nodes, highlighting a propagation pattern across the network as in [raj2012]. Extending directly the model in [schiratti2015nips] may lead to an explosion of the number of parameters proportional to the mesh resolution. At infinite resolution, the parameters take the form of a smooth continuous map defined on the image domain. In this paper, we propose to constrain these maps to belong to a finite-dimensional Hilbert Space, penalizing high frequency variations. In practice, these maps are generated by the convolution of parameter values at a sparse set of control nodes on the network with a smoothing kernel. The number of control nodes, whose distribution is determined by the bandwidth of the kernel, controls the complexity of the model regardless of the mesh resolution. Furthermore, the propagation of non-normalized signal could not adequately be modeled by the same curve shifted in time as in [schiratti2015nips]. We introduce new parameters to account for smooth changes in the profiles of changes at neighbouring spatial locations.
We introduce a mixed-effect generative model that learns a distribution of spatiotemporal trajectory from series of repeated observations. The model evaluates individual parameters (time reparametrization and spatial shifts) that enables the reconstruction of individual disease propagation through time. This non-linear problem is tackled by a stochastic version of the EM algorithm, the MCMC-SAEM [kuhn2005, allassonniere2010] in a high-dimensional setting. It considers fixed-effects describing a group-average trajectory and random effects characterizing individual trajectories as adjustment of the mean scenario. It is used to detect the cortical thickness variations in MRI data of MCI converters from the ADNI database.
2 Manifold-valued networks
In the following, we consider a longitudinal dataset of individuals, such that the th individual is observed at repeated time points . We assume that each observation takes the form of scalar measures referred to as a measurement map.
2.0.1 Manifold-valued measurements distributed on a fixed graph
Let be a non-oriented graph where is a set of vertices of a mesh in and is a subset of pairs of vertices defining the graph edges. We assume that is a common fixed graph such that the th coordinate of each measurement map corresponds to the vertex . As the graph corresponds to measurements spatially distributed on a mesh, the edges embed a spatial configuration. Therefore, any edge is valued with , a geodesic distance on the graph, defining a distance matrix such that for all , . Each measurement map produces a network , i.e. a fixed graph with one-dimensional values associated to each vertex and with distances associated to each edge.
We consider that the measurements of the patients at each node corresponds to observations of a signal function at particular time points. Thus, the function describes the evolution of the signal over the whole network. We assume that each signal function is continuous and that measurement map lies in a space defined by smooth constraints, as expected for bounded or normalized observations (eg. volume ratios, thickness measures, SUVR). Therefore, the space of measurements is best described as a Riemannian manifold [docarmo1992, lee2003], leading to consider each function as a one-dimensional geodesically complete Riemannian manifold such that all spatial observation is a point in the product manifold . It follows that for each , is a manifold-valued network.
2.0.2 Spatial smoothness of the propagation
Besides the temporal smoothness of the propagation, we expect the signal to be similar for neighbour nodes. We consider that each node is described by
parameters that parametrize the signal trajectory. In order to ensure smooth variations of the parameters values at neighbouring nodes, we assume that they result from the interpolation of the parameter values at a sparse subset of uniformly distributed nodes, called control nodes. For each parameter , potentially estimated at each node, the control nodes define a parameter evaluation function encoding for all the nodes: , , and where the are the new model parameters and is a Gaussian Kernel such that , being the geodesic distance on the graph and the kernel bandwidth.
This convolution guarantees the spatial regularity of the signal propagation. Moreover this smooth spatial constraint enables a reduction of the number of parameters, reducing the dimensional complexity from independent parameters at each node, to parameters only at the control nodes.
3 The statistical model
3.0.1 A propagation model
Given a set of manifold-valued networks , the model describes a group-average trajectory in the space of measurements, defined by a geodesic that allows to estimate a typical scenario of progression. Individual trajectories derive from the group-average scenario through spatiotemporal transformations: the exp-parallelization and the time reparametrization.
First, to describe disease pace and onset specific to each subject, we introduced a temporal transformation, called the time-warp, that is defined, for the subject , by where is the reference time-point in the space of measurements. The parameter corresponds to the time-shift between the mean and the individual age at onset and is the acceleration factor that describes the pace of an individual, being faster or slower than the average. This time reparametrization allows to reposition the dynamics of the average scenario in the real time-line of the th individual.
allows to translation the observations in the space of measurements, from the mean scenario to individual trajectories, encoding a variation in the trajectory of changes across the nodes of the graph. This exp-parallelization is handled by a family of individual vectors, called space-shifts. As shown on Figure 1 (left), the orange dots refer to individual observations in the space of measurements. The group-average trajectory estimated from the longitudinal measurements corresponds to the blue line. The space shifts characterize a spatial shift perpendicular to that describes the velocity of the mean scenario.
Finally, the parameters allow the reconstruction of the individual trajectories from the mean scenario of propagation.
Given a noise , the mixed-effect model writes, for a arbitrary vertex function :
3.0.2 Parameters estimation with the MCMC-SAEM algorithm
To reconstruct the long-term scenario of the disease propagation, we estimate the parameters of the group-average trajectory using a maximum likelihood estimator. The random-effects are considered as latent variables, whose distributions characterize the variability of the individual trajectories. Due to the non-linearity in Eq. (1), we use a Stochastic Approximation Expectation Maximization [delyon1999]
coupled with a Monte-Carlo Markov Chain sampler (MCMC-SAEM)[kuhn2004]. Let be the current estimation of the parameters and the current iterate of the Markov chain of the latent variables. The algorithm alternates between a simulation step, a stochastic approximation step and a maximization step, until convergence [allassonniere2010]. The simulation uses an adaptive version [atchade2006] of the Hasting Metropolis within Gibbs sampler to draw from . This algorithm was chosen as it handles non-linear mixed effects models [kuhn2005] with proven convergence and consistent estimations in practice.
3.0.3 Model instantiation
As many measurements correspond to positive values (eg. the cortical thickness, volume ratios), we consider in the following the open interval as a one-dimensional Riemannian manifold equipped with a Riemannian metric such that for all and for all , . With this metric and given , is a geodesically complete Riemannian manifold whose geodesics are of the form where , . These parameters are represented on Figure 1 (right) at two nodes where the decrease of the signal varies spatially. For identifiability reasons, we choose to fix the parameters among the nodes, leading to a shared parameter such that for all . As can be arbitrarily chosen in , we fix defined in Section 3.0.1. Considering the interpolation functions introduced in Section 2.0.2 and the fact that the parameters are , it leads to define and
Finally, the model defined in (1) rewrites:
such that and
4 Experimental results
We used this model to highlight typical spatiotemporal patterns of cortical atrophy during the course of Alzheimer’s Disease from longitudinal MRI of 154 MCI converters from the ADNI database, which amounts for 787 observations, each subject being observed 5 times on average. We aligned the measures on a common atlas with FreeSurfer [reuter2012] so that the measurement maps are distributed on the same common fixed-graph which is constituted of 1827 nodes that map the surface of the left pial surface. Out of the vertices, we selected 258 control nodes uniformly distributed over the surface. They encode the spatial interpolation of the propagation. The distance matrix is defined by a geodesic distance on .
4.0.2 Cortical thickness measurements
We used the model instantiation defined in Section 3.0.3 to characterize the cortical thickness decrease. Multiple runs of 30.000 iterations (
4hours) of this MCMC-SAEM lead to a noise standard deviationmm with of the data included in mm. The mean spatiotemporal propagation, described on the first three rows of the Figure 2 as the cortical thickness at respectively 65, 68, 71 and 74 years old shows that the primarily affected area is the medial-temporal lobe, followed by the temporal neocortex. The parietal association cortex and the frontal lobe are also subject to important alterations. On the other side, the sensory-motor cortex and the visual cortex are less involved in the lesion propagation. These results are consistent with previous knowledge of the Alzheimer’s Disease effects on the brain structure. As the model is able to exhibit individual spatiotemporal patterns with their associated pace of progression, the fourth and fifth rows of the Figure 2 represent consecutively the effect of the parallel shifting and of the time reparametrization on the cortical thickness atrophy. The figure 2(a) shows the real cortical thickness of a subject and the reconstruction predicted by the model. The relative error and its histogram are represented on Figure 2(b). The reconstruction is coherent to the real observation, the remaining error represents essentially an unstructured noise that we precisely try to smooth out.
5 Discussion and perspectives
We proposed a mixed-effect model which is able to evaluate a group-average spatiotemporal propagation of a signal at the nodes of a mesh thanks to longitudinal neuroimaging data distributed on a common network. The network vertices describe the evolution of the signal whereas its edges encode a distance between the nodes via a distance matrix. The high dimensionality of the problem is tackled by the introduction of control nodes: they allow to evaluate a smaller number of parameters while ensuring the smoothness of the signal propagation through neighbour nodes. Moreover, individual parameters characterize personalized patterns of propagation as variations of the mean scenario. The evaluation of this non-linear high dimensional model is made with the MCMC-SAEM algorithm that leads to convincing results as we were able to highlight areas affected by considerable neuronal loss: the medial-temporal lobe or the temporal neocortex.
The distance matrix, which encodes here the geodesic distance on the cortical mesh, may be changed to account for the structural or functional connectivity information. In this case, signal changes may propagate not only across neighbouring locations, but also at nodes far apart in space but close to each other in the connectome. The model can be used with multimodal data, such as PET scans, introducing numerical models of neurodegenerative diseases that could inform about the disease evolution at a population level while being customizable to fit individual data, predicting stage of the disease or time to symptom onset.
This work has been partly funded by ERC grant No678304, H2020 EU grant No666992, and ANR grant ANR-10-IAIHU-06.