A model of brain morphological changes related to aging and Alzheimer's disease from cross-sectional assessments

TL;DR

This paper introduces a deformation-based generative model that jointly captures brain morphological changes due to aging and Alzheimer's disease, providing interpretable biomarkers for disease progression and aging effects.

Contribution

It presents a novel spatio-temporal framework that models normal aging and Alzheimer's disease progression simultaneously using MRI data, generating biomarkers and morphological trajectories.

Findings

Markers differentiate clinical conditions despite high variability.

Model reveals aging effects around ventricles and AD effects in temporal lobe and hippocampus.

Markers indicate accelerated aging in Alzheimer's patients.

Abstract

In this study we propose a deformation-based framework to jointly model the influence of aging and Alzheimer's disease (AD) on the brain morphological evolution. Our approach combines a spatio-temporal description of both processes into a generative model. A reference morphology is deformed along specific trajectories to match subject specific morphologies. It is used to define two imaging progression markers: 1) a morphological age and 2) a disease score. These markers can be computed locally in any brain region. The approach is evaluated on brain structural magnetic resonance images (MRI) from the ADNI database. The generative model is first estimated on a control population, then, for each subject, the markers are computed for each acquisition. The longitudinal evolution of these markers is then studied in relation with the clinical diagnosis of the subjects and used to generate possible morphological evolution. In the model, the morphological changes associated with normal aging are mainly found around the ventricles, while the Alzheimer's disease specific changes are more located in the temporal lobe and the hippocampal area. The statistical analysis of these markers highlights differences between clinical conditions even though the inter-subject variability is quiet high. In this context, the model can be used to generate plausible morphological trajectories associated with the disease. Our method gives two interpretable scalar imaging biomarkers assessing the effects of aging and disease on brain morphology at the individual and population level. These markers confirm an acceleration of apparent aging for Alzheimer's subjects and can help discriminate clinical conditions even in prodromal stages. More generally, the joint modeling of normal and pathological evolutions shows promising results to describe age-related brain diseases over long time scales.

Figures

Tables

Table 1. Table 1: Socio-demographic and clinical information of the study cohort. Standard deviations are shown in parentheses.

group	CN-	CN+	MCIs-	MCIs+	MCIc	AD
Number of subjects	108	69	96	120	228	203
Age at baseline	73.4 (5.6)	74.5 (6.5)	71.1 (7.7)	73.5 (6.6)	73.8 (7.1)	74.5 (7.7)
Gender (female)	47.2%	56.5%	47.9%	37.5%	41.7%	48.3%
Education (years)	16.4 (2.6)	16.2 (2.7)	16.2 (2.8)	16.4 (2.7)	16.0 (2.9)	15.0 (2.9)
ADAS13 at baseline	9.0 (4.0)	8.6 (5.0)	12.0 (4.9)	14.0 (5.4)	19.9 (6.7)	31.4 (7.3)

Table 2. Table 2: Specified mean regional atrophy for the healthy and the pathological evolutions. The goal is to get simple but realistic atrophy patterns. It should be noted that the atrophy is specified using the divergence of the SVF in the area. The local volume changes are computed using a spatio-temporal integration scheme.

brain area	mean pathological (in %)	mean healthy (in %)
white matter	1.0	0.8
cortex	3.0	0.4
hippocampi	5.2	1.0
amygdalae	5.2	1.0
entorhinal cortex	6.5	0.7
temporal poles	6.2	0.6
other areas	0.0	0.0

Equations

$$I_k=T_0\circ {\varphi }_k+{ϵ}_k.$$

$$\mathcal{T}=\left\{T_0\circ \Phi ({\lambda }_{MA},{\lambda }_{DS})\text{ for }{\lambda }_{MA},{\lambda }_{DS}\in \mathbb{R}\right\}.$$

$$I_k=T_0\circ \Phi ({\lambda }_{MA}^k,{\lambda }_{DS}^k)\circ {\varphi }_k+{ϵ}_k,$$

$$\varphi (x)={\varphi }_1(x)={\int }_0^{1}v({\varphi }_t(x))dt\text{with }{\varphi }_0=\text{Id}.$$

$$\mathcal{T}=\left\{T_0\circ \mathrm{Exp}({\lambda }_{MA}{\mathrm{v}}_A+{\lambda }_{DS}{\mathrm{v}}_D)\text{ for }{\lambda }_{MA},{\lambda }_{DS}\in \mathbb{R}\right\}.$$

$$I_k=T_0\circ \mathrm{Exp}(({\lambda }_{MA}^k-{\lambda }_0){\mathrm{v}}_A+{\lambda }_{DS}^k{\mathrm{v}}_D)\circ {\varphi }_k+{ϵ}_k,$$

$$\mathrm{Exp}(({\lambda }_{MA}^k-{\lambda }_0){\mathrm{v}}_A+{\lambda }_{DS}^k{\mathrm{v}}_D)\circ {\varphi }_k=\mathrm{Exp}(({\lambda }_{MA}^k-{\lambda }_0){\mathrm{v}}_A+{\lambda }_{DS}^k{\mathrm{v}}_D+w_r^k).$$

$$w^k=({\lambda }_{MA}^k-{\lambda }_0){\mathrm{v}}_A+{\lambda }_{DS}^k{\mathrm{v}}_D+w_r^k.$$

$$v^k=s_{MA}^k{\mathrm{v}}_A+s_{DS}^k{\mathrm{v}}_D+v_r^k.$$

$$\underset{{\mathrm{v}}_A,{\mathrm{v}}_D}{\min }\sum _{k\in {\mathcal{G}}_h}\Vert v^k-{\mathrm{v}}_A{\Vert }^2+\sum _{k\in {\mathcal{G}}_{ad}}\Vert v^k-{\mathrm{v}}_A-{\mathrm{v}}_D{\Vert }^2.$$

$${\hat{\mathrm{v}}}_A$$

$${\hat{\mathrm{v}}}_D$$

$$\underset{w_t\in \mathcal{T},w_r\in {\mathcal{T}}^{\perp },w^{\prime }}{\min }\text{Sim}(I,T,\mathrm{Exp}(w_t+w_r))+\text{Dist}(w_r,w^{\prime })+\text{Reg}(w^{\prime }).$$

$$w^k\cdot {\mathrm{v}}_A$$

$$w^k\cdot {\mathrm{v}}_D$$

$$E(\Vert {\hat{\mathrm{v}}}_A{\Vert }^2)$$

$$=\frac{1}{|{\mathcal{G}}_h{|}^2}\sum _{i,j\in {\mathcal{G}}_h}E(<s_{MA}^i{\mathrm{v}}_A+s_{DS}^i{\mathrm{v}}_D+v_r^i|s_{MA}^j{\mathrm{v}}_A+s_{DS}^j{\mathrm{v}}_D+v_r^j>)$$

$$=\frac{1}{|{\mathcal{G}}_h{|}^2}\sum _{i,j\in {\mathcal{G}}_h}E(s_{MA}^i{s}_{MA}^j)\Vert {\mathrm{v}}_A{\Vert }^2+E(s_{DS}^i{s}_{DS}^j)\Vert {\mathrm{v}}_D{\Vert }^2+E(s_{MA}^i{s}_{DS}^j+s_{DS}^i{s}_{MA}^j)<{\mathrm{v}}_A|{\mathrm{v}}_D>+E(<v_r^i|v_r^j>)$$

$$E(\Vert {\hat{\mathrm{v}}}_A{\Vert }^2)$$

$$=\left(1+\frac{\mathrm{Var}(s_{MA}^h)}{|{\mathcal{G}}_h|}\right)\Vert {\mathrm{v}}_A{\Vert }^2+\frac{\mathrm{Var}(s_{DS}^h)}{|{\mathcal{G}}_h|}\Vert {\mathrm{v}}_D{\Vert }^2+\frac{E(\Vert v_r^h{\Vert }^2)}{|{\mathcal{G}}_h|}$$

$$E(\Vert {\hat{\mathrm{v}}}_D{\Vert }^2)$$

$$=E(\Vert {\hat{\mathrm{v}}}_A{\Vert }^2)+E\left(\Vert \frac{1}{|{\mathcal{G}}_{ad}|}\sum _{k\in {\mathcal{G}}_{ad}}{v}^k{\Vert }^2\right)-2E\left(<\frac{1}{|{\mathcal{G}}_{ad}|}\sum _{k\in {\mathcal{G}}_{ad}}{v}^k|{\hat{\mathrm{v}}}_A>\right)$$

$$=\left(1+\frac{\mathrm{Var}(s_{MA}^{ad})}{|{\mathcal{G}}_{ad}|}\right)\Vert {\mathrm{v}}_A{\Vert }^2+\left(1+\frac{\mathrm{Var}(s_{DS}^{ad})}{|{\mathcal{G}}_{ad}|}\right)\Vert {\mathrm{v}}_D{\Vert }^2+2<{\mathrm{v}}_A|{\mathrm{v}}_D>+\frac{E(\Vert v_r^{ad}{\Vert }^2)}{|{\mathcal{G}}_{ad}|}$$

$$+E(\Vert {\hat{\mathrm{v}}}_A{\Vert }^2)-2(\Vert {\mathrm{v}}_A{\Vert }^2+<{\mathrm{v}}_A|{\mathrm{v}}_D>)$$

$$=\left(\frac{\mathrm{Var}(s_{MA}^{ad})}{|{\mathcal{G}}_{ad}|}+\frac{\mathrm{Var}(s_{MA}^h)}{|{\mathcal{G}}_h|}\right)\Vert {\mathrm{v}}_A{\Vert }^2+\left(1+\frac{\mathrm{Var}(s_{DS}^{ad})}{|{\mathcal{G}}_{ad}|}+\frac{\mathrm{Var}(s_{DS}^h)}{|{\mathcal{G}}_h|}\right)\Vert {\mathrm{v}}_D{\Vert }^2+\frac{E(\Vert v_r^{ad}{\Vert }^2)}{|{\mathcal{G}}_{ad}|}+\frac{E(\Vert v_r^h{\Vert }^2)}{|{\mathcal{G}}_h|}$$

Full text

A model of brain morphological changes related to aging and Alzheimer’s disease from cross-sectional assessments

Raphaël Sivera111Corresponding author at: Epione Research Project, INRIA Sophia-Antipolis, 2004, route des Lucioles, 06902 Sophia-Antipolis, France, raphael.sivera@inria.fr

Université Côte d’Azur, Inria Sophia Antipolis, Epione Research Project, France.

Hervé Delingette

Université Côte d’Azur, Inria Sophia Antipolis, Epione Research Project, France.

Marco Lorenzi

Université Côte d’Azur, Inria Sophia Antipolis, Epione Research Project, France.

Xavier Pennec

Université Côte d’Azur, Inria Sophia Antipolis, Epione Research Project, France.

Nicholas Ayache

Université Côte d’Azur, Inria Sophia Antipolis, Epione Research Project, France.

Abstract

In this study we propose a deformation-based framework to jointly model the influence of aging and Alzheimer’s disease (AD) on the brain morphological evolution. Our approach combines a spatio-temporal description of both processes into a generative model. A reference morphology is deformed along specific trajectories to match subject specific morphologies. It is used to define two imaging progression markers: 1) a morphological age and 2) a disease score. These markers can be computed regionally in any brain region.

The approach is evaluated on brain structural magnetic resonance images (MRI) from the ADNI database. The model is first estimated on a control population using longitudinal data, then, for each testing subject, the markers are computed cross-sectionally for each acquisition. The longitudinal evolution of these markers is then studied in relation with the clinical diagnosis of the subjects and used to generate possible morphological evolutions.

In the model, the morphological changes associated with normal aging are mainly found around the ventricles, while the Alzheimer’s disease specific changes are located in the temporal lobe and the hippocampal area. The statistical analysis of these markers highlights differences between clinical conditions even though the inter-subject variability is quite high. The model is also generative since it can be used to simulate plausible morphological trajectories associated with the disease.

Our method quantifies two interpretable scalar imaging biomarkers assessing respectively the effects of aging and disease on brain morphology, at the individual and population level. These markers confirm the presence of an accelerated apparent aging component in Alzheimer’s patients but they also highlight specific morphological changes that can help discriminate clinical conditions even in prodromal stages. More generally, the joint modeling of normal and pathological evolutions shows promising results to describe age-related brain diseases over long time scales.

Keywords: aging, Alzheimer’s disease, brain morphology, deformations, spatio-temporal model, imaging biomarkers.

1 Introduction

Age-related diseases are a growing public health concern with the aging of the population. For this reason, a precise description of aging would be useful to predict and describe the evolution of these diseases. In complement to the chronological age, i.e. the time elapsed since birth, one would like to estimate a biological age that reflects the current physiological, functional or structural status of an organ relatively to the aging changes. However there is no unified way to describe aging in a clinical context since aging is a complex process which affects every part of the body with specific mechanisms and specific rate. As a consequence multiple theories of aging have been proposed [Medvedev, 1990], leading to the definition of surrogate age variables based on the quantification of biological changes.

1.1 Modeling brain morphological aging

In this paper we focus on the aging of the brain based on the study of its shape evolution. The brain is not exempt from aging and a decline of cognitive processing speed, working memory, inhibitory function, and long-term memory is generally observed. This decline has been associated with neural activity changes [Park and Reuter-Lorenz, 2009] and it was also shown to be directly correlated with structural changes such as brain atrophy, cortex thinning and decrease of white matter integrity [Rosen et al., 2003, Rodrigue and Raz, 2004].

The normal brain morphology has been studied in image-based studies from the development stage to the most advanced ages. Measurements of brain structures (volumes, cortical thickness, etc.) have been performed for wide age ranges and the statistical analysis of the evolution of these measurements helps in providing an initial understanding of the normal brain shape evolution across life span [Good et al., 2001, Long et al., 2012]. These descriptions have been used to estimate models characterizing brain aging in order to highlight differences across brain areas [Hutton et al., 2009, Sowell et al., 2003]. The inverse problem, i.e. how to associate an age to a brain image, was also addressed. Models have been designed to estimate the chronological age [Cole et al., 2018] from an image but they can also be used to characterize abnormal evolutions. For instance, a mean brain age gap estimate was highlighted for Alzheimer’s patients [Franke et al., 2010]. More generally, these surrogate brain age estimates have been associated with an increase of risk factors for several age-related disorders such as cardio-vascular diseases [DeCarli et al., 2005, van Velsen et al., 2013]. In a longitudinal setting, a brain age measurement could be used to compare the evolution of several clinical conditions (see Figure 1).

1.2 Disease progression modeling

Brain aging is often associated with the development of neurodegenerative pathologies. For example, it is estimated that one in three people over 85 have Alzheimer’s disease (AD), the most common form of dementia [Alzheimer’s Association and others, 2017]. This disease comes with its own specific apparent brain morphological changes [Ohnishi et al., 2001] and computer-aided diagnosis techniques using neuroimaging features have shown promising results to classify and to predict clinical evolutions [Davatzikos et al., 2009, Klöppel et al., 2012, Schmitter et al., 2015]. Longitudinal studies provide us with multiple acquisitions at different times for every subject but the disease affects patients over several decades, starting even before the first symptoms occur, and few studies follow a significant number of subjects over such long times. Progression models have been developed to describe the global evolution and to put in relation the individual trajectories that could only be observed a limited number of times. They have been used to model the progression of biomarkers [Fonteijn et al., 2012, Donohue et al., 2014] but also directly brain shape [Cury et al., 2016] or spatio-temporal patterns in brain images [Koval et al., 2017, Schiratti et al., 2017]. These models produce good representations of the disease progression and can combine a variety of available biomarkers for patient monitoring [Bilgel et al., 2016, Lorenzi et al., 2019].

1.3 Toward a joint model of brain aging and disease progression

The morphological aging of the brain and disease progression have generally been modeled independently. However we know that the structural features used in diagnosis (e.g cortical thickness or atrophy patterns) are also generally related to age. Indeed aging and neurodegenerative diseases involve intertwined processes with entangled consequences. Surrogate age measurements have been used to support the disease characterization [Franke and Gaser, 2012] or to put aside the aging part in order to focus on the disease specific changes. Lorenzi et al. [2015] proposed to model the normal aging evolution to separate the contribution of aging from a remainder that is not explicitly modeled. This remainder is then used to describe the pathological evolution. However, this method does not propose an intrinsic model of the disease progression making it difficult to describe and to characterize disease specific changes.

In this study, we propose a generative model of the brain morphological evolution that jointly takes into account the normal aging and the disease effects. Our model is based on the approach proposed by Lorenzi et al. [2015] and gives a deformation-based description of subject trajectories. It extends the original approach by explicitly modeling the disease specific brain morphological evolution. In addition to the apparent morphological age computed in the proposed approach, it allows us to get a disease score, thus providing two morphological imaging biomarkers accounting for the progression of the two main ongoing processes: normal aging and Alzheimer’s disease. These biomarkers are estimated cross-sectionnally from a single structural MRI but the model estimation exploits longitudinal data to match as accurately as possible the morphological evolution of the subjects.

We introduce in section 2 the generative model used to represent the brain morphology. Section 3 focuses on the estimation of the model parameters and on the inverse problem solved to compute the morphological age and the disease score of a subject. Experimental results are then presented in section 4 in order to evaluate our model and parameter estimation procedure. We illustrate how the model helps to describe the evolution of subjects at different disease stages using the ADNI database. Then we show how the two proposed markers can help to follow the evolution of elderly patients. Finally, limitations and perspectives are discussed in section 5.

2 Definition of the generative model

In the sequel, we quantify differences between morphologies by spatial deformations that can be estimated from magnetic resonance images (MRI) through non-linear image registration. A deformation represents either a morphological difference between the anatomies of two subjects or a longitudinal evolution of one subject-specific anatomy. Therefore deformation based frameworks are well suited to define a parametric model of the morphology. Our approach involves a population template morphology which is parametrized by two progression markers: the morphological age and the disease score.

In section 2.1, we expose the main ideas behind the definition of our reference parametric space. Then in section 2.2, we explain how the morphological evolutions are modeled in our framework. Finally in section 2.3, we show how we use the deformation framework to model individual morphologies relatively to this reference.

2.1 A space of reference morphologies

In deformation based morphometry, a single morphology is classically used to approximate and represent a population. For a set of images $\{I_{k}\}$, we define a common reference image $T_{0}$, called template. The difference of morphology between the subject $k$ and the template is modeled with a spatial deformation $\phi_{k}$ and an intensity noise $\epsilon_{k}$ is added in the subject space accounting for local intensity variability. Therefore the images are modeled as follows:

[TABLE]

In our approach, we want to take into account two major processes that affect the brain morphology over time: the normal aging and the disease evolution. To do so, we model the effects of these processes on the template using a deformation $\Phi$ parametrized by two variables that measure the progress of each process: the morphological age $\lambda_{MA}$ and the disease score $\lambda_{DS}$. The two variables $\lambda_{MA}$ and $\lambda_{DS}$ can be seen as time variables and are scaled to correspond to years of standard evolution. In this model, $T_{0}\circ\Phi(\lambda_{MA},\lambda_{DS})$ represents the template morphology after $\lambda_{MA}$ years of normal aging and $\lambda_{DS}$ years of normalized disease progression.

In the ideal case, the morphological age is equal to the chronological age. Therefore, $T_{0}\circ\Phi(t,0)$ represents the morphology of a $t$ years old healthy subject. Similarly $T_{0}\circ\Phi(t,t^{\prime})$ would be the typical morphology of a AD patient of age $t$ with a disease duration of $t^{\prime}$ years. If one is able to associate an age $\lambda_{0}$ to the image $T_{0}$, and assumes that this is the image of a healthy subject, then it is natural to enforce $\Phi(\lambda_{0},0)=\text{Id}$.

The parametric subspace of images generated like this will be used as a reference:

[TABLE]

An initial model of the subject images is then:

[TABLE]

where $\lambda_{MA}^{k}$ and $\lambda_{DS}^{k}$ are the subject morphological age and disease score, while $\phi_{k}$ encodes differences specific to the subject morphology. Of course multiple options are available to combine the longitudinal deformation $\Phi$ with the subject specific changes $\phi^{k}$. Here we write this operation as a right-composition but this choice will be discussed in section 2.2 and a similar but symmetrical operation will later be used.

Also, in the reference space, trajectories parametrized by time $t\mapsto T_{0}\circ\Phi({\lambda_{MA}(t)},{\lambda_{DS}(t)})$ give possible morphological evolutions where the morphological age and the disease score can be seen as reparametrization of time. In particular, it defines two archetypal trajectories (i.e. ideal models of evolution): the normal aging template trajectory $t\mapsto T_{0}\circ\Phi(t,0)$ and the disease specific template trajectory $t\mapsto T_{0}\circ\Phi(\lambda_{0},t)$.

Here, we assume that each progression of the two major processes (aging and disease) can be described with only one parameter. This implies that the evolution of healthy aging is similar for each subject, following the normal aging template trajectory, even if the speed of aging may vary from one subject to the other. Similarly, the disease progression is described using a single trajectory and we combine both template trajectories to model pathological evolutions.

2.2 Modeling the morphological evolution of the template

In a simplified approach, the template trajectories are assumed to be geodesics in an appropriate deformation space. Geodesics define continuous paths that can be easily parametrized and constrained to allow regularity in time. They can be used to interpolate between two anatomies or to approximate more complex trajectories [Christensen et al., 1994, Wang et al., 2007]. In this paper, we use the stationary velocity field (SVF) framework [Arsigny et al., 2006] for its ability to describe complex and realistic diffeomorphic (smooth and invertible) brain deformations in a straightforward manner [Lorenzi et al., 2013]. In this framework, the observed anatomical changes are encoded by diffeomorphisms which are parametrized with the flow of SVFs. Within this setting, the metric between deformations is not chosen a priori even if we need a regularization criterion for the registration. To compute a deformation $\phi$ we integrate trajectories along the vector field $v$ for a unit of time.

[TABLE]

This relationship is denoted as the group exponential map $\phi=\operatorname{Exp}(v)$.

By writing $\Phi(\lambda_{MA},\lambda_{DS})=\operatorname{Exp}(\lambda_{MA}\mathrm{v}_{A}+\lambda_{DS}\mathrm{v}_{D})$, we propose a linear model in the SVF space (i.e. the space of the parameter of the deformations) parametrized by two SVFs $\mathrm{v}_{A}$ and $\mathrm{v}_{D}$. In particular, the two template trajectories are then separately parametrized: $\mathrm{v}_{A}$ controls the normal aging template trajectory and $\mathrm{v}_{D}$ the disease specific template trajectory.

For each subject, the processes are meant to be intertwined and this can be modeled in different ways depending on the parametrization of the trajectories, for instance a right or a left composition. The proposed linear combination of the parameters provides us a middle ground. Indeed in the SVF setting, the relationship between composition and the linear combination of SVFs is given by the Baker-Campbell-Hausdorff formula [Bossa et al., 2007] and the linear combination of the SVFs is equivalent to alternate between right and left composition with infinitesimal steps.

To sum up, the longitudinal deformation $\Phi$ modeling the effects of aging and the disease on a reference morphology $T_{0}$ is parametrized by two SVFs: $\mathrm{v}_{A}$ and $\mathrm{v}_{D}$. This ideal model generates a surface $\mathcal{T}$ of possible images describing the evolution of the template morphology:

[TABLE]

2.3 Individual morphological variability and generative model

An individual image is modeled as follows:

[TABLE]

where the choice of the intensity noise $\epsilon_{k}$ is implicitly related to the registration similarity metric. To specify the constraint on $\phi_{k}$, we define a subject specific residual SVF $w^{k}_{r}$ ($r$ stands for residual) such that:

[TABLE]

In this formula, $\operatorname{Exp}(w^{k}_{r})$ is approximately equal to $\phi_{k}$ given the first order of the BCH equation between composition and linear combination of SVFs. Moreover, we wish to have the subject specific deformation to encode what cannot be described using the template trajectories. That is why we impose $w^{k}_{r}$ to be orthogonal to both $\mathrm{v}_{A}$ and $\mathrm{v}_{D}$.

As we can see in Figure 2, the model parametrized by $\mathrm{v}_{A}$ and $\mathrm{v}_{D}$ allows us to characterize the subject morphology with two scalar variables, the morphological age $\lambda_{MA}^{k}$ and the disease score $\lambda_{DS}^{k}$, and a SVF $w^{k}_{r}$ for the subject-specific part. The orthogonality constraint makes the description of the subject uniquely defined. We denote by $w^{k}$ the subject-to-template deformation SVF:

[TABLE]

3 Estimation of the model parameters

The model parameters are of two kinds: the population parameters ($T_{0}$, $\mathrm{v}_{A}$ and $\mathrm{v}_{D}$) and the subjects parameters ($\lambda_{MA}^{k}$, $\lambda_{DS}^{k}$ and $w^{k}_{r}$). To tackle this joint estimation problem in a computationally efficient way, several assumptions are made:

When available, the longitudinal evolution in the template space can be approximated by the transported deformation estimated in the subject space. Parallel transport algorithms are commonly used in the geometrical analysis of longitudinal data. The use of geodesic parallelograms is in general an efficient way to bring individual trajectories in a common reference space [Lorenzi and Pennec, 2013]. In practice, it allows us to work only with intra-subject deformation to estimate the model population parameters. It simplifies the optimization and is also more stable as the intra-subject variability is in general smaller than the inter-subject one. 2. 2.

We also assume, while estimating the population parameters, that the aging speed and the disease progression speed are constant for all the subjects in the training set. 3. 3.

Intra-subject deformations are relatively small and smooth. Consequently, the registration regularization has less impact on the estimated deformation. This consideration allows us to estimate these longitudinal evolutions independently of the population model.

These assumptions allow us to efficiently decompose the problem of the parameter estimation. First, subjects with longitudinal data are processed independently and the intra-subject evolutions are modeled in the subject space. Then the population parameters ($T_{0}$, $\mathrm{v}_{A}$ and $\mathrm{v}_{D}$) are estimated using only intra-subject longitudinal evolutions. Finally, the subjects’ parameters are estimated.

3.1 Estimation of the template trajectories $\mathrm{v}_{A}$ and $\mathrm{v}_{D}$ in a given template space

In this section we suppose that we know $T_{0}$ and that we can compute the subject-to-template deformation $w^{k}$ for a reference time point. We also consider that we have longitudinal data for every subject.

First we address the inverse problem of estimating the intra-subject evolution parameters with the framework proposed by Hadj-Hamou et al. [2016]. Images are preprocessed, rigidly aligned to the MNI-152 template and then longitudinally registered. Intra-subject deformations between follow-up images and the baseline image are computed using non-linear registration. The resulting intra-subject model in the subject’s space is estimated using ordinary least square regression in the tangent space of SVFs. It is equivalent to the assumption that the deformation noises are centered, uncorrelated and have equal variance in the space of SVFs.

Then for a given template $T_{0}$ and subject-to-template deformation $w^{k}$, the intra-subject model can be transported using parallel transport in the template space to get $v^{k}$. This deformation can be decomposed along the template trajectories giving a morphological aging rate (noted $s_{MA}^{k}$), a disease progression rate (noted $s_{DS}^{k}$) and an orthogonal component (noted $v_{r}^{k}$):

[TABLE]

These progression rates are strongly related to the evolution of the morphological age and of the disease score (in a ideal euclidean case $s_{MA}^{k}=d\lambda_{MA}^{k}/dt$ and $s_{DS}^{k}=d\lambda_{DS}^{k}/dt$) and they will be used to normalize the speed of evolution. The estimation is done on two groups of subjects: a group $\mathcal{G}_{h}$ composed by healthy subjects and a group $\mathcal{G}_{ad}$ composed by patients diagnosed with AD. We assume that each healthy subject of $\mathcal{G}_{h}$ is aging at normal speed $s_{MA}^{k}=1,\forall k\in\mathcal{G}_{h}$ and does not have any evolution toward the disease $s_{DS}^{k}=0,\forall k\in\mathcal{G}_{h}$. Similarly, each patient of $\mathcal{G}_{ad}$ has a normal morphological aging rate $s_{MA}^{k}=1,\forall k\in\mathcal{G}_{ad}$ and a constant unit disease progression rate $s_{DS}^{k}=1,\forall k\in\mathcal{G}_{ad}$. Finally the subject specific components are assumed to be centered, uncorrelated and to have a fixed variance. The maximum likelihood problem writes:

[TABLE]

The solution of the optimization problem is explicit:

[TABLE]

We should however note that $\|\hat{\mathrm{v}}_{A}\|$ (resp. $\|\hat{\mathrm{v}}_{D}\|$) is a biased estimator of $\|\mathrm{v}_{A}\|$ (resp. $\|\mathrm{v}_{D}\|$). We detail the bias estimation in the Appendix 6.2.

3.2 Estimation of the template morphology $T_{0}$

The population specific template morphology is computed using the algorithm proposed by Guimond et al. [2000] by alternating the registration of subject images to the template and the recomputation of the template intensities. However, in our approach, the subjects’ images do not need to be registered to $T_{0}$ directly but to their projection on the template space. To tackle this problem, we propose an iterative procedure where we register the image to its current projection on the reference space. Algorithm 1 details the procedure with simplified notations in the general case where $w$ parametrizes the deformation between the image $I$ and $T$. The reference linear subspace of SVF is denoted $\mathcal{T}$ and $w$ is decomposed accordingly $w=w_{t}+w_{r}$ with $w_{t}\in\mathcal{T},w_{r}\in\mathcal{T^{\perp}}$. The registration regularization should only be applied to the residual part $w_{r}$. In the context of the LCC-demons registration algorithm, it boils down to the following minimization problem [see Lorenzi et al., 2013]:

[TABLE]

The idea is to alternate between the optimization and the projection on the constraints.

Since we do not have any theoretical guarantee on the convergence of the algorithm, the stability and the convergence will be evaluated empirically. As the template estimation also involves iterative search, we can combine both algorithms for a faster optimization. The projection coordinates are kept from one iteration to the next and the images are registered to their estimated projections in the template space. The deformation update $u$ is then computed and finally the new atlas image $T$ and the estimated projections are updated (see Algorithm 2).

In this work, the intra-subject models are transported to the template space to update the template trajectories parametrizing $\mathcal{T}$ at each iteration (using the approach described in the previous section 3.1). $T$ is initialized using the MNI-152 template and the convergence is manually assessed comparing the template for successive iterations. At convergence, we obtain the template image $T_{0}$ and both template trajectories $\mathrm{v}_{A}$ and $\mathrm{v}_{D}$.

3.3 Estimation of the subject’s parameters

When the population parameters are learned, the estimation of the individual parameters for a new subject is relatively simple. The deformation $w^{k}$ is computed by registration between a subject image and the template using Algorithm 1, and then linearly decomposed, $w^{k}=(\lambda_{MA}^{k}-\lambda_{0})\mathrm{v}_{A}+\lambda_{DS}^{k}\mathrm{v}_{D}+w^{k}_{r}$, by solving the following linear system:

[TABLE]

In practice, the estimation is not exact because we work with the noisy estimators $\hat{\mathrm{v}}_{A}$ and $\hat{\mathrm{v}}_{D}$. The linear decomposition can also be computed locally or by using any voxel weighting scheme for the scalar product. When longitudinal data is available, this estimation is independently done for each time point.

4 Results

4.1 Experiments with synthetic data

We first evaluate our approach using synthetic data in order to assess the accuracy and the reproducibility of the biomarkers estimation. Realistic longitudinal MRIs are simulated using the software proposed by Khanal et al. [2016]. The simulation algorithm relies on a biophysical model of brain deformation and can be used to generate longitudinal evolutions with specific atrophy patterns. In this context, local atrophy is measured by the divergence of the stationary velocity field.

4.1.1 Simulated dataset

In this controlled experiments we choose to simulate two populations that are characterized by their atrophy patterns and that respectively emulate healthy controls and AD patients. Atrophy of the aging brain and the effect of AD have been extensively studied [Pini et al., 2016] and the atrophy measurements may vary depending on the methodology and the population studied. In this experiment, we choose to prescribe piecewise-constant atrophy map with constant value in brain areas delimited by the segmentation provided by FreeSurfer [Fischl et al., 2002]. For every subject, the atrophy value of a region is sampled around a fixed population mean with an additive Gaussian noise of relative standard deviation of 5% . The healthy population is designed to have a small atrophy in the whole brain while the patients have a stronger atrophy especially in the hippocampal areas and the temporal poles. The means are chosen to give the order of magnitude of a one year evolution accordingly to what was reported in Fjell et al. [2010] for healthy aging and in Carmichael et al. [2013], with an additional scaling for the pathological evolution. We detail the exact regional values in Table 2 in the appendix (see 6.1).

Structural MRIs of 40 healthy subjects from the ADNI database are taken as input to the simulations. For every subject, deformations are simulated for both pathological and healthy settings. The deformation extrapolated 5 times is then applied to the original image in order to simulate 5 years long evolutions. We then have two matched populations of 40 pairs of images.

4.1.2 Model estimation

Individual longitudinal deformations are computed using registration, and the reference anatomy and the template trajectories are built using our framework. The divergence fields associated with these template trajectories can be compared to the prescribed atrophy. Figure 3 shows the average atrophy maps in the estimated template anatomy.

The estimated atrophy patterns are smoothed versions of the simulated ground truth. This effect was already observed [Khanal et al., 2016]. First of all, the registration algorithm is unaware of the underlying simulation model and is unable to localize precisely the atrophy in homogeneous areas. Moreover, the spatial regularization of the registration and the parametrisation using SVFs also contribute to the smoothness of the estimated atrophy patterns. This is particularly visible in small (hippocampus) or thin regions (cortex). We therefore have a consistent bias when the atrophy measurements are integrated over the regions (see Figure 14 in appendix). Indeed, the local atrophy is affected by neighboring regions evolving in the opposite direction (the ventricles or the CSF for example). However, we can see that the ratio between the pathological and healthy cases is conserved in every region. It was already noted that quantitative estimation using registration can be biased but can be more reproducible than the segmentation based approaches [Hadj-Hamou, 2016, Cash et al., 2015].

4.1.3 Imaging biomarkers estimation

The morphological age and the disease score were computed for each image. By construction there is no difference at t=0 (exactly the same images in the two groups). We compare the simulated differences at t=5 and the evolution of the cross-sectional assessments for each subject (see Figure 4a).

In this experiment, the initial anatomical variability is important, indeed the standard deviation at t=0 is equal to 48 years for MA(0) and 24 years for DS(0). At t=5, a difference is visible for the disease score but it is still diluted in the inter-subject variability. It is also possible to extrapolate the evolution to determine how many years of evolution are needed in order to get a significant difference between the healthy and the disease groups. For a significance level of 0.05, the disease score would be significantly discriminant after 13 years of evolution. This figure for synthetic data highlights the slow time-pace of the disease and the interest in modeling and extrapolating the evolutions.

Looking at the evolution of these cross-sectional biomarkers in Figure 4b), we see that the measures are relatively stable despite the large inter-subject variability and that the difference measured between two different time points gives a good estimate of the longitudinal evolution. In practice the estimations are slightly biased, for example the increase of morphological age $\Delta MA=MA(5)-MA(0)$ is expected to be equal to 5 for both populations while the the mean of the estimation is equal to 3.92 for the healthy group and to 3.57 for the patients group. More importantly, the standard deviation is small in comparison to the standard deviation of the cross-sectional measurement ($\sigma_{\Delta MA}=1.13$ while $\sigma_{MA(0)}=47.8$). For the change of disease score $\Delta DS$, the mean is equal to $-0.6$ for the healthy group and $2.8$ for the patients group and it approximates the ideal expected values (respectively 0 and 5). The variance is also very small with respect to the cross-sectional one. In particular the difference between healthy and diseased subjects is clearly observed for the longitudinal evolution: the difference between the means is equal to 3.2 standard deviation for $\Delta DS$. To summarize, the cross-sectional measurement gives a relatively stable assessment of the aging and disease progression and the markers’ evolution is strongly associated with the clinical diagnosis. Our generative model is able to explain most of the independently simulated changes.

4.2 Experiments on ADNI data

Longitudinal T1 sequences were obtained from the Alzheimer’s Disease National Initiative (ADNI) database. The ADNI was launched in 2003 with the primary goal to test whether magnetic resonance imaging (MRI), positron emission tomography (PET), other biological markers, and clinical and neuropsychological assessment could be combined to measure the progression of mild cognitive impairment and early Alzheimer’s disease333see www.adni-info.org for more information.. Subjects are classified according to the evolution of their cognitive diagnosis. Three diagnoses are possible at each time point: normal, mild cognitive impairment (MCI) and Alzheimer’s disease. The subjects are also sub-classified according to the positivity of the beta-amyloid 1-42. We then have 6 distinct sub-groups: CN- (cognitively normal with negative $A\beta$), CN+ (cognitively normal with positive $A\beta$), MCIs- (MCI stable during the study time-window with negative $A\beta$), MCIs+ (MCI stable with positive $A\beta$), MCIc (MCI converter to AD) and AD (diagnosed with Alzheimer’s disease starting from the beginning). The table 1 sums up the demographic description of the population.

We estimate our template morphology and the template trajectories on a subset of subjects. In order to form this training set, we randomly selected 30 subjects from the CN- group and 30 from the AD group. To reduce the variability associated with the estimation of the model, these subjects were selected among the ones with strictly more than one followup acquisitions. In the following we distinguish between the training set of 30+30 subjects used to build the model and the remaining testing set (with in particular 78 CN- subjects and 173 AD subjects).

4.2.1 Estimation of the normal aging and the disease-specific template trajectories

The template anatomy is an average of the healthy subjects anatomies, so its age corresponds to the mean group age $\lambda_{0}=73.46$ y.o. The result of the estimation is shown in Figure 5. The estimated normal aging template trajectory is characterized mainly by ventricular expansion caused by the atrophy of the surrounding regions. Disease specific changes are widespread in the brain with a strong emphasis on the temporal areas.

The local atrophy can be measured by the divergence of the velocity field. The atrophy associated with each template trajectory is shown in Figure 6. For the normal aging, we see a well spread and mild atrophy pattern in the whole brain. The disease specific atrophy is particularly strong in the temporal area and is mainly located near the cortical surface. The precise localization of the atrophy is always difficult with a morphometrical approach but the atrophy patterns are similar to what was already observed in the past for healthy subjects and AD patients [Hadj-Hamou, 2016].

These morphological evolutions can be compared to the normal aging model and the mean residual for the AD subjects in Lorenzi et al. [2015]. The deformation characteristics and the magnitude of the atrophy are really similar. Our anatomy looks sharper which may be partially explained by the use of the template plane in the estimation. More importantly, differently from the work of Lorenzi et al. [2015], the joint estimation of a disease model in addition to the normal aging one provides us with a direct comparison between both processeses.

The generative model can also be used to directly visualize the modeled morphologies. In particular the reference template plane $\mathcal{T}$ described in Figure 2 can be sampled to shows the evolution of the template morphology in the two main directions. In Figure 7, we choose to represent the evolution over 20 years in both directions which is comparable to the longitudinal span of our data-set. Indeed the IQR of baseline age in the training set is equal to 9.3 years (and the total span is 29.9 years wide). The overlaid difference of intensity is used to show the changes at tissue boundaries. The global atrophy and the expansion of the ventricles is clearly visible for the aging evolution. The pathological changes are associated with smaller structures but the shrinking of the hippocampi, the atrophy of the temporal lobes and also the widening of the sulci (related to the cortical thickening) are visible.

4.2.2 Intra-subject variability of our progression markers

In Lorenzi et al. [2015], a morphological age similar to our measure was shown to be correlated with the chronological age and also that advanced AD stages were associated with “morphologically older” brains. To go further, we want to show that our proposed model represents also the aging at the individual level. For multiple acquisitions of a same subject, an aging measurement is expected to increase smoothly with time. If the subject is healthy, we can expect a linear increase with a slope of 1. We would also like to see an increase of disease score for the patients, while for the healthy subjects this marker should be stable and close to 0. Therefore, the disease score is expected to specifically characterize the pathology (high and increasing disease score). The morphological age is expected to be more associated with the global evolution of the population over time and to be independent of the clinical condition.

The morphological age and the disease score are computed for each subject at each time point. Figure 8 shows the evolution of these cross-sectional measurements. First, at the population level, subjects are generally associated with a morphological age similar to their chronological age even though its variability is quite high. Second, for each subject, the evolution is mostly linear and the morphological age steadily increases. Third, the disease score is almost constant for healthy subjects and steadily increases for the AD subjects. Finally, we note that the AD subjects look older, age faster and, most of all, have a higher disease score than the healthy ones.

A linear random effects model can help us to quantify these observations. The model is fitted, for both morphological age and disease score, with fixed effects on age and sex and a random intercept and slope for each subject. The focus is set on the analysis of the regression for the CN and the AD groups. For each coefficient of the regression, the confidence intervals are given for one standard deviation of the estimation.

The model is first fitted to the morphological age measures in the CN group leading to a coefficient of $0.26\pm 0.11$ for the fixed effect of age while the mean subject slope is $0.10\pm 0.02$. Both are significantly positive. In comparison to the same model without the random slope the relative improvement brought by the intra-subject linear evolution is significant by a large margin (p-value inferior to $10^{-6}$ for the likelihood ratio test). The regression has also a positive (but not significantly) coefficient for male subjects ($1.81\pm 1.2$) meaning that male morphologies looks older (similar to a $7$ years shift). Concerning the disease score, we also observe a relatively good fit of the linear model. The evolution is generally slower with $0.12\pm 0.1$ for the fixed effect of age and $0.12\pm 0.01$ for the mean individual slope.

For the AD group, the linear model is also well adapted (p-values inferior to $10^{-6}$ for the likelihood ratio test). The main remark is probably that the intra-subject slopes are in average more important than for healthy subjects (around $0.52\pm 0.06$ for the morphological age and $0.71\pm 0.05$ for the disease score) while the fixed effect related to age of $0.23\pm 0.07$ (for the morphological age) and $0.14\pm 0.04$ (for the disease score) are more similar to the one observed previously.

4.2.3 Cross-sectional discriminating power

We want to study the relation between the observed disease progression and the proposed markers. We start with a discriminant analysis using only the first image available for each subject.

Figure 9 shows the distribution of the estimations for each group. These results are related to global differences in brain shape observable cross-sectionally between clinical groups. We see a gradual increase of both markers towards more advanced disease states. Significant differences in morphological age and in disease score are observed between the control group CN-(train) and both the MCIc and AD groups. Moreover, the difference between the MCI stable and the MCI converters is stronger for the morphological age while the disease score better differentiates MCIs- and MCIs+. As such, the morphological age is more associated with the general cognitive degradation while the disease score seems more correlated with more AD-specific biomarkers.

We also perform a simple linear classification task between the MCIs and MCIc groups using this two cross-sectional markers. A SVM linear classifier is fitted to the full data-set to perform the binary classification444The classification task was performed using the SVC module of the scikit-learn python library [Pedregosa et al., 2011]. The error penalty weights are adjusted between the two classes to balance the trade-off between false positive and false negative rate. The mean classification accuracy using a 10-folds cross-validation scheme is equal to $0.59$. The linear decision function is equivalent to the projection on the SVF $\mathrm{v}_{A}-0.003\mathrm{v}_{D}$, so the differences between MCIs and MCIc subjects is, in our model, only associated with the aging trajectory while the disease specific changes do not seem to have an impact before the conversion. For comparison, for the same experiment between CN and AD, the linear classifier corresponds to the projection on $\mathrm{v}_{A}+0.49\mathrm{v}_{D}$ so approximately $(v_{bn}+v_{ad})/2$, i.e. the mean evolution of the whole population. Of course, in both cases we do not reach the performance of state-of-the-art dedicated algorithm but it allows us to see how both markers are associated with the diagnosis. Moreover this discriminant approach could be extended by using information in a subset of targeted areas.

4.2.4 Regional analysis of the progression

In the context of Alzheimer’s disease, specific morphological changes are known to be non-uniformly distributed. This spatial information can be taken into account in our model using a regional segmentation. In this section we focus on the temporal area which is often associated with Alzheimer’s disease [Double et al., 1996]. Using the AAL atlas [Tzourio-Mazoyer et al., 2002], we segment the temporal lobe of our template anatomy. The mask is then used to compute the regional morphological age and disease score for each subject. These markers only encodes the morphological differences in this specific area visible at the time of acquisition (i.e. at baseline in this case). Results are shown in Figure 10.

The region is more adapted to the disease score than to the morphological age model. Indeed, for a healthy subject the deformations in this area are really small. However the choice of a disease-adapted region is improving the performance of the disease score. It is now able to capture early specific changes and the difference between CN- and CN+ is significant.

4.2.5 Longitudinal evolution of the markers

To explore more in details the longitudinal evolution of these markers, a linear model is fitted to the individual evolutions. The intercept can be interesting as it aggregates the measure at every time point and helps reduce the noise but more importantly the slope can be very informative. Results are shown in Figure 11 for the whole brain markers.

A progressive evolution, from CN- to AD, is visible for the morphological age with subjects evolving faster and faster. Concerning the disease score, the evolution is almost negligible for CN- and MCIs- and relatively slow for CN+ and MCIs+ while the changes are clearly visible for MCIc and AD. Significant differences are visible between healthy subjects (CN-) and MCIc or AD subjects or even between MCI stable and MCI converters, but also between MCIs- and MCIs+ (or more generally between subjects with negative amyloid or positive amyloid marker). This may indicate that our measures are able to capture the global progression of the disease. The changes are larger for diagnosed patients but similar patterns of evolution are observed in the early stages of the disease. A significant difference is also observable between the CN- and CN+ group for the temporal disease score slope.

The disease score evolution is close to zero for the healthy subject and close to one for the pathological one but more generally, the slopes are in average smaller than their expected values. For example the average disease score evolution in the AD group is equal to $0.82$ and this discrepancy is particularly important for the morphological age slope of the CN- group that is only equal to $0.33$ (instead of $1$). This bias was already observed previously and can be in part explained by the estimation procedure (see  6.2).

4.3 Generating diagnosis driven morphological evolution

One of the main advantages of our model is its ability to be generative. From a pair of biomarkers (morphological age, disease score) we can generate a corresponding morphology and to deform a specific subject anatomy in the directions defined by the template trajectories. In this section, the model is used to generate plausible morphological evolutions of a subject, for several diagnosis condition, and compare them to the observed one.

4.3.1 Modeling the group evolution

In what follows, each subject is associated with a pair (position, speed) in the marker plane using the coefficient of the linear regression. For each group, a continuous vector field is regressed from these observations. This field estimates the local average speed at every point. In order to keep a local consistency in the markers space, we use a kernel ridge regression with an RBF kernel. The spatial scale is set to 10 years, for both the morphological age and the disease score axis, to get large scale patterns despite the high inter-subject variability. The regularization weight, which does not seem to have a large effect on the result, is set to 1. Results are presented in Figure 12 for the CN (i.e CN-, CN+ together), MCIs, MCIc and AD groups. The figure is centered on the high data density domain (as expected the extrapolation can be less reliable in lower density sectors).

Differences in amplitude, i.e. speed of evolution, and orientation are clearly visible between the groups and are in agreement with the linear regression results shown in the previous section. In particular we see a progressive amplitude increase from the CN group to the AD group.

These diagrams also help to describe the variability within the same clinical group. For the CN group we can distinguish between the low disease score and low morphological age area (in the bottom left) where in average the changes are negligible, and the rest where there is a slow horizontal evolution. This pattern suggests that the healthier and younger subjects are morphologically stable and do not show the same visible aging process. The MCI stable evolution is relatively uniform and in average with slightly larger amplitudes but overall similar to the CN one.

The MCI converters however show a stronger and more vertical evolution. We should also note that subjects with high morphological age and low disease score (bottom right) seem to follow a different, more horizontal evolution implying a fast morphological aging but less important disease specific changes. The AD group confirms this trend and in fact MCIc and AD look very similar. The mean evolution is strong and more vertical. A main evolution is visible from bottom left to top right with a sightly more horizontal part in the middle giving this global tangent-like aspect. Beside, a horizontal evolution, similar to what was observed for the MCIc model, is also visible in the bottom right. This difference of evolution suggests a possible stratification of the disease in two sub-categories.

4.3.2 Generating a subject specific evolution

The regression model can then be used to simulate the evolution of a given subject. We choose here to model the evolution of MCIc subject in order to predict the changes subsequently observed around the time of diagnosis. Several evolutions were computed for the subject 0361. This subject is chosen among the MCIc subjects because it has the longest time interval, here 8.5 years, between the first and the last acquisition. From the starting point, the markers at t=0, we integrate a trajectory using the speed given by the regression model at each point. Each clinical group is associated with its own model and then a different trajectory is computed for each diagnosis.

These markers changes can be directly translated in brain images to visualize the morphological evolution. Figure 13 shows the results for the end point of the trajectories for the CN and AD models. Images are generated by deforming the baseline image using the simulated deformation transported in the subject space. A bias correction is applied to the markers’ estimation using the value estimated on the synthetic dataset (the measured changes only correspond to 80% of the expected value).

For both subjects, the cortical atrophy, in particular in the temporal lobe, is clearly visible in the real image and to a lesser extent in the AD simulation. The expansion of the ventricles is always clearly visible, but even in the AD simulation, the volume change is inferior to the observation. Subject 0566 is, from the baseline, generally more atrophied and by comparison it makes some real changes less visible (in particular the sulci widening). This difference in baseline anatomy is interesting here for two reasons: first the deformations match correctly the anatomical structures (the hippocampi for example) despite the anatomical variability, second it generated two different evolutions for the same diagnosis. For the AD simulation, the deformation seems stronger in the temporal lobe and relatively weaker in the ventricles for the second case (0556). These evolutions are learned from population trends and even if they are not predictive for a particular subject, they are not aberrant in comparison with the real evolutions. Overall, the morphological changes simulated looks realistic even if they do not match perfectly the observed changes. Other aspects of the evolution are hard to quantify and often poorly documented. For example the evolution of the shape of the ventricles of subject 0361 is different in the three images. It may be related to different spatial distribution of the degeneration in the brain and inherent mechanical constraints. We also observe a global motion towards the bottom of the temporal lobe and a local rotation in the image for the real case and the AD simulation.

5 Discussion and perspectives

In this work, we proposed a novel deformation-based approach to measure the progression of normal and pathological processes from their effects on brain morphology. In the context of Alzheimer’s disease, it provides a simple description of the brain morphological evolution for elderly patients using only two degrees of freedom: an aging measurement and a disease score. The advantages come from three main properties.

First, we disentangle the aging and the disease progression using interpretable image-based biomarkers. Second, these markers are cross-sectional assessments and are consistent for intra-subject longitudinal analyses. They can be seen as alternative aging measurements compatible with ongoing biological processes. In particular the disease specific evolution appears to be associated with a positive amyloid marker even in prodromal stages. Third, we show that the markers and the generative model can be used in a personalized image simulation setting. It allows us to generate smooth and realistic evolutions for several diagnosis conditions.

Biological (or here morphological) age estimates were proven to be interesting to analyze the patient condition; here, in addition, the disease score is used to get a simple marker of the disease progression. The joint modeling gives a more complete description of the disease progression than a brain-age metric. The evolution is not seen as a simple accelerated aging process or the divergence from the normal evolution. On the contrary, it seems possible to capture the general worsening with the morphological age while the disease score measures an additional evolution that is more specific to the disease. Both patterns of evolution appeared to be related with the development of Alzheimer’s disease and this approach provides an intuitive interpretation and a simple decomposition of the morphological changes observed in mass-univariate morphometric approaches. Further analysis should be conducted to analyze the relation of these complementary patterns to clinical and cognitive variables.

Moreover the ideas behind the decomposition and the longitudinal/cross-sectional estimation of the parameters could be used in a different setting: the geodesic description is much more general than the SVF framework applied to structural MRIs. Other geodesic parametrizations could be more practical in another setting and even for the same image data, an LDDMM approach, for example, would generally concentrate the deformation more on the high intensity gradient interfaces while the SVFs model more spatially diffuse deformations.

Several approaches have been developed to model the longitudinal shape evolution using geodesic modeling [Singh et al., 2013]. Here we suggest that decoupling the estimation of individual trends using parallel transport could allow a more accurate and non-linear description of the disease while keeping a similar hierarchical mixed-effect structure.

As could be expected, a single variable is not enough to precisely describe the morphological changes that can occur while aging. Likewise, the inter-subject variability limits the accuracy of the modeling of the evolution of individual subjects. Some limitations come from the error and the approximation in the estimation of the model, others are related to aspects that are not taken into account and would require to modify the approach.

5.1 Approximations in the model estimation

One of the limitations of our approach is the difficulty to accurately estimate the markers. For example, some time points look like outliers when we perform the estimation for successive images of the same subject. This might be due to MR distortions that cannot be completely corrected by the registration. As a result, some artifacts are still visible in the estimated deformations. A better understanding of the effects of these distortions on the registration results would be useful to improve our additive deformation noise model. A first step could be to work with data where the ground-truth deformation is more controlled (scan-rescan images for example).

In this work, we used an orthogonal projection using the $L_{2}$ scalar product in the SVF space to define the subject specific deformation. This choice is arbitrary and using another region of interest could change the result. A possible alternative would be to match the subject with the closest morphology in the model for a metric more adapted to deformations. We could also try to decorrelate the two information, using ICA for example.

We saw, using synthetic and longitudinal data, that we were under-estimating the evolution speed. This bias is partially explained by the estimation procedure that does not take into account the uncertainty on the population parameters of the model (see 6.2). This problem cannot be easily solved for two reasons: first the bias caused by the estimation of the template anatomy using the healthy controls groups is hard to quantify, second an unbiased estimation would be far more complex and would not be possible without the full knowledge of the training set for every new subject. Another source of error can be the registration algorithm: the inter-subject registrations are larger than the intra-subject ones and no transitivity is guaranteed.

5.2 Limitations of the description

Probably the main limitation of our approach is the inter-subject variability of the markers. The markers are very stable in longitudinal but the anatomical variability is not completely absorbed by the biomarker estimation and they are sensitive to the subject differences. The use of an explicit or even implicit model of the shape of a normal brain would help to makes them more specific to time related changes.

We could also note that the template anatomy is estimated using healthy subjects only. This modeling choice may introduce a bias towards the healthy population. We argue that the explicit modeling of the disease in our model could help mitigating this issue. Moreover, this choice makes the origin of the morphological age and the disease score simpler to define. Finally, pulling all the subjects could also enable a more accurate estimation [Marchewka et al., 2014] but the relatively large and homogeneous populations involved does not make it critical here.

In this context, the use of a single reference anatomy to parametrize the template space could also be discussed. Here, for example, it introduced a bias toward a certain age because of the way we composed deformations. A multi-atlas approach could be a better solution if a single anatomy is not enough to describe our population [Blezek and Miller, 2007, Sabuncu et al., 2009]. And we could do something similar to what was proposed by Rohé et al. [2016] to intricate the registration and the template subspace prior.

Regarding the estimation of the template trajectories using individual longitudinal SVFs transported in the template, potential bias associated with the asymmetrical use of the baseline image (in processing and in modeling) has been highlighted [Reuter et al., 2012] and the imbalance in number of time points or total longitudinal span should be taken into account for a more reliable estimation. However, Hadj-Hamou et al. [2016] has shown that the asymmetry caused by the non centrality of the time point in the longitudinal sequence is not completely relevant in the stationary velocity field framework because the SVF are expressed in Eulerian coordinate and should be identical at all time points.

5.3 Perspectives

The use of segmentation to compute the progression markers in multiple regions would be another way to extend the description. We showed that using a segmentation of the temporal lobes could tighten the link between the morphological markers and early clinical conditions. However the question of the regional interactions is not addressed in this work and the spatial analysis of brain deformation remains a research topic.

Our model of the morphological evolution is generic and of low dimensionality. Consequently it only partially captures the changes related to AD and, in practice, the subject specific field $w^{k}_{r}$ is actually encoding changes that are related to the disease and its variability but are not currently modeled. To enhance the description, the model should integrate the observed and latent clinical condition of the subjects. A description that goes beyond the global cognitive diagnosis would be interesting to pursue in the temporal description of the morphological changes in order to better describe the disease progression and capture this evolution from a healthy state to a pathological one.

Beside, the observed longitudinal evolution is also influenced by various factors such as sex, genetic, other pathologies or even the image acquisition protocol. The ability to handles these covariates has shown its utility to better describe the brain morphological evolution [Muralidharan et al., 2016].

These extensions would also improve the generative aspect. Coupling our approach with a proper disease progression model, and using a mixture model for the trajectories, would enable the generation of morphological trajectories in a more diverse setting to explore and sample the range of possible evolutions.

Acknowledgment

This work has been supported by the Université Côte d’Azur 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) (National Institutes of Health Grant U01 AG024904) and DOD ADNI (Department of Defense award number W81XWH-12-2-0012). 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.

Authors have no conflict of interest regarding this article

6 Appendix

6.1 Validation on a synthetic dataset

The regional values of atrophy set are given in Table 2. For each subject, the values are sampled around these means with a 5% standard deviation.

The comparison with the values obtained after simulation or estimated through registration are shown in Figure 14. The estimation can be biased by the spatial regularization and the loss of information in intensity homogeneous areas. The relative changes is however similar between the two populations and the method is by consequence adapted to compare the two evolutions.

6.2 Bias on the estimated template trajectories

The norms of the SVFs parametrizing the template trajectories have an effect on the normalization of the individual biomarkers estimation. We estimate here the bias on the the norm of the estimation relatively to the estimated norm.

For the normal aging trajectory we have:

[TABLE]

Assuming that the subjects are independent and identically distributed:

[TABLE]

And similarly for $\|\hat{\mathrm{v}}_{D}\|$:

[TABLE]

The terms $\operatorname{Var}(s_{MA})$, $\operatorname{Var}(s_{DS})$ and $\|v_{r}\|^{2}$ are related to the individual variability that is not modeled and to the noise in the estimated deformations. We empirically estimate the variances in the training population and we get that, for our training dataset, $\|\hat{\mathrm{v}}_{A}\|^{2}\approx 1.65\|\mathrm{v}_{A}\|^{2}$ and $\|\hat{\mathrm{v}}_{D}\|^{2}\approx 1.21\|\mathrm{v}_{D}\|^{2}$. We should however note that the same subjects from $\mathcal{G}_{h}$ are used to estimate the template anatomy $T_{0}$ and the normal aging template trajectory $\mathrm{v}_{A}$. Therefore, we are certainly underestimating the bias coming from the intra-subject morphological variability for this population.

This bias has a direct influence on the markers estimations. Indeed for a subject $k$, the markers are solution of the linear system 9 involving $\|\mathrm{v}_{A}\|^{2}$, $\|\mathrm{v}_{D}\|^{2}$ for which the biased estimators are used.