Keywords: Neurodegeneration, la maladie d'Alzheimer, la modélisation biophysique et la simulation de Neurodegeneration, Alzheimer's disease, biophysical modelling and simulation of atrophy, longitudinal MRIs simulation, biomechanical simulation, synthetic longitudinal images Atrophy Map . . . . . . . . . . . . Atrophy Maps . . . . . 2.4.2 Simulating Realistic Atrophy Patterns . . . . . . . . . . . . . 2.4.3 Simulating Large Atrophy With Multiple Time-steps . . . . . 2.4 Parameters . . . . . . . . . . . . . . . . . . . . . B.2.2 Minimization . . . .

Dans cette thèse, nous avons développé des outils pour simuler des images longitudinales réalistes de cerveau présentant de l'atrophie ou de la croissance. Cette méthode a été spécifiquement élaborée pour simuler les effets de la maladie d'Alzheimer sur le cerveau. Elle se fonde sur un modéle de déformation du cerveau qui décrit les effets biomécaniques d'une perte de tissue due à une carte d'atrophie prescrite. Nous avons élaboré une méthode pour interpoler et extrapoler les images longitudinales d'un patient en simulant des images avec une carte d'atrophie spécifique au sujet. Cette méthode a été utilisée pour interpoler des acquisitions temporelles d'Images par Résonnance Magnétique (IRM) de 46 patients souffrant de la maladie d'Alzheimer. Pour ce faire, des cartes d'atrophie sont estimées pour chaque patient, d'après deux acquisitions IRM temporelles distinctes. Les IRM cliniques présentent du bruit et des artefacts. De plus, les acquisitions longitudinales présentent des variations d'intensité d'une image à l'autre. Nous avons donc élaboré une méthode qui combine le modèle de déformation du cerveau, ainsi que les différentes images cliniques disponibles d'un patient donné, afin de simuler les variations d'intensité des acquisitions longitudinale. Pour finir, les outils de simulation d'images réalistes développés au cours de cette thèse sont mis à disposition en open-source.

1.1 Alzheimer's Disease (AD)

In 1907, Alois Alzheimer reported a study of one of his patients with "an unusual illness of the cerebral cortex" [START_REF] Alzheimer | über eine eigenartige erkankung der hirnrinde[END_REF], Stelzmann 1995]. Alzheimer detailed the symptoms of his patient suggesting that they were so unusual that the patient could not be classified as having one of the recognized illnesses. Indeed, Alzheimer was right when he said that he was dealing with a special illness, which is now named after himself, as Alzheimer's Disease (AD). Ever since this first reporting, several studies have been performed in the past century and in particular the last couple of decades. However, the exact mechanisms of AD and its causes are poorly understood, and there are still no cure to date. Progressive decline of cognitive functions of the human brain associated with the gradual death of neurons is a more general condition named dementia. Prominent symptoms of dementia are loss of memory, mood changes, and problems with communication and reasoning. AD is the most common form of dementia contributing to 60 -70% of cases [WHO 2015]. It primarily affects the aging population and the risk increases with age. The increment of longevity in most parts of the world due to the advancement of medical science and consequently better health care can cause a significant rise in the incidence of AD in the coming years. Thus, it is imperative to improve our understanding of AD to alleviate the growing socio-economic impact of AD in the world. [Bird 2008], which shows extracellular senile plaque (also called neuritic plaque) in the lower left corner and NFTs in the upper right corner. They are the most important characteristics of AD along with widespread brain atrophy and progressive cognitive impairment.

Characteristics of AD

The primary characteristics of AD before external cognitive symptoms appears are [Braak 1995]:

• the presence of neuronal atrophy (Figure 1.1),

Characteristics of AD 3

• abnormal accumulation of Amyloid-beta (Aβ) as senile plaques in extracellular matrix (Figure 1.2),

• abnormal accumulation of hyperphosphorylated tau protein as neurofibrillary tangles (NFTs) inside neurons (Figure 1.2).

Although the exact cause of AD is not known, there has been some hypotheses to explain the cause of AD, and they mostly focus on either the Aβ plaques or taupathy leading to NFTs. One such hypothesis is the amyloid cascade hypothesis where the accumulation of amyloid beta peptide (Aβ) as amyloid plaques leads to AD [Hardy 1992, Hardy 2002]. Many researchers support this hypothesis but some others consider tangle formation due to the aggregation of tau protein to be the major cause [START_REF] Di | Alzheimer's disease: biological aspects, therapeutic perspectives and diagnostic tools[END_REF]]. The possible causes of AD can be classified into three groups of deficiencies shown in Figure 1.3: cellular, molecular and genetic imbalances [START_REF] Herrup | The case for rejecting the amyloid cascade hypothesis[END_REF]]. showing the risk factor genes where PSEN2, PSEN1 and APP are emphasized to indicate their involvement in familial AD (fAD). See [START_REF] Herrup | The case for rejecting the amyloid cascade hypothesis[END_REF]] for detailed description.

The brain shape changes due to atrophy can be observed in longitudinal timeseries structural magnetic resonance images (MRIs). This thesis focuses on developing plausible biophysical model of brain shape changes as a consequence of neuronal deaths and to use the model to simulate longitudinal images that can reproduce the changes we observe in real images.

Before detailing the context and motivation of the thesis, we briefly present some of the established biomarkers that are based on the three primary characteristics of AD mentioned above. These are the most widely studied aspects of AD. However, it is worth noting that there are also recent studies proposing other hypotheses and possible mechanisms [START_REF] Pisa | Different Brain Regions are Infected with Fungi in Alzheimer's Disease[END_REF], Lopategui Cabezas 2014[START_REF] Hong | Complement and microglia mediate early synapse loss in Alzheimer mouse models[END_REF], Chung 2015]. As these studies are still at a nascent stage, and since it is not the objective of the thesis to completely understand the mechanisms of AD, they are not discussed in this chapter.

Biomarkers for AD

The entities whose presence, activity or concentration indicate the presence or severity of a disease are known as the biomarkers for that disease. Several in-vivo markers of Alzheimer's pathology are available which aid in the diagnosis of AD [Dubois 2010]. These biomarkers have different relevance at different stages of AD progression. Jack et al [Jack 2010, Jack 2013] describe the temporal evolution of AD biomarkers and their interrelationships amongst themselves, and with the onset and progression of clinical symptoms. The staging of the disease follow a progression from normal to mild cognitive impairment (MCI) to eventually dementia. Figure 1. 4 shows the staging of the disease and the dynamics of biomarkers for AD as proposed by Jack et al in [Jack 2010, Jack 2013

Senile plaques

One major school of thought in AD links the cause of AD to Aβ pathophysiology supporting the amyloid cascade hypothesis proposed by Hardy et al. [Hardy 1992].

Figure 1.4: A hypothetical model of dynamic biomarkers of AD, reproduced from [Jack 2013] with permission. This model incorporates tau and Aβ pathology as independent processes trying to reconcile the following two seemingly contradictory existing hypotheses: i) Several AD biomarkers study suggesting Aβ pathophysiology occurs first followed by tau-related neurodegeneration. ii) Autopsy data suggesting AD like tauopathy precedes Aβ deposition [START_REF] Braak | The pathological process underlying Alzheimer's disease in individuals under thirty[END_REF]]. The model thus suggests that in many individuals the first AD pathophysiological process to arise is the subcortical taupathy which is detectable only by immunostaining methods. However, the tauopathy alone does not lead to AD and Aβ pathophysiology arising later independently from pre-existing tauopathy accelerate the spread of NFTs through unknown mechanisms.

The main idea of this hypothesis is that the deposition of Aβ, a peptide fragment of a membrane protein called amyloid precursor protein (APP), is the root cause of AD which directly results in the NFTs, neuronal deaths, vascular damage, and dementia. Aβ peptides aggregate to form oligomers1 that subsequently produce depositions of extracellular macroscopic Aβ plaques (also known as senile plaques).

There are only three genes known to lead to the early onset in familial form of AD (fAD): APP, PSEN1 and PSEN2 [Bird 2008]. Due to the role of these three genes in the sequence of senile plaques formation from the Aβ peptides, the amyloid cascade hypothesis has received stronger support [START_REF] Herrup | The case for rejecting the amyloid cascade hypothesis[END_REF]]. Although lots of studies have been based on this hypothesis, it is important to note that there are also studies disagreeing with it [START_REF] Herrup | The case for rejecting the amyloid cascade hypothesis[END_REF]]. The existing debate is in whether Aβ is the primary cause of AD or not. However, with the growing body of evidence from several studies, it is universally accepted that Aβ is strongly correlated to AD and it will remain as an important component in AD research.

Both CSF and imaging biomarkers targeting Aβ pathophysiology are available and have been widely used in AD research. The established CSF biomarker is the decreased level of amyloid-β 1-42 (Aβ42) in AD [Anoop 2010], while the most used imaging biomarker is the increased level of Pittsburgh Compound-B (PiB) seen in Positron Emission Tomography (PET) scans of AD patients [Johnson 2012]. Aβ 42 is 42-amino-acid soluble form of Aβ and is believed to be lowered when it gets converted to the insoluble Aβ plaques. This view is consistent to the findings of inverse relationship between CSF Aβ42 levels, and PIB levels which binds with Aβ plaques [Fagan 2006]. Figure 1. 4 shows the dynamics of these biomarkers where we see that biomarkers targeting Aβ pathophysiology are sensitive well before the cognitive impairment begins. Aβ related biomarkers are already in saturation phase for demented patients. Figure 1.5 shows the distribution of amyloid image levels in normal, MCI and AD patients.

Tau pathology

Microtubules are small tubular structures present inside human cells playing important role in maintaining the structure of cells and also in several other cellular processes. These microtubules are critical structures for stable neuronal morphology serving as tracks for long-distance transport, providing dynamic and mechanical functions, and controlling local signaling events [Kapitein 2015]. The major protein associated to microtubule in healthy neurons is tau which promotes the as-sembly of tubulin into microtubules and aids in the stabilization of their structure [Weingarten 1975].

The abnormal hyperphosphorylation of tau proteins and its aggregation into bundles of filaments can lead to a class of neurodegenerative diseases called as taupathies. AD is the best known tauopathy with the aggregation of hyperphosphorylated tau proteins leading to the deposition of neurofibrillary tangles (NFTs) within neurons. The presence of NFTs is one of the two major hallmarks of AD, the other being the senile plaques described in the previous section. The progression of cognitive impairment in AD follows more closely the spatio-temporal spread of NFTs rather than the senile plaques [Braak 1995, Jack 2013].

CSF biomarkers targeting tau proteins are available just like with Aβ, but the imaging biomarkers targeting tau are very recent developments and finding the best PET tracers that can bind with high sensitivity and selectivity to phosphorylated tau is an ongoing research topic [Villemagne 2015]. High levels of total tau (t-Tau) and of phosphorylated tau (p-Tau) in CSF along with low levels of Aβ42 have been shown to be highly predictive of MCI to AD converters [Shaw 2009]. Since the tau is an intracellular protein binding to microtubules, in the CSF of healthy subjects they are found in low numbers. In taupathies these normal tau proteins become phosphorylated and get dissociated from microtubules, eventually contributing to the formation of NFTs [Goedert 1993]. Tau protein could be released into CSF during the process of formation of NFTs and the subsequent disruption of neuronal architecture and cell deaths [Formichi 2006]. Hence the elevated levels of t-Tau and p-Tau in CSF can correlate with the the onset of neurodegeneration in AD [Anoop 2010]. There are several tau imaging ligands currently under research [Villemagne 2015] with fluroine-18 ( 18 F) isotopes based PET tracers being the most promising [Murray 2014].

Neurodegeneration

Main imaging biomarkers targeting neurodegeneration that are well established are FDG-PET scan and structural MRI; functional MRI (fMRI) and diffusion weighted imaging (DWI) have also been increasingly explored in their potential use as an early biomarker for AD [START_REF] Sperling | The potential of functional MRI as a biomarker in early Alzheimer's disease[END_REF][START_REF] Acosta-Cabronero | [END_REF][START_REF] Colliot | Magnetic resonance imaging for diagnosis of early Alzheimer's disease[END_REF].

fMRI is a noninvasive imaging technique which can be used to probe functional integrity of brain networks. It provides an indirect measure of neuronal activity obtained from measuring changes in blood oxygen level-dependent (BOLD) MR signal. It has been increasingly used to analyse the integrity of brain networks supporting memory and other cognitive domains in aging and early AD [START_REF] Sperling | The potential of functional MRI as a biomarker in early Alzheimer's disease[END_REF], Johnson 2012]. Task fMRI activation studies are usually not feasible in demented patients; for these patients only resting state fMRI may be feasible, although the head motion of these patients during scan still make it challenging to obtain better data. Furthermore, high variability and limited reproducibility studies of BOLD fMRI signals across subjects mean more studies are required exploring the utility of fMRI in studying AD progression. There is also a need of studies that Top two rows show that Aβ deposition could be either present or absent in NCs. In second row tau deposits in NC is seen only in hippocampal region. MCI and AD patients have progressively increased depositions of both tau and Aβ. Although Aβ seems to be deposited almost everywhere in the brain regions, tau deposits are more concentrated in the cortical regions. MMSE=Mini Mental State Examination.

track the evolution of changes in the fMRI activation patterns during the progression continuum from preclinical to clinical AD. Finally, the contribution of atrophy observed in sMRI to changes observed with fMRI must also be evaluated.

DWI generates contrasts on images based on the measurements of the random Brownian motion of water molecules within a voxel of tissue. Diffusion tensor imaging (DTI) and white matter tractography based on DWIs have been used to study white matter integrity in the human brain. Most of the studies on AD using DTI focus on measuring differences in white matter integrity of AD patients against normal controls. Although there has been several such studies, it is still unclear whether it can be considered as a reliable biomarker for AD [START_REF] Acosta-Cabronero | [END_REF]. Cerebral metabolic rate of glucose can be measured by PET scan using the tracer fluorine-18 (F-18) fluorodeoxyglucose (FDG). In AD, the reductions in cerebral metabolism is well established and is one of the two biomarkers targeting neurodegeneration. The glucose metabolism measures with FDG-PET scan have been used with good accuracy to distinguish AD from both NCs and other dementias [Mosconi 2007[START_REF] Ballard | [END_REF]. FDG-PET has also been used to predict MCI to AD converters with better accuracy than structural MRI [Yuan 2008].

Longitudinal structural MRIs have been widely studied as an imaging biomarker for AD. At present, brain atrophy is measured from the high-resolution acquisitions with MRI scanners of 1.5T or 3T magnets. The best established and validated atrophy assessment methods are based on T1-weighted MRIs [Frisoni 2010]. Progressive death of neurons or neurodegeneration leads to structural changes in the brain which can be observed in strucutral MRI. Structural changes seen in MRI correlates well with the cognitive impairment [Jack 2013]. By this time, Aβ abnormality is already saturated. Thus, for monitoring the impact of disease modifying drugs, the ability to track and predict structural changes in MRI can play an important role.

Figure 1.7 shows an example of the changes seen in the brain structure of an AD patient from the patient's two brain scans acquired two years apart. In the figure, we can see hippocampal volume loss, medial temporal lobe atrophy and ventricular expansion.

Context and motivation 10

Context and motivation

There has been several clinical trials and disease modifying drug development efforts in the past three decades [Schneider 2014]. Since the external symptoms appear several years after the changes seen in imaging, longitudinal structural MRIs can play an important role in the development of disease modifying drugs. So far, structural MRIs have primarily been used for estimating local volume changes in individual AD patients; these measurements have been used to test or formulate hypotheses on the temporal dynamics of atrophy in AD [Whitwell 2005, Sabuncu 2011[START_REF] Dubois | Revising the definition of Alzheimer's disease: a new lexicon[END_REF]].

An interesting alternative avenue consists in modeling the tissue loss process in order to compare (in a forward modeling setting) different hypotheses for the prediction of patient-specific time series MRIs. The development of reliable models of brain shape changes to predict accurately patient-specific longitudinal MRIs can have far reaching consequences. For instance, the patient-specific evolution of the brain structures predicted by the model could be useful in monitoring drug effects in AD patients by comparing the prediction against the observed brain changes.

For these models to be practically useful, they must be able to reliably predict the evolution of brain structures of the patients from available data such as spatial distribution of Aβ and tau proteins coming from multi-modal imaging sources, genetic information etc. These models could be validated by comparing the predicted MRIs against the observed real MRIs. Since the disease mechanisms are not well known, it is very challenging to develop such a comprehensive model that can predict time series of MRIs following AD patient's real evolution of brain shape. As far as we know, there has not been any studies which attempt to develop realistic models in this direction.

The effect of atrophy of the brain due to the death of neurons is well characterized using longitudinal MRIs. Recent developments in computational anatomy using longitudinal image registration allow computing statistics on the longitudinal brain changes in population of patients. For instance, Lorenzi et al. showed that deformation trajectories can be reliably computed from time series of structural MRIs in AD as well as in normal subjects [Lorenzi 2011]. However, the mean deformation trajectory is parameterized by a very large number of parameters which are difficult to understand.

One interesting direction to explore in this context is to develop plausible biomechanical models of brain atrophy that can explain the results observed from the images. Biomechanics based modeling have also been used to model the observed shape evolution during the development of the fetal brain. For instance, the role of the mechanical forces on the determination of the brain shape is well-illustrated by Tallinen et al.; using a biophysical model, the authors reproduced remarkably well the formation of the gyral and sulcal convolutions in a 3D-printed layered gel mimic of the developing smooth fetal brain [Tallinen 2016]. There is evidence that endogenous mechanical forces at the cellular level influence brain structure and function [Tyler 2012]. Although the detailed mechanisms of these interaction still deserve further investigation [Tyler 2012, Mueller 2015], it is clear that they play a role at the macroscopic level which is the scale where we observe changes in the structural MRIs. Thus, we believe that when developing macroscopic models of brain shape changes due to neurodegeneration, it is important to emphasise the link between model parameters and the assumptions about the real mechanobiology of neurodegeneration.

A plausible biomechanical model of brain atrophy could be used to simulate many different scenarios of evolution along the course of the disease. The resulting simulated longitudinal image sequences could then be used to evaluate the accuracy of the spatial location and timing of the changes along the course of the disease. Such a model might also be helpful in studies attempting to disentangle potentially different mechanisms of structural changes in AD and their relationship in time.

Biomechanical tissue deformation models have been used for the validation of non-rigid medical image registration. One of the first applications was on the breast MRIs by Schnabel et al., who proposed a validation framework using the images simulated with physically plausible, biomechanical tissue deformations [Schnabel 2003].

In the context of brain morphometry tools targeting neurodegeneration, a model that can generate realistic synthetic MRIs with prescribed local volume changes can be useful in the development, validation, and evaluation of: i) atrophy estimation algorithms ii) data-driven disease progression models. In addition to simulating the appearance of morphological changes in the brain, these synthetic time-series MRIs should also produce realistic variation present in real longitudinal images such as independence of noise and MRI acquisition artifacts. There are also several machine learning algorithms that attempt to automatically classify AD patients in large databases [START_REF] Cuingnet | Automatic classification of patients with Alzheimer's disease from structural MRI: A comparison of ten methods using the ADNI database[END_REF]]. Simulating a large database of realistic ground truth images could be helpful in the training of such classification algorithms too.

Structure of the Manuscript and Contributions

In this thesis, we develop a biophysical model of brain deformation with prescribed atrophy, which is based on plausible assumptions of brain tissue loss biomechanics in AD. Using the model of brain deformation as a core component, we also develop a framework to simulate realistic longitudinal brain images with atrophy/growth and intensity changes.

The manuscript is organised based on three main ideas as follows:

• Chapter 2 presents the first contribution of the thesis: development of a carefully designed biomechanics-based tissue loss model to simulate deformations with prescribed atrophy. For a given baseline brain MRI, the model yields a deformation field imposing the desired atrophy at each voxel of the brain parenchyma, while allowing the CSF to expand as required to globally compensate for the locally prescribed volume loss. The model allows simulating synthetic MRIs by prescribing complex patterns of atrophy. The chapter also presents a pipeline that allows evaluating atrophy estimation algorithms by simulating longitudinal MRIs with complex subject-specific atrophy patterns.

This chapter was published as a journal article in NeuroImage [Khanal 2016b], which builds upon a conference paper published in MICCAI [START_REF] Khanal | A Biophysical Model of Shape Changes due to Atrophy in the Brain with Alzheimer's Disease[END_REF]].

• Chapter 3 presents the second contribution of the thesis: A simple approach to interpolate or extrapolate new subject specific time-point MRIs by using the biophysical brain deformation model developed in Chapter 2. Atrophy rates per year in a large number of segmented brain regions are estimated from the available two extremal time-point scans of each patient. Middle time-point atrophy map is then predicted by assuming linear progression; the predicted atrophy map is prescribed to the biophysical model of brain deformation and a middle time-point image is simulated. The approach is evaluated on the MRIs of 46 AD patients. This chapter is adapted from our published work in [Khanal 2016c].

• Chapter 4 presents our third contribution where we propose a method to add realistic variations of intensities and noise to synthetic longitudinal image sequences as observed in the real ones. The proposed method uses the brain deformation model and different available images of a subject to add realistic variations of intensities in the synthetic longitudinal images. Simulated examples of realistic longitudinal MRI sequences are shown which contain the variations present in real longitudinal images, such as independence of noise and MRI acquisition artifacts. The software developed during the thesis to simulate realistic longitudinal brain images with our brain deformation model is named Simul@atrophy and will be soon released in an open-source repository. The chapter also explains various options available in Simul@atrophy.

The open-source release of the software is the final contribution of the thesis. This chapter is under review in NeuroImage [Khanal 2016a].

Finally, Chaper 5 concludes the thesis by summarizing the contributions and providing some perspectives. Abstract We propose a framework for developing a comprehensive biophysical model that could predict and simulate realistic longitudinal MRIs of patients with Alzheimer's Disease (AD). The framework includes three major building blocks: i) Atrophy generation ii) Brain deformation iii) Realistic MRI generation. Within this framework, this paper focuses on a detailed implementation of the brain deformation block with a carefully designed biomechanics-based tissue loss model. For a given baseline brain MRI, the model yields a deformation field imposing the desired atrophy at each voxel of the brain parenchyma while allowing the CSF to expand as required to globally compensate for the locally prescribed volume loss. Our approach is inspired by biomechanical principles and involves a system of equations similar to Stokes equations in fluid mechanics but with the presence of a non-zero mass source term. We use this model to simulate longitudinal MRIs by prescribing complex patterns of atrophy. We present experiments that provide an insight into the role of different biomechanical parameters in the model. The model allows simulating images with exactly the same tissue atrophy but with different underlying deformation fields in the image. We explore the influence of different spatial distributions of atrophy on the image appearance and on the measurements of atrophy reported by various global and local atrophy estimation algorithms. We also present a pipeline that allows evaluating atrophy estimation algorithms by simulating longitudinal MRIs from large number of real subject MRIs with complex subject-specific atrophy patterns. The proposed framework could help understand the implications of different model assumptions, regularization choices and spatial priors for the detection and measurement of brain atrophy from longitudinal brain MRIs.

Introduction

Alzheimer's Disease (AD) is one of the most common types of dementia. It is a neurodegenerative disease that progresses gradually over several years with the accumulation of neurofibrillary tangles (NFTs) and amyloid-β (A-β) plaques [Braak 1991]. These microscopic neurobiological changes are followed by the progressive neuronal damage that leads to the atrophy of the brain tissue. The atrophy or the volume changes of brain tissue is a macroscopic change that structural Magnetic Resonance Imaging (MRI) can estimate in different brain regions [START_REF] Dubois | Revising the definition of Alzheimer's disease: a new lexicon[END_REF]].

There is no treatment of AD so far, partly because the exact mechanisms of the disease are not well known. Nevertheless, there has been several clinical trials and disease modifying drug development efforts in the past three decades [Schneider 2014]. Since the external symptoms appear several years after the changes seen in imaging [Frisoni 2010], longitudinal images can play an important role in the development of disease modifying drugs. So far, structural MRIs have primarily been used for estimating local volume changes in individual AD patients; these measurements have been used to formulate hypotheses on the temporal dynamics of AD.

An interesting alternative avenue consists in modeling the tissue loss process in order to compare (in a forward modeling setting) different hypotheses for the prediction of patient-specific time series MRIs. The ability of developing realistic individual models of brain shape changes to predict patient-specific longitudinal MRIs can have far reaching consequences. For instance, the patient specific AD trajectories predicted by the model could be useful in monitoring drug effects in AD patients by comparing them against the observed brain changes.

It is nevertheless very challenging to develop a comprehensive model that can predict realistic synthetic time series of MRIs following AD patient's trajectory. Modeling neurodegeneration is a complex task requiring a hierarchy of models accounting respectively i) for how and where neuronal death occurs, ii) for its effects on brain shape changes, and iii) for the subsequent brain appearance in longitudinal MRI. In Figure 1 we show a breakdown of this complex process in three major modeling blocks which represents, at a very high level, the comprehensive modeling and simulation of realistic longitudinal MRIs in AD. The first block abstracts the multi-scale models of neuronal death at the cellular level into a macroscopic map of how the atrophy spreads spatially and evolves temporally at each voxel of the brain MRI.

Figure 2.1: High level systems diagram for modeling and simulation of longitudinal MRIs in AD patients. Spatial and temporal distribution of neuronal deaths is represented in Atrophy Generation block which causes the brain shape changes represented in Brain Deformation block. This deformation along with the MRI acquisition conditions variability result intensity change in time series structural MRI of AD patients. The error in predicted follow-up from the actual observed follow-up MRI could also be used to optimize for the parameters of the developed models using a feedback system as shown above.

Knowing the patterns of local neuronal deaths and local volume loss is just one aspect of the problem; we also need to model the consequences of neuronal loss on brain shape changes. This is represented by the block Brain Deformation in Figure 2.1. We believe that biomechanics of brain tissue does play an important role in the way brain's shape change as a result of local volume loss, and this topic is going to be one of the main subjects of this paper.

Finally, time-series of structural MRIs capture the brain shape changes but also contain additional noise, partial volume effects and image acquisition artifacts. This is also an important aspect to consider when modeling and simulating the appearance of change in longitudinal MRIs for AD patients. This part is shown in Realistic MRI generation. Furthermore, a proper optimization framework might also be necessary to estimate the patient specific parameters of the models if we are to perform model personalization. This is represented by a feedback loop in Figure 2.1.

A number of atrophy simulators [Smith 2003[START_REF] Camara | [END_REF], Karaçali 2006[START_REF] Pieperhoff | Detection of structural changes of the human brain in longitudinally acquired MR images by deformation field morphometry: methodological analysis, validation and application[END_REF], Sharma 2010] have been proposed in the literature. These simulators address either just the Brain Deformation or both the Brain Deformation and Realistic MRI generation blocks in Figure 2.1. They propose different methods to simulate time-series images with a desired volume change. All of these simulators were developed with the objective of evaluating atrophy estimation algorithms. We can broadly distinguish two major approaches used in such simulators: Jacobian based, and biomechanical models.

In Jacobian based methods [Karaçali 2006[START_REF] Pieperhoff | Detection of structural changes of the human brain in longitudinally acquired MR images by deformation field morphometry: methodological analysis, validation and application[END_REF], Sharma 2010], the desired level of atrophy is set at each voxel, and the deformation that best approximates the prescribed level of atrophy is found. Optimization of the deformation involves regularization to enforce the smoothness of the transformation and topology preservation. These simulation approaches have a number of limitations, which prevent their use and generalization in modeling oriented applications. The main issues that we identified are the following:

Plausibility and interpretation. The modeling assumptions and the regularization parameters of the energy minimization cannot be easily linked to the biophysical and mechanical process of tissue deformation. The choice of certain regularizations such as topology preservation can also have some undesirable side effects such as making it difficult to simulate the opening up of sulci.

Spatially varying tissue properties. Brain tissue and CSF are considered to respond to the volume change with the same law which is not the case in reality. Indeed, while neuronal loss in brain tissue is a gradual process, the CSF is replaced three to four times with the production of about 500 ml to 600 ml per day [Damkier 2013]. Jacobian-based approaches with uniform tissue properties are thus limited to explore questions such as: do different brain regions such as brain stem, cerebellum, cortex etc. respond with physical deformation in the same way to the neuronal deaths and local volume loss? Can we have parameters with a physical meaning for different brain tissue types that change the deformation we get even for exactly the same atrophy pattern? If tissues respond differently to the same amount of volume loss in brain, these models cannot accurately model the resulting shape changes and on the appearance of time-series MRIs unless the regularization is made spatially varying.

Skull invariance. In AD the brain deforms but the skull is rigid and hence the deformation model should not allow skull to move. The skull invariance is not imposed in [Karaçali 2006]; In [Sharma 2010], as the authors show, imposing skull invariance results in larger error in the obtained Jacobian near the skull. Since the cortical surface lying near the skull is an important area for AD, it is desirable not to have error in the obtained Jacobian in these areas. Finally, when only volume loss is prescribed, as seems to be the case in the evaluation experiments of [Sharma 2010, Sharma 2013], it is not clear which regions of the brain expand to compensate for the volume loss since the volume within the skull must be constant when skull invariance is imposed. The spatial distribution of the resulting non-zero error in the desired vs. obtained Jacobian map is not easy to control in this case. Biomechanical models generate tissue deformation based on biomechanical principles. As far as we know, the only model proposed so far for AD application other than the one we present here was a thermoelastic one [Smith 2003[START_REF] Camara | [END_REF]. In this thermoelastic model, one defines the volume changes in particular structures and tissues of a meshed brain by assigning different thermal coefficients. Thermoelastic model of tissue deformation is solved using Finite Element Method (FEM) to obtain a deformation field. To simulate time series of images, the deformation field interpolated from the mesh to input baseline image is used. An important limitation of this method is that it requires estimating regional thermal coefficients based on the desired volume changes which makes it difficult to prescribe complex voxel-wise atrophy patterns accurately. Although different tissue types can be differently modeled by considering tissue-specific values of thermo-elastic constants, the meaning of these parameters are difficult to link to the AD process. Furthermore, the variability of the resulting brain deformation depending on the choice of the tissue-specific parameters has not been investigated. Finally, FEM involves moving back and forth from voxels of patients MRI to reference labeled 3D mesh which creates numerical difficulties and inaccuracies in the model personalization.

In [START_REF] Khanal | A Biophysical Model of Shape Changes due to Atrophy in the Brain with Alzheimer's Disease[END_REF]] we proposed a proof of concept for a new biomechanics-based tissue loss model that addresses the limitations of the previous simulators discussed above. This biophysically plausible model of brain deformation due to atrophy is constrained to fit a prescribed atrophy rate at each voxel of the parenchyma. In this work, after analyzing in detail the modeling assumptions, we provide a thorough derivation of the mathematical formulation and of the numerical implementation. There is evidence that endogenous mechanical forces at the cellular level influence brain structure and function [Tyler 2012]. Although the detailed mechanisms of these interaction still deserve further investigation [Tyler 2012, Mueller 2015], it is clear that they play a role at the macroscopic level which is the scale where we observe changes in the structural MRIs.

Our model thus emphasizes, for the first time, the link between model parameters and assumptions on the real mechanobiology of neurodegeneration at a macroscopic scale. The presented experiments provide a better insight on the role of different biomechanical parameters of the model, and show that different assumptions about the atrophy process can lead to different deformations even for the same input tissue atrophy. Furthermore, we use the proposed model to study the interrelationship of various spatial atrophy patterns and how they affect the image intensity appearance. The proposed model implements the Brain Deformation block of Figure 2.1 and provides a mathematically solid and flexible framework to allow the future implementation of more complex modeling assumptions about neurodegeneration in the Atrophy Generation block.

Concerning the realistic MR image generation block, previous works in [START_REF] Camara | [END_REF], Sharma 2010, Sharma 2013] provide an interesting framework for adding different kinds of intensity noise on the simulated images for the benchmarking of atrophy measurements tools. Even though this is a desirable component of a generic atrophy simulation tool, intensity noise accounts only for a small part of the variability of atrophy measurement tools. Indeed, it has been shown that the largest variability in the atrophy measurements is due to the individual variability of the brain anatomy and atrophy pattern, as well as to the wrong modeling hypothesis [Sharma 2013, Rohlfing 2006]. Therefore, in this work we focus on the development of a framework that can exactly prescribe any complex pattern of atrophy in order to simulate a wide range of patient specific brain changes.

In the following section we present the detailed assumptions and the development of our biophysical model of brain deformation due to atrophy. Section 2.3 provides the implementation details and describes how follow-up images with any desired atrophy can be simulated from any input brain MRIs. Section 3.3 shows some examples of such simulations. It also presents experiments that provide an insight into the role of different model parameters on the model outputs. In Section 2.5, we study how local and global atrophy estimation algorithms perform when a same amount of global volume changes are prescribed in two completely different ways: uniform volume changes exclusively in gray matter vs. uniform volume changes exclusively in white matter. We also present qualitative analysis of the impact of varying model parameters on the results of local atrophy estimation method for the same prescribed atrophy. Section 2.6 presents a framework to evaluate atrophy measurement algorithms and is illustrated by assessing the atrophy measurements in various brain structures by using representative segmentation based and registration based estimation algorithms.

Biophysical Model of Brain Deformation

The proposed model is based on a series of basic assumptions motivated by the following anatomical and biophysical notions. Human brain is enveloped by a set of membranes called meninges and the CSF that lies between the skull and the brain. Due to neurodegeneration in AD, the brain shrinks along with some structural readjustment of the parenchyma. This process can be seen as a deformation of the brain parenchyma along with its mass and volume loss. The CSF volume increases to compensate the tissue volume loss while the skull remains rigid without any deformation. It is important to note that the CSF production is at a much smaller time-scale (hours) compared to the tissue atrophy (months).

In the following sections we explicit the mathematical details of the tissue loss model based on these basic assumptions.

Impact of Loss of Volume on Conservation Equation

When modeling the deformation in AD, classical continuum mechanics formulation cannot be directly applied because conservation of mass does not hold due to the presence of atrophy. To model the brain shape changes due to neuronal deaths in AD, it is important to take note of the observations in longitudinal brain MRIs of AD patients and also compare against other neurodegenerative diseases. In diseases like Creutzfeldt-Jakob disease, no gross brain shape changes are reported and the imaging only shows hyperintense signals on T2-weighted images [Johnson 1998]. However, this is not the case in AD and longitudinal MRIs show a remarkable decrease of brain volume instead [Frisoni 2010] without any "holes". That means the tissues should restructure as the neuronal deaths increase with time. This leads us to a basic assumption in the proposed model that after the death of neurons, remodeling of the tissue occurs such that the tissue density remains constant while both the mass and volume decrease. This assumption of incompressible material but with mass loss leads us to the conservation law given by (see Appendix B for the derivation):

∇ • u = -a (2.1)
where ∇ • u is the divergence of a displacement field u associated to the deformation of the brain during a period of time t and a = a(x) is the volume loss per unit volume at position x during this time.

Constrained Minimization of the Elastic Energy

To explicitly model the neuronal loss and tissue remodeling at the microscopic level, one requires biochemical and cellular physiological knowledge in detail. These mechanisms along with the spatial and temporal evolution of the cell loss are not well known for AD. The proposed model abstracts the phenomenon that evolves during several months or years in the brain at a macroscopic scale. It is based on the assumption that atrophy creates internal stress which results in the deformation minimizing a strain energy. In other words, the brain parenchyma deforms with the prescribed atrophy by minimizing the strain energy. Using Saint Venant-Kirchoff model for an elastic material, this can be expressed as the minimization of:

R(u, p) = µtr(E(u) 2 ) + λ 2 (tr(E(u))) 2 -p (∇ • u + a) (2.2)
where p is a Lagrange multiplier, µ and λ are Lamé constants, and E is Langrangian Green strain defined as: E = 1 2 ∇u + ∇u T + ∇u T ∇u . By taking a sufficiently small time step ∆t, this deformation could be reasonably modeled as being linear elastic. For example, for a 2% global atrophy rate per year, we have ∆t = 1 year, and the actual atrophy after one year is a = 0.02. This linear elastic assumption is done for a small time step only because remodeling occurs to eliminate the internal stress induced by the death of neurons.

Under linear elastic assumptions, minimizing the energy in equation (2.2) is equivalent to solving the following system of equations.

µ∆u -∇p = (µ + λ)∇a

∇ • u = -a (2.3)
where ∆u is a component-wise Laplacian of u. This system of equations is very similar to the Stokes flow equation in fluid dynamics [Batchelor 2000]. The difference is in the non-zero divergence term which corresponds to the loss of matter in the tissue. The right hand side of the first equation of this system can be seen as a force term. That means the gradient of the prescribed volume loss acts as the force term, f = (µ + λ)∇a, that moves the tissue for the structural remodeling. The Lagrange multiplier p can be interpreted as a virtual pressure whose algebraic values can be seen as sources and sinks of matter. The second equation of this system of equations has a mass source term, -a.

Modeling CSF Region

The timescale of CSF production is hours [Damkier 2013], which is much smaller compared to the time scale of tissue remodeling due to atrophy. Thus the CSF should be allowed to expand as required when the brain deforms due to the prescribed atrophy. This expansion should automatically adjust and compensate for the total loss of volume prescribed in the parenchyma. For this, we release the strict constraint present in the second equation of the system of equations (2.3) by re-introducing the variable p as: ∇•u+kp = 0, where k is the compressibility with units of Pa -1 . Now, the pressure in the CSF adapts to the expansion of CSF as required to compensate the prescribed volume loss in the parenchyma. Furthermore, since all of the CSF is considered as the fluid circulating in the brain and being constantly produced (for the timescale of months), the notion of structural readjustment due to the internal stress is not relevant. Thus the force term in system of equations (2.3) can also be set to zero. This leads us to the following system of equations for the CSF region:

µ∆u -∇p = 0 ∇ • u + kp = 0 (2.4)
If we take the gradient on both sides of the second equation above, we have ∇p

= -(1/k)∇ (∇ • u).
Replacing ∇p in the first equation of this system of equations, we get: µ∆u + (1/k)∇ (∇ • u) = 0. This is the same equation as that of elasticity but with no external force and with k = 1/(µ + λ). The CSF deformation is constrained to be such that the total volume change is equal and opposite to the total prescribed tissue volume change, and the CSF-tissue interface has a continuity in u and p. How different regions of CSF contribute to compensate for the global volume loss depends on the choice of the value of k. This is further discussed in Section 2. 4.4 The combined systems of equations for both regions are as follows:

Brain Parenchyma µ∆u -∇p=(µ + λ)∇a ∇ • u = -a CSF µ∆u -∇p=0 ∇ • u + kp =0 (2.5)
Boundary conditions: As the skull is rigid and no deformation occurs in the skull, we set Dirichlet boundary conditions with zero displacement at the skull. This means that there is no flux of matter across the brain-skull interface. This interface is the skull boundary that completely encloses the CSF and brain parenchyma. In other words, setting u = 0 at the boundary ∂Ω results in ∂Ω u • n ds = 0 where n is the vector normal to the surface ds of the boundary. Using divergence theorem we can re-write this surface integral as a volume integral Ω ∇ • u dΩ = 0. Thus, with the Dirichlet boundary conditions we set, the sum of integrals of ∇ • u in the parenchyma and the CSF should equal zero. Since we constrain the divergence in the parenchyma region with the prescribed atrophy, the system will find u in the CSF region such that integral of ∇ • u over the CSF region is opposite to the integral of ∇ • u in the parenchyma. This is how the CSF expands to compensate the volume loss prescribed in the tissue.

Material Parameters µ and λ: The Lamé parameters are related to the material's Young's modulus (E) which describes its response to mechanical stress in the following way:

λ = νE (1 + ν)(1 -2ν) , µ = E 2(1 + ν) (2.6)
where ν is a Poisson's ratio.

The first Lamé parameter λ does not have a direct physical meaning but is related to the compressibility. However, there is a strict incompressibility constraint with ∇•u = -a in our model. The system adapts the value of p based on the chosen value of µ and the input a. Unlike in standard elasticity, λ does not appear in the coefficient term of u and is only present in the force term (µ + λ)∇a. It weights the stress exerted due to the gradient of the prescribed atrophy. Its impact on the solution is explained with experimental results in Section 2.4.4.

The second Lamé parameter µ, also known as shear modulus or modulus of rigidity, is related to the stiffness of the material. However there is no consensus on the stiffness of the brain with widely varying estimates in the range of 0.1 to 16 kPa [Tyler 2012]. Furthermore, the stiffness of the brain tissue is reported to reduce in normal aging [Sack 2011] and particularly in AD [START_REF] Murphy | [END_REF]]. In [START_REF] Nogueira | [END_REF]], authors estimate equivalent mechanical stress related to brain atrophy in AD by using strains computed from the deformation of the brain in longitudinal images and with µ = 2.2 kPa taken from [START_REF] Murphy | [END_REF]]. The effect of choosing different values of µ in the proposed model is further explained in 2.4.4.

Implementation of the Biomechanical Tissue Loss Model

This section describes the major implementation steps required to simulate synthetic follow-up images from an input baseline brain MRI using the proposed model of brain deformation.

Skull Stripping and Brain Segmentation

The model minimally requires segmentation of at least three regions of the input MRI:

i. Brain parenchyma (Gray matter/White matter region)

ii. CSF region iii. Skull and outside region Skull stripping separates the brain region from the skull and other outside regions. This enables us to impose Dirichlet boundary conditions of zero displacement in the skull and outside regions of the input image. Similarly, the segmentation of the brain into GM/WM and CSF enables us to numerically solve the system of equations (2.5). Any skull stripping algorithm, and any segmentation algorithm that can extract GM/WM and CSF can be used. The choice of algorithms/softwares used for skull stripping and segmentation in this work are detailed in later sections.

Input Prescribed Atrophy Map

To simulate the desired atrophy, a voxel-wise atrophy map in the GM/WM region must be provided as the input to the model. This corresponds to a in equation (2.3). This atrophy map is thus a scalar image with desired values of volume changes in GM/WM regions, and zeros in all other regions. This image should be in the same space as the input MRI. The model can be initialized with any desired atrophy pattern, either at the regional level, or voxel-wise. Regional atrophy, such as in the whole brain, in specific brain tissues, or in regions of interests (ROIs) such as the hippocampi, can be prescribed by using any reliable anatomical segmentation tools. Otherwise, desired atrophy patterns can be generated at the voxel level by using other brain morphometry methods such as non-linear image registration (tensor based morphometry).

In the following sections we explicitly detail the brain morphometry tools chosen for each experiment. These computational methods are among the state-of-art instruments for brain morphometry, and are freely available on the web.

Staggered Grid Discretization and Finite Difference Method

For computing the deformation field that would warp the input image, equation (2.5) must be solved numerically. The computational domain for this equation is obtained from the input MRI using skull stripping and segmentation as described above. We use Finite Difference Method (FDM) with staggered grid discretization to solve the system of PDEs (2.5). Using staggered grid with proper placing of the pressure and velocity variable ensures stability in the solution. FDM is chosen instead of FEM to avoid brain meshing and the complexity of transporting computed variables from mesh to image at each iteration. This allows us to solve the system in a grid that is of the same size as the input image where the grid fits naturally to the image. This also makes it easier to obtain the partition of the computational domain into different regions directly from the skull stripping and segmentation step.

For typical brain MRIs of 1mm 3 resolution, this computational problem size becomes so large that direct solvers are impractical due to memory limitations. The system of equations (2.5) is similar to Stokes flow equation which is a saddle point system [START_REF] Benzi | [END_REF]]. It needs a suitable combination of an iterative solver and a preconditioner to solve it. We use a Schur factorization to split the equations into the momentum equation and the pressure equation. Each of these equations is solved using different iterative solvers. Our implementation uses composable solvers for multiphysics with PETSc library [Balay 2013] using fieldsplit preconditioner, an approach detailed in [Brown 2012] with an example for Stokes flow solver with Schur complement factorization. The momentum equation is preconditioned with hypre which is an algebraic multigrid preconditioner and can be called from the PETSc interface. The implemented system is run using distributed computing in a locally available cluster.

Using the Model to Explore Different Atrophy Patterns

Simulating Images with Any Desired Atrophy Maps

Figure 2.2 shows a simulation example where a 20% hippocampal atrophy is prescribed in addition to a uniform global 2% atrophy. For this case, we used ROBEX [Iglesias 2011] for skull stripping; FSL FAST [Zhang 2001] to segment the skull stripped image into GM, WM and CSF regions; and FSL FIRST [Patenaude 2011] to segment hippocampi and amygdalae to obtain the region of enhanced atrophy. This is a simple example to illustrate that any desired atrophy can be prescribed at voxel level once the desired segmentation of the input image is obtained. Solving the model using equation ( 2

Simulating Realistic Atrophy Patterns

We can also simulate more realistic atrophy patterns in different brain structures by taking the values reported in literature. For instance, in Figure 2.3 the prescribed atrophy is derived from a table in [Carmichael 2013] that reports a mean two year atrophy of amnestic MCI patients in 35 different cortical regions. Bottom row of the figure is the computed atrophy values from the obtained deformation field. As should be the case, the divergence map is the same in the tissue while in CSF the region expands to compensate for the volume loss. We used the FreeSurfer segmentation tool [Fischl 2002] for the whole brain segmentation and to extract the regions of interests (ROIs) used in [Carmichael 2013]. The whole brain segmentation with recon-all step of FreeSurfer includes skull stripping. However, FreeSurfer does not segment the sulcal CSF and only ventricular regions of CSF are segmented. Once the FreeSurfer segmentation is obtained, the sulcal CSF regions are added using the following approach:

Step 1. Binarize FreeSurfer segmentation image into brain tissue vs background.

Step 2. Get a distance map of the binary image which approximates the Euclidean distance of the foreground (tissue) in number of voxels.

Step 3. In the FreeSurfer segmentation image, label as sulcal CSF the voxels that correspond to a distance less than the threshold β in the distance map.

The chosen value of β corresponds to the width of CSF region desired outside of the cortical surface. In the following experiments, β was set to 2 voxels. [Carmichael 2013]. The computed atrophy map is the negative of the divergence obtained from the solution of the model when solved for the prescribed atrophy. This is shown in the bottom row. The negative values of the computed atrophy map (in blue) corresponds to the expansion of CSF to compensate the volume loss prescribed in the brain tissues. We can see here that the resulting atrophy of the parenchyma is exactly the one prescribed at the voxel level while the expansion of the CSF is not uniform.

Simulating Large Atrophy With Multiple Time-steps

The proposed model allows simulating very large atrophy. In Figure 2.4, we see several slices a baseline and the simulated follow-up obtained by prescribing large atrophy of 90% in Hippocampi and cortical gray matter, and 40% atrophy in white matter. The results in the figure shown was simulated with a single time-step. Usually such large atrophies are not so common. However, we could use the model to simulate large atrophies that span over several years. In this case, the linear model with single time-step may not capture the associated deformation realistically. We can incorporate nonlinearity by solving the system of equations of the model multiple times and composing the deformation field obtained at each time-steps as explained below:

• Get the displacement field u 0 by solving the model using the initial atrophy map a 0 , baseline image I 0 , and the segmented label image L 0 as input.

• For each time step t = 1 to n:

-Warp a t-1 , I t-1 and L 0 using u t-1 • u t-2 ... • u 0 to get a t , I t and L t respectively.

-Solve for u t using a t , I t and L t as input.

In Figure 2.5, we simulate six time steps starting from an initial atrophy map measured from two MRIs of AD patient 2 years apart with the baseline age of 60 yrs.

Although conceptually there is no limitation in this method to simulate very large atrophy over very long periods of time, there are nevertheless technical difficulties.

If the material parameters µ is set to be discontinuous with big jump in the values of µ in different regions of the brain, the solver converges very slowly. Thus, for these cases of highly discontinuous parameters, it can become challenging to simulate large number of time steps. Similarly, when simulating multiple time steps, the atrophy map and the label image must be updated by warping them with the displacement field of the previous time step. Since the label images are required to use nearest neighbor interpolation, this can result in some of the tissue atrophy values to leak into the nearby CSF regions during the warping of the atrophy map. So, for the next step, we modify the atrophy map by redistributing uniformly all the non-zero atrophy from these CSF voxels to the nearest tissue voxels in the 3x3 neighborhood.

We qualitatively assessed the results shown in Figure 2.4, by asking to an expert neurologist of the Nice Resource & Research Memory Centre (Nice, France), to assess the plausibility of the simulated atrophy progression. According to the clinical evaluation, the pattern of morphological changes shows realistic CSF expansion at multiple scales (ventricular enlargement and sulcal widening) and cortical thinning, along with a plausible pattern of whole brain shrinkage.

Role of Different Model Parameters

Once the atrophy is prescribed and the segmentation of the input brain image is obtained, the remaining parameters that can be varied in the model are λ, µ and k. In this section we present the role of these parameters in the model. All of the numerical values presented in this work are with the following units: λ, µ, p in kPa, k in kPa -1 and u in mm.

Impact of Changing λ

In standard elasticity, once the shear modulus µ is fixed, λ is linked to the compressibility of the material as its value depends on µ and ν. However, in the presented model the deformation field must satisfy the incompressibility constraint strictly and consequently λ does not have the same usual meaning. It appears only in the equation for the brain tissue as a scaling factor to the force produced by the gradient of the prescribed atrophy. Thus, the choice of λ affects the equivalent force exerted by the gradient of atrophy and consequently the deformation field obtained from the simulation. In Figure 2.6, we show the results of varying λ for the same prescribed tissue atrophy and the same values of µ and k.

The figure shows that setting λ too large makes the deformation field unrealistically follow-up MRIs with large uniform atrophy of 90% in cortical gray matter and Hippocampus, and of 40% in white matter. These large atrophies were simulated in a single time-step in this case. We can see that the model is able to simulate realistic pattern of widening and opening of sulci, the narrowing of gyri, and at the same time the cortical surface does not move unrealistically farther away from the skull.

Impact of Changing µ

The shear modulus µ is the stiffness of the material and changing its values affects the deformation field obtained from the model. It appears on both the left hand side and the right hand side of Eq. 2.5. On the left side of Figure 2.7, we see the effect of varying µ in tissue and/or in CSF. When µ is same everywhere, changing it equally everywhere does not have a big effect on the deformation field but it scales up the pressure. However, when the value of µ is set differently in the tissue and the CSF, we see that the pattern of the deformation field also changes. In the last row of the figure, we see that the displacement fields are bigger near the regions where the input atrophy was non-zero. In this case the tissue is much stiffer than the CSF and for the same volume change the nearby tissue deforms to compensate for the prescribed volume changes. Different brain structures could have different material properties. However, there has not been a consensus on how these properties differ in the brain [Cheng 2008]. Varying the values of µ in the brain structures also produces different results for the same prescribed atrophy. For example, Figure 2.8 shows the difference in the simulated image between the same µ in all brain tissue vs. µ in brain stem 100 times more than other brain structures. It could be interesting to explore the impact of having different stiffness in various brain regions, or to optimize for these parameters by using multiple time-point images. This is further discussed in Section 4.5.

Impact of changing k

The compressibility coefficient k is always zero in the brain tissue since we have a strict incompressible constraint. In the CSF, the choice Figure 2.7: The prescribed atrophy is a uniform atrophy only in the cortex. Left: Effect of varying µ with k = 1 kPa -1 and λ = 0 kPa. µ 1 and µ 2 denote the shear modulus of the brain tissue and the CSF respectively. Right: Effect of varying k with µ 1 = µ 2 = 1 and λ = 0. Setting very small values for the compressibility k in CSF results in its uniform expansion everywhere while allowing it to by highly compressible with large values of k results in its more local expansion to compensate for the nearby tissue volume loss. The pattern of the deformation field does not change much when µ is changed in the same way in both the CSF and tissue. However, making tissue stiffer compared to the CSF makes the tissue deform differently even with the same volume loss.

of the value of k determines its ability to expand locally. On the right side of Figure 2.7, we see that setting large value of k allows the CSF to expand more locally in response to the nearby local tissue loss, while very small values of k minimizes its expansion by distributing the same value of expansion everywhere. In a certain range of the values of 1/k, in the same or up to three to four orders smaller than µ, it does not affect much the deformation pattern. This is the range where the variation of the redistribution of the CSF expansion in different parts of the brain is not large enough to have a big impact on the deformation field. However, if the value of k is made much larger, then all of the volume loss in the tissue will be compensated by the nearby CSF regions and can even overcompensate such that the CSF in the ventricles shrink instead of expanding. This is shown in the last row of Figure 2.9.

Investigating the Relationship Between Image Appearance and Atrophy Patterns

One important question investigated in this section concerns the non-trivial effects induced by a given atrophy pattern on the appearance of the simulated images. It is important because the atrophy estimation algorithms depend on the image intensity or the appearance of the follow-up images, and the proposed model is an important instrument for studying this type of question. For this purpose we initially propose a toy example to show that different atrophy patterns can lead to very similar changes in the appearance of binary segmentations of a region of interest (ROI) (2.5.1), to finally provide an application in brain atrophy detection while using state-of-art brain morphometry tools (2.5.2).

A Synthetic Example with Binary Image

We consider a hypothetical scenario in which the segmentation of the desired brain anatomical region of interest (ROI) is provided. Our model is then applied to the anatomical region by prescribing two different atrophy patterns consisting in the same amount of global atrophy, but with different spatial distribution. In formal terms, let a g be a desired global atrophy of the brain and V g be the global brain volume. If we want α i as the desired fraction of the total atrophy to be contributed by region R i with volume V i , we have:

a i = α i a g V g V i .
where a i is the uniform atrophy needed to be prescribed in the region R i .

The experimental scenario is tested on a cylindrical ROI, on which we generate two longitudinal deformations with different atrophy properties: the first atrophy pattern is concentrated in the medial axis, while the second one is more uniformly distributed in the volume.

We note in Figure 2.10 that since there is not enough texture in the input image, the resulting simulated images look very similar although the underlying volume changes are very different. This is experimentally verified by non-linearly registering the simulated follow-up images to the cylindrical baseline with the LCC-logDemons algorithm [START_REF] Lorenzi | [END_REF]]. Figure 2.11 shows that the apparent volume changes detected by the registration algorithm are very similar. Furthermore the spatial pattern of the estimated volume changes is different than both of the original input atrophy patterns. This illustrates the dependence of the registration algorithm on an implicit model (i.e. spatial regularity assumptions) when inferring the volume changes from the intensity images. Figure 2.12 shows results of a very well known atrophy estimation algorithm boundary shift integral (BSI) [Freeborough 1997] on images with same global atrophy but distributed very differently in the brain. We generated two simulated images from a single baseline MRI by prescribing global 4% atrophy either only in GM or only in WM. Figure 2.12 shows the brain edge movement for these two cases Figure 2.13: Estimations of volume changes using LCC-Demons registration algorithm for a real baseline and two different simulated follow-up images. The first and the third columns show the log Jacobian determinants of the transformation that was used to simulate the images. Both of them have the same global atrophy of 4% but distributed differently: exclusively in white matter (1st column), exclusively in cortical gray matter (third column). For each of these underlying volume changes we show results of the LCC-Demons using two different regularization schemes. Second and second-last columns show results when using harmonic regularization while the third and last columns show result when using traditional Gaussian smoothing regularization. The Log-Jacobians of the registration results show that they do not exactly match the actual prescribed Log-Jacobians but they more or less capture the underlying atrophy patterns. This is expected since the registration algorithm is unaware of the underlying model that generated the volume changes. When using the harmonic regularization (Reg1), the Log-Jacobian maps are sharper while using Gaussian (Reg2) results in more diffused maps.

obtained by running the generalized BSI [Prados 2015], part of NifTK software tools [nif 2016]. We see that the brain edge movement reported are similar in most areas of the brain although the underlying atrophy patterns that generated the follow-up are very different. So from a shape analysis perspective when looking at only GM-CSF interface, we are not able at all to differentiate between gray and white matter atrophy. This is quantitatively confirmed by the estimated volume change: in both cases we get the same overall volume loss in terms of total volume, and percentage change. The percentage brain volume change (PBVC) reported for the GM atrophy case was -2.63% while for the WM case is -2.72%. The baseline (follow-up) total volume estimated are 1032 ml (1011 ml) and 1032 ml (1000 ml) respectively. This shows as expected that gBSI can only estimate global brain volume change and cannot localize the atrophy to be in either gray matter or white matter.

Registration methods are usually used when one needs an estimate of local atrophy. Figure 2.13 shows the results of using a non-linear registration, the LCC- Figure 2.14 shows the results of LCC-logDemons from the images simulated using exactly the same pattern of tissue atrophy but with different model parameters. The figure shows that the LCC-logDemons in general finds well the underlying volume change for both sets of parameters. However, on the left part of the figure we see that it estimates non-zero volume change in the region (shown by the axis) where no underlying volume changes were actually present. Similarly we could find for any other registration algorithm different model parameters for which the registration will produce mislocalized atrophy patterns, in particular in the areas where there are less texture such as the white matter. When estimating atrophy patterns from real observed time-series images, it is not trivial to study the relationship between the real underlying unknown deformation and the estimation provided by the algorithms such as non-linear registration. The ability to simulate different images with same atrophy patterns but with different, biomechanically inspired, model parameters 2.6. Simulating Complex Patterns of Patient-specific Atrophy for the Evaluation of Atrophy Measurement Algorithms 37

allows us to study the behavior of estimation algorithms under various assumptions. This could help in the future to make more informed and biologically motivated modeling choices in the development of atrophy estimation algorithms with spatially varying regularization. The approach can be summarized in the following three steps:

Simulating

2.6. Simulating Complex Patterns of Patient-specific Atrophy for the Evaluation of Atrophy Measurement Algorithms 38

1. Measure atrophy using real baseline and real follow-up image (RBRF).

2. Prescribe the measured atrophy and simulate a follow-up image.

3. Measure atrophy again using the real baseline and simulated follow-up (RBSF).

The performance of atrophy measurement tools can thus be studied by comparing the discrepancy of the measured atrophy in the first and third step. The effect of noise and image acquisition artifacts is also an important element that could be considered while simulating the follow-up. However, this is outside the scope of the presented work, since it would require the study of reliable simulation methods to generate image artifacts such as bias field, ringing, and motion effects.

Experimental Setup

The experiments in this section uses following values for the model parameters: µ = 1 kPa, λ = 0 kPa, k = 1 kPa -1 . The value of the shear modulus µ is in the range reported in the literature [Tyler 2012]. We used FreeSurfer as the atrophy measurement tool for AMT1 in Figure 2.15. FreeSurfer is publicly available, is widely used to study longitudinal changes in different brain regions and can segment large number of cortical, sub-cortical and white matter regions of the brain. Then we made two separate measurements of the atrophy from simulated images: i) using FreeSurfer as AMT2 in Figure 2.15 ii) using LCC-logDemons as AMT2.

We used T1 structural MRI of 46 Alzheimer's patients each having multiple time-point images in the range of 2 weeks to 2 years from the Miriad dataset [Malone 2013]. For each of these 46 subjects following processing steps are performed:

Step 1. Create a subject specific template using all the available time-points. This uses longitudinal stream of FreeSurfer to create an unbiased subject specific template image [Reuter 2012].

Step 2. Get whole-brain FreeSurfer segmentation of the extremal time-point image.

The first time-point corresponds to real baseline (RB) while the last timepoint corresponds to real follow-up (RF) as shown in Figure 2.15.

Step 3. For each segmented region:

• Get the volumes reported by the segmentation in RB and RF images:

V 0 and V 1 .
• Compute the atrophy from the obtained volumes:

a r = (V 0 -V 1 )/V 0 .
This results in the RBRF atrophy map of Figure 2.15.

Step 4. Simulate follow-up image (SF) from the RB image and the RBRF atrophy map.

Step 5. Get whole-brain segmentation of SF using the previously created subjectspecific template.

Simulating Complex Patterns of Patient-specific Atrophy for the Evaluation of Atrophy Measurement Algorithms 39

Step 6. Similar to step 3, get the volume measurements in SF: V s 1 , and atrophy estimates as a s = (V 0 -V s 1 )/V 0 . This results in the RBSF atrophy map of Figure 2.15 For LCC-logDemons, we registered real baseline images with their corresponding simulated follow-up images and computed the average Jacobian determinants of the resulting deformation field in each of the ROIs provided by FreeSurfer. The volume change measure to compare against the FreeSurfer measurements is obtained by taking the mean of divergence of the stationary velocity field (SVF) obtained from the registration in each ROI. We used default parameters of the LCC-logDemons for all the subjects.

It is important to note that the scheme used for divergence computation when discretizing the grid to solve the model must be compatible to the divergence computed from the atrophy estimation algorithm. We illustrate this issue with a figure in .. for 2D. 3D case naturally extends this. The divergence computed by many registration algorithms from the image domain is not the same as the divergence we compute in the staggered grid formulation. In order to remove bias due to this difference, we use the scheme compatible to the divergence computed .. Figure 2.16: The figure shows nine slices of the real baseline, simulated follow-up and the real follow-up images of one of the subjects generated using the pipeline shown in Figure 2.15. The real follow-up image is aligned to the real baseline image using a rigid registration for visualization purpose. This allows visual comparison of the discrepancies between the simulated follow-up and the real follow-up. For this particular subject, the ventricular expansion in the simulated follow-up seems to be less than the real follow-up. It is expected that the simulated and the real followup do not match exactly since the atrophy prescribed to simulate the follow-up comes from an atrophy estimation algorithm which does not use the same modeling assumptions as our model.

Results

In [Reuter 2012], developers of the FreeSurfer longitudinal stream present test-retest reliability of the FreeSurfer segmentation by using 115 pairs of same day scans of healthy controls. The discrepancies in the volumes measured in two scans of the same patient on same day gives an idea on the variability of FreeSurfer segmentation. As a dimensionless measure of variability, they compute the absolute symmetrized percent volume change (ASPVC) of a structure with respect to the average volume:

ASPVC := 100 |V 2 -V 1 | 0.5(V 1 + V 2 )
Since we simulate follow-up images which should ideally have the same volume as their corresponding real follow-up images in the selected regions, we use the same ASPVC measure as in [Reuter 2012] to compare the FreeSurfer volume measurements on the real and simulated follow-up image pairs. Figure 2.17 shows mean and standard deviation of the ASPVC for two different sets of regions. The regions on the left are the twenty regions with lowest mean ASPVC while on the right are the same regions for which the results of the test-retest reliability study are available 2.6. Simulating Complex Patterns of Patient-specific Atrophy for the Evaluation of Atrophy Measurement Algorithms 41

in [Reuter 2012]. For the regions presented in [Reuter 2012], we find that the mean ASPVCs from our real-simulated image pairs are in the same order as that of the test-retest real image pairs. The results show that the mean ASPVCs in most regions when using longitudinal FreeSurfer stream in our real-simulated image pairs are in between the results of cross-sectional and longitudinal stream runs presented in [Reuter 2012]. However, it should be noted that the result of the test-retest study is not directly comparable to our real-simulated study because the datasets used are different; the dataset used in [Reuter 2012] is not available in the public domain. This might also have resulted in the increased variability of the computed ASPVC. For instance, a study in [Wenger 2014] shows that the FreeSurfer reliability on hippocampal volume measurements is non-uniform across different age groups.

In particular, the study shows that the volume measurements in older age groups are not as reliable as in younger groups. Since in AD patients, the structural changes are more pronounced than in normal ageing, it is possible that the reliability will be worse in AD patients compared to normal ageing. In order to ascertain this effect, further test-retest study is required with several datasets of different age groups.

Similarly, the amount and pattern of atrophy prescribed is different for each of the subjects. FreeSurfer volume measurements on the simulated images might be impacted differently when images are simulated with varying anatomy and atrophy patterns. Finally, the ventricles and other CSF regions are not constrained to have exactly the same volume change as the one measured from the real follow-up. These factors could also have increased the variability in the real vs. simulated repeat volume measurements. We can also compare the distribution of atrophy estimates of the population of AD patients from real follow-up images (a r ) with that from simulated follow-ups (a s ). The simulated follow-up images are obtained by warping the corresponding input baseline images. Thus they have two important differences from the real follow-up images: i) Image noise in the simulated images are highly correlated to the noise in their corresponding baseline images, while the noise in the real followup images are independent from the noise in the baseline images. ii) The simulated images are obtained by resampling the baseline image and hence are smoother than the real follow-up images. Both of these factors can be expected to reduce the variability of the measured atrophy in the population when using simulated followup images instead of real follow-up images. Figures 3.6 and 3.7 show the atrophy estimates in the MIRIAD dataset using real follow-ups and using the simulated follow-ups. We see that, as expected, variability in atrophy estimates in the population is reduced remarkably when using the simulated images. Most regions show a trend of underestimation of atrophy but again there are some regions such as the white matter of cerebellum in Figure 3.7, and Pallidum and cortical Cuneus in Figure 3.6 where we observe an overestimation of the underlying atrophy with FreeSurfer. The observed variability in atrophy measures of different regions in the proposed controlled scenario could also be related to the non-uniform reliability of the atrophy estimation algorithms in different brain regions. For instance, non-uniform reliability of FreeSurfer segmentation for differ-Figure 2.18: For each region, the box plots on the left show the estimated atrophy from the real follow-up images (RBRF atrophy map in Figure 2.15) using FreeSurfer, while the one on the middle and right are from the simulated follow-up images (RBSF atrophy map in Figure 2.15) using FreeSurfer and LCC-logDemons respectively. The brain structures shown are the regions on the left of Figure 2.17. We see reduced variability in the estimated atrophy of structures in the population when simulated follow-up is used. * signifies that the average atrophy of the region for the population measured from the RBSF is significantly different from the measurements obtained from RBRF(p < 0.01, two sided paired t-test). The blue * (bottom) is for LCC-logDemons while the orange one (top) is for FreeSurfer. We see reduced variablity in atrophy estimation of the population when using simulated follow-up images.

ent structures can also be seen in the result of test-retest repeat scan experiments presented in [Reuter 2012].

Discussion

In this paper we focused on the development of a brain deformation model from a known map of local volume changes, and on the study of the impact of local tissue loss on brain shape changes.

In section 2.4.4 we showed that by changing the model parameters we obtain different deformation fields even for the same input atrophy. Setting different values Figure 2.19: Same as Figure 3.6 but the brain structures shown are the regions on the right of Figure 2.17 (20 regions with least ASPVC). We again see that the atrophy estimated in the population from the simulated images has reduced variability compared to the one from real images. The atrophy estimations in real and simulated images are significantly different only for two regions when using FreeSurfer, but for nine regions when using LCC-logDemons of shear modulus in various brain regions can produce different deformations as seen from an example in Figure 2.8 where we set higher value of shear modulus in brain stem. It will be interesting in the future to study if one can optimise the values of µ to obtain even more realistic morphological changes in very specific structures as expected by the neurologists, such as the rotation of temporal poles in the coronal view or the movement of inferior part of the temporal lobe farther away from the cerebellar tentorium. If this set of stiffness parameters corroborates well with the values reported in the literature, it could provide better understanding of the response of various brain structures to the local volume loss.

The framework presented in 2.15 could also be used to calibrate the model parameters for a particular AMT. In Figures 3.6 and 3.7 we see that LCC-logDemons does not perform the same way as FreeSurfer. The result could be explained by noticing that the prescribed atrophy used in simulating images are piecewise continuous since it was computed from the segmented volumes from FreeSurfer but the model used by LCC-logDemons promotes smoothly varying Jacobian determinants. These results were obtained by using only one set of model parameters; by defining a suitable parameters optimization strategy we could obtain a minimum discrepancy in the atrophy measurements from the real and simulated follow-up images. Once the model is calibrated to a particular AMT this way, the pipeline in 2.15 can be used to study the relative bias of other AMTs.

One interesting future work concerns the optimization of the atrophy parameters to best explain the observed longitudinal images. This is akin to registration methods where one finds a best explanation of the observed changes based on some implicit models of regularization. In this case we have an explicit biophysical model of deformation with the advantages described earlier in the paper. There are two major challenges that need to be taken care of in this regard. The first one concerns the very large number of parameters: since the atrophy is prescribed in every voxel, the number of parameters equals the number of voxels in the image, or the number of considered regions of interest. The second issue is the computational time for solving the model. It requires from around thirty minutes to few hours in a locally available cluster computing resource using 80 cores (depending on the choice of model parameters, and the load in the cluster from other users) to solve the model for brain MRIs of around 1 mm resolution. This means special efforts will be required to develop the optimization framework in a computationally feasible manner for the given image size and number of selected atrophy parameters.

The experiments and results presented in section 2.6 are for the illustrative purpose of using the framework in different scenarios. These experiments are not a full fledged benchmarking of the atrophy estimation algorithms as it is not the primary focus of this paper. Extensive analyses are required for proper evaluation and to find out the discrepancies in the atrophy measurement from the real and simulated follow-up images. For instance, atrophy in each ROIs can be distributed differently such as being concentrated towards the middle of the ROIs or in the boundaries. This could impact the contrast in the edges of each of these ROIs and consequently affect the atrophy estimations. Similarly, simulated follow-up images were obtained by resampling the intensity of the baseline images and are in the same physical space as their corresponding baseline images. This way of simulating images could possibly have failed to reproduce some of the information available in the real follow-up images (such as high intensity contrast) that are used by the atrophy estimation algorithms. Answering such questions with conclusive evidence requires additional sets of experiments studying the origins of bias and variability, and the impact of different ways of simulating ground truth images on the estimation algorithms. These experiments should also be carried out on a number of additional atrophy estimation algorithms than the ones we presented in this paper and it will be the focus of one of our future studies.

In [Sharma 2013], the authors create a database of simulated images from 18 MRIs by simulating uniform hippocampal atrophies in the range of 1-14% with step size of 1%. For each pair of atrophy value and patient image, a number of simulated images are created by degrading the simulated atrophies with independent Gaussian noise. The authors estimated bias in the atrophy estimation of hippocampus using the simulated database and also developed a framework to provide confidence intervals of the atrophy estimation. The nature and magnitude of bias computed were based on the database containing simulations of the images with atrophy only in a single region. The framework presented by the author and the related database can be enriched by using the model we proposed in this work. We can simulate more images for the same prescribed atrophies in a particular region by varying the model parameters, and by varying the atrophy patterns in other brain regions. This allows studying the nature of bias due to the variation in the underlying model of deformations and in the presence of complex atrophy patterns in multiple brain regions. In [Carmichael 2013], the authors use all the cortical regions segmented from FreeSurfer to identify patterns of coordinated atrophy distributed in Gray Matter of aMCI patients. The objective of that study was to explore the distributed network account in AD by studying how different groups of cortical regions are correlated to best explain the longitudinal change. Estimating bias in the measurements from atrophy estimation algorithms in the presence of atrophy in large number of structures simultaneously can be useful in assessing results of such studies too.

Towards an Integrative Multimodal Model

The anisotropic nature of the brain parenchyma due to fibers could have an impact on the way it deforms due to atrophy. Since not much is known about this, the proposed model can be useful as it allows such an exploratory study. It has parameters µ and λ where this anisotropic information could be introduced, for e.g. from DWI images. For the same atrophy map, the effect of anisotropy on the brain deformation is an interesting question to explore.

Reliably simulating neurodegeneration due to AD and its trajectory in structural MRIs is quite challenging as we need accurate models for all three major blocks shown in Figure 2.1. The most difficult part is to generate accurate patterns of atrophy and its evolution with time. As we have seen from the examples in Section 2.5 with differential patterns of atrophy producing similar images, the atrophy estimation algorithms would benefit from a biologically motivated prior on the assumed model for regularizations. Accurate atrophy generation models require more information from other sources in addition to the structural MRIs. In addition to the research in biology pertaining to AD, perhaps a progress in other imaging modalities could also potentially provide information on the spread of imminent neuronal deaths. For instance Aβ plaques seem to occur very early at the beginning of atrophic process [Chetelat 2010]. Studies such as brain's structural connection breakdown on AD patients using Diffusion Imaging [START_REF] Stebbins | Diffusion tensor imaging in Alzheimer's disease and mild cognitive impairment[END_REF][START_REF] Daianu | [END_REF]], or functional connectivity breakdown along with the structural connectivity [Filippi 2011] could also provide better insight in the future. Similarly, there is ongoing research in developing good tracers to bind to tau proteins and to image in vivo the neurofibrillary tangles (NFT) [James 2015] in AD patients. In the future we might be able to exploit such data to propose basic hypotheses of spatial atrophy distribution using multi-modal images. This could be valuable in developing suitable models for the Atrophy Generation block.

Conclusions

We developed a biophysical brain deformation model that describes the consequence of the neuronal deaths and atrophy on the brain shape changes at macroscopic scale. The model is inspired from biomechanical principles, and treats the brain parenchyma and the CSF differently to account for the fact that the CSF is produced at a very short time-scale of hours compared to the slow process of tissue atrophy taking months. We were able to achieve different deformations of the brain even with exactly the same atrophy in the brain tissue by varying the model parameters. Any desired atrophy can be prescribed at the voxel level and simulate realistic deformations of a patient specific MRI. This ability to both i) accurately prescribe complex patterns of atrophy at each voxel, and ii) to treat different tissue types differently in accordance to their biomechanical properties, was very difficult with the previous models of atrophy simulators existing in the literature. The proposed model could be used in testing hypotheses about the distribution of brain atrophy and in exploring the interaction of mechanical response of different brain tissues to neurodegeneration. It can also be a valuable tool to understand better the interrelationship between the underlying brain deformations corresponding to specific atrophy patterns, longitudinal MRI appearance, and the bias of various atrophy estimation methods due to the modeling error.

2. Data used in the preparation of this article were obtained from the MIRIAD database. The MIRIAD investigators did not participate in analysis or writing of this report. The MIRIAD dataset is made available through the support of the UK Alzheimer's Society (Grant RF116). The original data collection was funded through an unrestricted educational grant from GlaxoSmithKline (Grant 6GKC).

3. This work benefited from the use of the Insight Segmentation and Registration Toolkit (ITK), an open source software developed as an initiative of the U.S. National Library of Medicine and available at www.itk.org.

4. The multi-platform configuration tool CMake was used for configuring ITK and facilitating its use from our project. CMake was partially funded by the U.S. National Library of Medicine as part of the Insight Toolkit project.

CMake is an open source system and it is freely available at www.cmake.org.

Introduction

Alzheimer's Disease (AD) is one of the most common types of dementia. It is a neurodegenerative disease that progresses gradually over several years with the accumulation of neurofibrillary tangles (NFTs) and amyloid-β (A-β) plaques [Braak 1991]. These microscopic neurobiological changes are followed by the progressive neuronal damage that leads to the atrophy of the brain tissue. The atrophy or the volume changes of the brain tissue is a macroscopic change that structural Magnetic Resonance Imaging (MRI) can estimate in different brain regions. Many different methods have been proposed to estimate atrophy in some particular regions of brain that are known to be affected in AD [START_REF] Dubois | Revising the definition of Alzheimer's disease: a new lexicon[END_REF]].

In addition to estimating specific brain structures with atrophy, longitudinal imaging data could also potentially be used to study the temporal inter-relationship of atrophy in different structures. For instance in [Carmichael 2013], authors estimated per-individual rates of atrophy in 34 cortical regions and in the hippocampi. Then they studied the groupings of these structures based on the correlation of the atrophy rates. In [Fonteijn 2012], authors modelled AD progression as a series of discrete events. The occurrences of atrophy in different parts of the brain were taken as different events along with clinical events. Without any prior to their ordering, the model finds most probable order for these events from the data itself. They used Bayesian statistical algorithms for fitting the event-based disease progression model. The objective of these kinds of studies is to understand how different regions of the brain interact during the neurodegeneration and find its evolution. Such studies can benefit from the availability of a large number of longitudinal images of AD patients. In this context, a model that can simulate many time-point images from a few available longitudinal images can be a valuable tool.

Atrophy simulators [Karaçali 2006[START_REF] Pieperhoff | Detection of structural changes of the human brain in longitudinally acquired MR images by deformation field morphometry: methodological analysis, validation and application[END_REF]][Smith 2003[START_REF] Camara | [END_REF]] have been proposed in the literature and used mostly for the validation of registration or segmentation methods [Camara 2007][Sharma 2010], or to estimate uncertainty in the measured atrophy [Sharma 2013].

The simulators in [Karaçali 2006[START_REF] Pieperhoff | Detection of structural changes of the human brain in longitudinally acquired MR images by deformation field morphometry: methodological analysis, validation and application[END_REF]][Sharma 2010] used a Jacobian based methods where the desired level of atrophy is specified at each voxel, and the deformation that best approximates the prescribed level of atrophy is found. Regularization is used in the optimization to enforce certain desired conditions such as topology preservation. The advantage of these methods is the ability to define atrophy maps at the voxel level. However, regularization parameters used to enforce topology preservation are generally difficult to relate to a plausible biophysical process of AD and can create difficulties in simulating opening of certain structures such as sulci. It is not trivial to consider different tissue behaviors in such approaches. In [Smith 2003[START_REF] Camara | [END_REF]], authors proposed a model of brain deformation based on thermoelasticity. Volume changes were defined in particular structures/tissues of a meshed brain by assigning different thermal coefficients. Synthetic images were simulated by warping real images with the displacement fields obtained from a Finite Element Method (FEM) based solution of the thermo-elastic model; it required the interpolation of the displacement field from the mesh to the image grid. Moving back and forth from voxels to meshes creates numerical difficulties and could introduce inaccuracies in the model personalization.

In [START_REF] Khanal | A Biophysical Model of Shape Changes due to Atrophy in the Brain with Alzheimer's Disease[END_REF]] we proposed a new biophysical model of brain deformation due to atrophy in AD that combines the advantages of the models mentioned in the previous paragraph. The mechanisms of neuronal deaths and its evolution are not well known for AD and are likely to be primarily guided by complex physiological processes. However we believe that the biomechanics of brain tissue might play an important role in determining the consequence of the neuronal deaths on brain shape changes. Our biophysical model presented in [START_REF] Khanal | A Biophysical Model of Shape Changes due to Atrophy in the Brain with Alzheimer's Disease[END_REF]] builds upon the assumptions that we relate to the biophysical process of tissue shape changes as the consequence of local volume loss. This model can be used to simulate time-series MRIs starting from a real input baseline MRI.

In this work we use our biophysical model developed in [START_REF] Khanal | A Biophysical Model of Shape Changes due to Atrophy in the Brain with Alzheimer's Disease[END_REF]] to present a framework that allows to interpolate or extrapolate patient specific unseen timepoint images from at least two available time-point images of the subject and to assess how closely these simulated trajectories follow real patient trajectories. We also improve the implementation of the boundary condition of the model by imposing zero deformation in the skull and all the regions outside of the skull. In [START_REF] Khanal | A Biophysical Model of Shape Changes due to Atrophy in the Brain with Alzheimer's Disease[END_REF]] the zero deformation was imposed at the image boundaries and not at the brain-skull boundary.

The following section briefly explains the assumptions and implementation of the biophysical model we presented in [START_REF] Khanal | A Biophysical Model of Shape Changes due to Atrophy in the Brain with Alzheimer's Disease[END_REF]], and in section 3.3 we present how we interpolate new images between two acquisition time points.

Biophysical Model of Brain Deformation due to Atrophy

The atrophy rate ã(x, t) at any position x at time t for a representative elementary volume of V (x, t) is defined as the negative rate of change of volume per unit volume:

ã = -1 V ∂V ∂t .
For any time ∆t that results in sufficiently small displacement, the amount of atrophy is a = ã∆t. Any deformation field that has atrophy a should satisfy the following equation:

∇ • u = -a, (3.1)
where u is the displacement of material particles during the ∆t.

We do not explicitly model the neuronal loss and tissue remodeling at the microscopic level which requires biochemical and cellular physiological knowledge in detail. We abstract the phenomenon that evolves over several months or years in the brain. In Creutzfeldt-Jakob disease, no gross brain shape changes are reported and the imaging only shows hyperintense signals on T2-weighted images [Johnson 1998]. However, this is not the case in AD and longitudinal MRIs show a decrease of brain volume instead without any "holes" [START_REF] Dubois | Revising the definition of Alzheimer's disease: a new lexicon[END_REF]. That means the tissues should restructure as the neuronal deaths increase with time. This leads us to a basic assumption in the proposed model that after the death of neurons, remodeling of the tissue occurs such that the tissue density remains constant while both the mass and volume decrease. We further assume that the atrophy creates internal stress which results in the deformation minimizing a strain energy.

Using Saint Venant-Kirchoff model for an elastic material, this can be expressed as the minimization of which results in the deformation of the tissue minimizing the strain energy.

R(u, p) = µtr(E(u) 2 ) + λ 2 (tr(E(u))) 2 -p (∇ • u + a) (3.2)
where p is a Lagrange multiplier, µ and λ are Lamé constants, and E is Langrangian Green strain defined as: E = 1 2 ∇u + ∇u T + ∇u T ∇u . By taking a sufficiently small time step ∆t, this deformation could be reasonably modeled as being linear elastic. For example, for a 2% global atrophy per year, we have ∆t = 1 year, and the atrophy during the year as a = 0.02.

Under linear elastic assumptions, minimizing the energy in Eq. ( 3.2) is equivalent to solving the following set of equations.

µ∆u -∇p = (µ + λ)∇a ∇ • u = -a (3.3)
where ∆u is a component-wise Laplacian of u. This equation is very similar to the Stokes flow equation in fluid dynamics. The difference is in the non-zero divergence term which corresponds the loss of mass and volume in the tissue. The momentum equation shows that the gradient of the prescribed volume loss acts as the force term that moves the tissue for the structural remodeling. The Lagrange multiplier p is some sort of virtual pressure whose algebraic values can be seen as the sources and sinks of fluid.

Modeling CSF Region

The timescale of CSF production is hours, which is much smaller compared to the time scale of tissue remodeling due to atrophy. To allow the CSF to expand as required when the brain deforms due to the prescribed atrophy, we release the strict incompressibility constraint in the system of Eqs. (3.3). Furthermore, the force term of the momentum equation in the system of Eqs. (3.3) is no longer required. Thus the combined equation for both the brain parenchyma and the CSF regions is:

µ∆u -∇p = (µ + λ)∇f ∇ • u + kp = -f (3.4)
where we have,

Biophysical Model of Brain Deformation due to Atrophy

51

• Brain parenchyma region: k = 0, and f = a

• CSF region: k = 1, and f = 0.

Boundary Conditions: Dirichlet boundary conditions with zero displacement is enforced at the skull.

Material Parameters µ and λ: The deformation model here corresponds to the structural readjustments due to cell loss, thus the Lamé parameters do not have the same usual meaning as during an elastic deformation of the material due to application of an external load/force. The voxel-wise volume change constraint and the boundary conditions i.e. the shape of the tissue-CSF and brain-skull interface have much more impact on the deformation of the brain parenchyma than any specific scalar values of µ and λ. In the present work these coefficients are set to 1 and 0 respectively.

Staggered Grid Discretization and Finite Difference Method

Eq. ( 4.1) requires a partition of the computational domain into different regions. These regions are obtained by using skull stripping and segmentation of the input baseline brain MRI. The solution of the PDE provides us a deformation field that is applied to the baseline image to generate simulated follow-up image. We use Finite Difference Method (FDM) with staggered grid discretization to solve the system of PDEs in (4.1). Using staggered grid with proper placing of the pressure and velocity variable ensures stability in the solution. FDM is chosen instead of FEM to avoid brain meshing and the complexity of transporting computed variables from mesh to image at each iteration. This allows us to solve the system in a grid that is of the same size as the input image where the grid fits naturally to the image. This also makes it easier to obtain the partition of the computational domain into different regions directly by using a skull stripping and a segmentation algorithm.

For typical brain MRIs of 1mm 3 resolution, this computational problem size becomes so large that direct solvers are impractical due to memory limitations. The system of Eqs. (4.1) is similar to the Stokes flow which is a saddle point system. It needs a suitable combination of an iterative solver and a preconditioner to solve it. We use a Schur factorization to split the equations into the momentum equation and the pressure equation. Each of these equations is solved using different iterative solvers. Our implementation uses composable solvers for multiphysics with PETSc library [Balay 2013] using fieldsplit preconditioner, an approach detailed in [Brown 2012] with an example for Stokes flow solver with Schur complement factorization. The momentum equation is preconditioned with hypre which is an algebraic multigrid preconditioner and can be called from the PETSc interface. The implemented system is run using distributed computing in a locally available cluster.

Simulating New Time-points from Two Available scans of MRIs 52

Simulating New Time-points from Two Available scans of MRIs

In this section, we propose a basic approach to simulate patient-specific atrophy pattern and generate new time-point longitudinal MRIs from two available scans. The method requires an atrophy estimator to estimate atrophy from two available timepoints. Once the atrophy is estimated, an atrophy predictor will predict the atrophy for the given patient for another time-point. This predicted atrophy is prescribed to our biomechanical model of brain deformation to produce a new synthetic image that corresponds to the input time-point. The schematic is shown in Figure 3.1 for generating a middle time-point image. By replacing the desired time-point t m by any other desired time-points, one can simulate images for different time-points.

Figure 3.1: Schematic illustrating how a mid-point synthetic image can be generated using: i) two available real time-point images, ii) an atrophy estimation algorithm, and iii) an atrophy predictor. If the real middle scan is available, the atrophy estimated from the real baseline and mid-point scan could be compared against the atrophy estimated from the real baseline and simulated mid-point scan. This approach could be adapted to interpolate or extrapolate any other desired timepoints by replacing t m by the one desired.

We use the Miriad dataset [Malone 2013] that has multiple time-point T1 structural MRIs of 46 Alzheimer's patients in the range of 2 weeks to 2 years. Since the dataset contains several time-point scans, we can compare the simulated intermediate time-point images to the corresponding real intermediate images. To prescribe personalized atrophy patterns we need an atrophy estimation for each sub-ject from the extremal time-points. We perform the whole brain segmentation using recon-all command in FreeSurfer [Fischl 2002]. For the segmentation, FreeSurfer's longitudinal stream [Reuter 2012] is used to create unbiased subject specific templates. This allows us to compare the volumes of large number of regions in the baseline and the follow-up images and estimate atrophy in each of these regions. These estimated atrophy can then be modified and prescribed to each of the baseline MRIs to predict intermediate time-point images. The setup of the experiment we performed is shown in Figure 3.2 and described as follows:

1. Find available extremal time-point scans: baseline I b and the final follow-up I f . Let t f be the time (in years) between the baseline scan and the final scan.

2. Find a mid-point scan I m that was scanned t m years after the first scan. This is found by finding t m that is closest to t f /2.

Use

FreeSurfer to estimate an atrophy map a f . This is a scalar image whose intensities are the atrophy estimated from FreeSurfer for all the segmented brain regions.

4. Simulate a follow-up image I m that corresponds to the mid-point scan I m by prescribing a m where a m = a f * t m /t f . 5. Run the FreeSurfer whole brain segmentation on this simulated image I m and compute volumes of all the segmented regions.

6. Compare FreeSurfer computed volumes of all regions of the images I m and I m .

Results and Discussion

Figure 3.3 shows the atrophy estimates for all the patients using FreeSurfer segmentation in all the regions that were used in [Carmichael 2013]. Since the volume changes estimated are for the population of AD patients, as expected, almost all of these regions have median of the volume change ratio less than zero. But, we can also see that there is variation in the volume changes of each of the regions in the population. This variation includes both the natural variation occurring in the population and the variation due to the variability of FreeSurfer. For example, the figure also shows that there are cases where volume expansion is reported in the cortical regions which are unlikely to be occurring in reality in the AD patients. For all the 46 patients, we also have the FreeSurfer estimates of volume changes in the simulated midpoint images in all the regions used above. Figure 3.4 shows the difference between volume change estimates in the real and simulated images. We can see that for most of the regions, the median difference is close to zero. However, there are also few regions where the difference is in the order comparable to the measured volume change ratio itself. In [Reuter 2012], absolute symmetrized percent volume change (ASPVC) was used to measure the test-retest reliability of the FreeSurfer segmentation by using 115 pairs of same day scans of healthy controls, where

ASPVC := 100 |V 2 -V 1 | 0.5(V 1 + V 2 ) .
The discrepancies in the volumes measured in two scans of the same patient on same day gives an idea on the variability of FreeSurfer segmentation. In our case, we compute ASPVC by considering V 1 and V 2 to be the volumes measured by FreeSurfer on the real and predicted mid-point images. Figure 3.5 shows mean and standard deviation of the ASPVC (from 46 patients) in two different sets of brain structures reported in [Reuter 2012]. Figures 3.6 and 3.7 show box plots of the estimated atrophy from the real midpoint images and predicted mid-point images side by side for all the regions for which ASPVC were shown in Figure 3.5. From the figures, variability of the atrophy in the population seems to be reduced when using predicted images compared to the real ones.

Higher variability in the difference seems to be mostly in the regions where there are higher variability in the atrophy estimates of the real mid-point images. Larger inter-subject variation of the difference between the atrophy estimate in the real mid-point image and the interpolated mid-point image could be due to several rea- [Reuter 2012] by computing FreeSurfer volume measurements on same day test-retest scan pairs of 115 healthy controls. The bar plots in orange, cross-run-mean and long-run-mean, are mean error of the Reuter regions shown on the right. cross-run-mean were computed using FreeSurfer cross sectional stream, while long-run-mean were computed using FreeSurfer longitudinal stream. ASPVC errors for the real vs. simulated mid-point images are in the same order as that of the mean error for real baseline repeat scans. It is worth noting that the dataset used in this work is not the same as the one used in [Reuter 2012], which is not publicly available. sons. One obvious issue is that the FreeSurfer segmentation with the longitudinal stream expects all the images that are to be segmented to be preprocessed in the same manner. However, in our case the interpolated mid-point image has undergone an extra resampling step while the real mid-point image has not. This extra resampling step is required because the interpolated mid-point image was obtained by warping the real baseline image with a displacement field. Furthermore, the choice of interpolation during the resampling step can also affect the volume measurements by FreeSurfer. We used trilinear interpolation for the resampling. The extra resampling step and the choice of interpolation does have an effect on the estimation of volumes. This has been shown, for instance in [Sharma 2010] for other segmentation based atrophy estimation techniques. Furthermore, the interpolated mid-point image has a noise (noise inherent in any MRI) that is highly correlated with the real baseline image. However, the noise in real mid-point image is not correlated to the baseline image. This also affects the atrophy estimation and hence contributes to the variability in the atrophy estimation difference. A detailed analysis must be done to find out the regions that are the most reliable ones in estimating volume changes for both the real and simulated images. The performance of the atrophy measurement tools on simulated images should be thoroughly evaluated to find out the best regions that we can rely upon to test how closely we predict volume changes in new time-point images.

In this case we have interpolated the intermediate time-point by linearly scaling the estimated atrophy. For a small time window of a couple of years this is reasonable but if we want to extrapolate for instance for several years we would need a non-linear model of atrophy progression. The presented framework allows one to compare the trajectory of brain shape changes with different models of atrophy progression. The ability to prescribe any desired atrophy at any time point allows one to introduce atrophy at different regions of brain at different times. This can be exploited in evaluating the methods proposed in studies such as [Fonteijn 2012] which order the events from time-series data.

Conclusions

We proposed a framework to generate patient specific multiple time-point images based on our biophysical model of brain deformation due to atrophy in AD. The used model is motivated from biomechanical principles and it models the consequence of tissue loss in brain shape changes. From the available two scans of MRI of a patient at two different time-points we estimated atrophy in large number of brain structures using FreeSurfer whole brain segmentation [Fischl 2002]. The derived atrophy patterns were linearly scaled and prescribed to the biophysical model to simulate the another time-point image. Using the MIRIAD dataset [Malone 2013] of 46 AD subjects with multiple time-points we compared the simulated time-point images against the actual time-point images. The future works will include building the most reliable methods to compare the volumes in simulated and real images. We will also explore the possibility of evaluating methods that study the temporal relationships, ordering and co-evolution of atrophy in different structures of the brain. Abstract This paper presents a simulator tool that can simulate large databases of visually realistic longitudinal MRIs with known volume changes. The simulator is based on a previously proposed biophysical model of brain deformation due to atrophy in AD. In this work, we propose a novel way of reproducing realistic intensity variation in longitudinal brain MRIs, which is inspired by an approach used for the generation of synthetic cardiac sequence images. This approach combines a deformation field obtained from the biophysical model with a deformation field obtained by a non-rigid registration of two images. The combined deformation field is then used to simulate a new image with specified atrophy from the first image, but with the intensity characteristics of the second image. This allows to generate the realistic variations present in real longitudinal time-series of images, such as the independence of noise between two acquisitions and the potential presence of variable acquisition artifacts. Various options available in the simulator software are briefly explained in this paper. In addition, the software is released as an opensource repository. The availability of the software allows researchers to produce tailored databases of images with ground truth volume changes; we believe this will help developing more robust brain morphometry tools. Additionally, we believe that the scientific community can also use the software to further experiment with the proposed model, and add more complex models of brain deformation and atrophy generation.

Introduction

Structural Magnetic Resonance Imaging (MRI) has been widely used for invivo observation of morphological changes over time in human brain. Atrophy or tissue volume loss measure from structural MRI is an established biomarker for neurodegeneration [START_REF] Dubois | Revising the definition of Alzheimer's disease: a new lexicon[END_REF]]. There is a large number of brain morphometry algorithms developed in the literature which estimate global or local atrophy from structural MRIs [START_REF] Wright | [END_REF], Freeborough 1997, Ashburner 2000, Smith 2002, Hua 2008]. Volume/atrophy measurements obtained from such algorithms have been used to test various clinical hypotheses about neurodegenerative diseases [START_REF] Wright | [END_REF][START_REF] Sepulcre | [END_REF][START_REF] Koch | [END_REF]. Similarly, comparison of different neurodegenerative diseases have also been performed based on these measurements [START_REF] Rosen | [END_REF], Whitwell 2005]. Since atrophy estimation is an inverse problem, the estimation algorithms require a model with certain parameters. The results obtained from such algorithms depend on model assumptions and the parameters used. Often, these assumptions are implicit and cannot be directly linked to the biophysical process of neurodegeneration. For instance, tensor based morphometry (TBM) encodes local volume changes by computing Jacobian determinants of the deformation field obtained from non-linear registration of longitudinal MRIs [Ashburner 2015]. Such methods contain model biases because TBM results depend on the choices of regularization used during the registration of images [START_REF] Ashburner | Symmetric diffeomorphic modeling of longitudinal structural MRI[END_REF]]. Estimating and correcting the bias present in such morphometry tools is important, especially for clinical applications.

In addition to tracking volumetric changes in specific brain structures, longitudinal imaging data can also be used to study the temporal inter-relationship of atrophy in different structures. For instance, [Carmichael 2013] studied the groupings of 34 cortical regions and hippocampi from the per-individual rates of atrophy estimates in these regions. In [Fonteijn 2012], authors defined AD progression as a series of discrete events. Along with other clinical events, the timings of atrophy in various brain structures were included in a set of discrete events. Without any prior to their ordering, the model finds the most probable order for these events from the data itself. They used Bayesian statistical algorithms for fitting the event-based disease progression model. The objective of these studies were to understand how different regions of brain evolve during the neurodegeneration.

In this context of increasing use of the atrophy measurements from longitudinal MRIs in testing or discovering clinically relevant hypotheses, it is important to study the bias and variability of the atrophy estimation algorithms. The actual volume changes in real longitudinal MRIs are not known. Thus, the evaluation and validation of atrophy estimation algorithms require generating images with known volume changes, called ground truth images.

A number of atrophy simulators have been proposed in the literature to produce ground truth MRIs [Smith 2003[START_REF] Camara | [END_REF], Karaçali 2006[START_REF] Pieperhoff | Detection of structural changes of the human brain in longitudinally acquired MR images by deformation field morphometry: methodological analysis, validation and application[END_REF], Sharma 2010, Khanal 2016b]. These simulators use a model that attempts to produce a deformation field with the specified volume changes in the input brain MRI. To produce realistic scenarios of noise and acquisition artifacts, some of these simulators also use a model to produce noise and artifacts in the simulated image.

Such simulators have been used for the validation of registration or segmentation based atrophy estimation algorithms [START_REF] Camara | [END_REF][START_REF] Pieperhoff | Detection of structural changes of the human brain in longitudinally acquired MR images by deformation field morphometry: methodological analysis, validation and application[END_REF], Sharma 2010], to estimate the bias in such algorithms, and also to estimate uncertainty in the measured atrophy [Sharma 2013]. These studies have estimated the bias by simulating simple atrophy patterns in a small number of brain regions, or uniform diffused global atrophies. However, real case scenarios could have a much more complex atrophy distribution occurring in many brain structures at the same time.

Noise and imaging artifacts have an important impact on the results obtained from atrophy estimation algorithms [START_REF] Pieperhoff | Detection of structural changes of the human brain in longitudinally acquired MR images by deformation field morphometry: methodological analysis, validation and application[END_REF][START_REF] Camara | [END_REF], Sharma 2010]. Thus, proper evaluation of atrophy estimation algorithms by using simulated ground truth images requires simulation of realistic variation in noise and intensity too. All the previous atrophy simulators have warped the input baseline image with the deformation field obtained from a model of brain deformation. Then, extra noise and artifacts are added on this warped image by using another artificial model. The intensity noise in structural MRIs has been shown to be governed by a Rician distribution where the noise is Gaussian in k-space [Gudbjartsson 1995]. Thus the Rician noise can be added in the simulated images as follows:

• Use two independent random variables following zero-mean Gaussian distribution to compute the real and imaginary parts of a complex number at each voxel.

• Add the magnitude of this complex number to the intensity of the simulated MRI at each voxel.

For example, [Sled 1998] used this approach to add noise in simulated MRIs that were used for the validation of intensity bias correction scheme they presented. Using the same approach, [START_REF] Camara | [END_REF]] added noise to the simulated ground truth images with atrophy. In addition to the Rician noise described above, other noise and artifacts have also been shown to affect the measurements of atrophy estimation algorithms [Sharma 2010[START_REF] Pieperhoff | Detection of structural changes of the human brain in longitudinally acquired MR images by deformation field morphometry: methodological analysis, validation and application[END_REF]]:

• Bias field inhomogeneity arising due to poor radio frequency (RF) coil uniformity.

• Geometrical distortions that are present due to the errors in gradient field strength and non-linearity of gradient fields in the MR scanner [Langlois 1999].

• Interpolation of intensities during various pre-processing steps of TBM based analysis framework (e.g., resampling of the images into a common template space).

Many other acquisition artifacts may not be simulated because we do not have faithful models. In this work, we develop a new framework to simulate longitudinal images with specified volume changes, and also the realistic intensity variations including the noise and acquisition artifacts. In the framework, we use our biophysical model of brain deformation [Khanal 2016b] to obtain a dense deformation field with specified volume changes. To obtain the realistic intensity variations, we use an approach introduced by [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF]] where the authors simulate visually realistic time series of cardiac images. Intensity variation in the simulated images of a patient is obtained by resampling the intensities from real images of the same patient taken at different times. Figure 4.1 shows a diagram of this framework. To implement this framework, we have developed an open-source atrophy simulator software called Simul@trophy.

Section 4.2 explains all the blocks of the framework shown in Figure 4.1. Starting from a small set of real scans, we show how longitudinal images with different atrophy patterns and realistic intensity variations can be simulated. Section 4.3 shows some simulation results using Simul@trophy, and also illustrates some potential applications of the simulator. In Section 4.4, we present some example simulations to illustrate some of the important points to consider when using Simul@trophy for different applications, such as evaluation of atrophy estimation algorithms, validation of data-driven disease progression models, training of brain morphometry algorithms based on machine learning etc.

Simulating Realistic Longitudinal Images with Atrophy/Growth 63

Figure 4.1: Pipeline to simulate synthetic images using Simul@trophy. Starting from a real baseline image of a subject, synthetic images with known volume changes can be generated. These synthetic images can follow intensity characteristics of either the input baseline or other images of the same subject. Pre-processing is required to generate an atrophy map and a segmentation image, which are fed as inputs to the brain deformation model. For a given set of parameters, the model computes a velocity field whose divergence is equal to the prescribed atrophy map at each voxel of the regions selected by using the segmentation image. Intensity simulator uses the output field to produce synthetic image whose intensity is resampled either from the input real baseline or from any other image as desired.

Simulating Realistic Longitudinal Images with Atrophy/Growth

We use the biophysical model presented in [START_REF] Khanal | A Biophysical Model of Shape Changes due to Atrophy in the Brain with Alzheimer's Disease[END_REF], Khanal 2016b] to generate dense deformation field with specified complex patterns of volume changes. This deformation field is then used to generate realistic synthetic longitudinal images with intensity variation, noise and artifacts, just like in real longitudinal images.

The major components of the simulation framework, as seen in Figure 4.1, are: i) Pre-processing ii) Brain deformation model iii) Realistic intensity simulator.

Simulating Realistic Longitudinal

Images with Atrophy/Growth 64

Pre-processing to generate a segmentation image and atrophy maps

A pre-processing step takes a real scan of a patient as an input baseline image, and generates the required inputs of the brain deformation model: a segmentation image and a specified atrophy map.

Segmentation Image

There are three labels in the segmentation image used by Simul@trophy (Figure 4.1):

• Label0: regions where no deformation should be prescribed,

• Label1: regions where certain volume changes are prescribed (the values of volume changes are provided with an input atrophy map),

• Label2: regions where the deformation model is allowed to adapt volume changes as required to compensate for the total volume change prescribed in regions with Label1.

Pre-processing usually starts with a brain extraction that excludes the skull and outside regions (also called skull stripping). Skull stripping is followed by a segmentation such that each voxel of the input image could be assigned to one of the three labels. For example, a typical pre-processing step that includes a segmentation of brain parenchyma and CSF would produce a segmentation image with the following labels:

• Label0: Skull and outside regions of the input image • Label1: Gray and white matter regions • Label2: CSF regions

Atrophy map

An atrophy map is a scalar image with desired values of volume changes in Label1 regions of the segmentation image, and zeros in all the other regions. It is defined at each voxel as follows:

a = V 0 -V 1 V 0 ,
where V 0 and V 1 are the volumes of the material lying in a voxel at time t 0 and t 1 respectively. Thus, regions with volume loss have positive values of a while the regions with volume expansion have negative values of a. An example atrophy map is shown in Figure 4.1.

In this work, we illustrate example simulations where two kinds of pre-processing steps were used to generate the atrophy maps:

Segmentation based atrophy map The user can set uniform values of atrophy in regions of interests (ROIs) of the 4.2. Simulating Realistic Longitudinal Images with Atrophy/Growth 65 brain. In this case, one must first perform a segmentation of all ROIs in which a non-zero value of atrophy is desired. Then, it is straightforward to create a scalar image having intensity values taken from a table, which contains the labels of ROIs and the corresponding desired atrophy values.

Registration based atrophy map The results of longitudinal non-rigid registration can be used to estimate local volume changes, for instance by computing Jacobian determinants of the displacement fields or by computing the divergence of the stationary velocity fields obtained from the registration. These local volume changes obtained from the registration based methods are usually smoothly varying in space and can be used to prescribe either:

• smoothly varying atrophy maps,

• or atrophy maps uniform in ROIs obtained by averaging, in each ROIs, the atrophy obtained above. Figure 4.2 shows two such atrophy maps with very different patterns, but having the same average regional volume changes.

A biophysical model of brain deformation with prescribed volume changes

Simul@trophy uses the biomechanics based model of brain deformation detailed in [Khanal 2016b]. The model abstracts the phenomenon that evolves during several months or years in the brain at a macroscopic scale. It is based on the assumption 4.2. Simulating Realistic Longitudinal Images with Atrophy/Growth 66 that atrophy creates an internal stress which results in the deformation minimizing a strain energy. In other words, the brain parenchyma deforms with the prescribed atrophy by minimizing the strain energy. By taking a sufficiently small time step ∆t, this deformation could be reasonably modeled as being linear elastic. For example, for a 2% global atrophy rate per year, we have ∆t = 1 year, and the actual atrophy after one year is a = 0.02. For a given segmentation image, the model yields a deformation field with the prescribed atrophy at each voxel of Label2 regions (e.g. brain parenchyma). Label1 regions (e.g. the CSF) will correspondingly adapt its volume to globally compensate for the prescribed volume changes in the Label2 regions. For a single time-step, the displacement field u is obtained by solving the system of Eqs 4.1, where Dirichlet boundary conditions of zero deformation are prescribed in Label0 regions.

Regions with: Label0

u = 0 Dirichlet boundary conditions      Label1 µ∆u -∇p=0 ∇ • u + kp =0      Label2 µ∆u -∇p=(µ + λ)∇a ∇ • u = -a      (4.1)
The system of Eqs. 4.1 shows that the incompressibility constraint is relaxed in Label1 regions, while it is strictly satisfied in Label2 regions. The impact of the choice of different values for the model parameters µ, λ and k are detailed in [Khanal 2016b]. For the same prescribed volume changes, we can obtain different deformation fields by varying these model parameters. In this work, we focus on generating ground truth images with known volume changes and not necessarily generating the exact evolution of the AD patients. Hence, we set the model parameters as follows unless specified otherwise: µ = 1 kPa, λ = 0 kPa, k = 1 kPa -1 .

Once the field u with the prescribed volume changes is obtained from the model as described above by using an input baseline image I b , we can simulate a synthetic follow-up image I s as follows:

• Let y = Φ sim (x) = u + x describe a mapping of a point x in physical space to another point y by applying the transformation corresponding to the dense deformation field Φ sim , or the displacement field u.

• Let Φ sim ⋆ I b describe an action of the diffeomorphism Φ sim on the image I b . Thus, the new synthetic image I s , obtained by warping I b with the deformation field Φ sim is given by:

I s = Φ sim ⋆ I b = I b • Φ -1 sim .

Adding realistic intensity variation to synthetic longitudinal MRIs

In realistic scenarios, longitudinal MRIs are taken at multiple scan sessions often with slightly different acquisition parameters or even with different scanners. For generating more realistic synthetic longitudinal MRIs, variations in intensity and noise present in real longitudinal MRIs must also be simulated. If multiple repeat scans of a subject are available, we can use them to simulate such variations in synthetic longitudinal sequences. Assuming that all the available scans of the subject are already aligned using affine registration, this section explains the proposed method of adding realistic variations in the intensity characteristics.

Starting from an input baseline image I b 0 of a subject, the previous sections explained how we can obtain a deformation field Φ sim from the brain deformation model, and use it to simulate a follow-up image

I s 0 = Φ sim ⋆ I b 0 .
I s 0 has the same intensity characteristics as I b 0 , and the intensity noise in I s 0 is strongly correlated to the noise present in I b 0 .

If I b 1 is another scan of the same subject taken on the same day, we can obtain a new simulated image by resampling the intensity from I b 1 , but still using the same Φ sim :

I s 1 = Φ sim ⋆ I b 1
The realistic variation of intensity and artifacts present between the two real scans I b 0 and I b 1 are now also present between the real baseline image I b 0 and the simulated follow-up image I s 1 .

The above approach assumes that the brain has not gone any morphological changes between the scan sessions of the two real images. If the scan time-points of the two images are too far apart to have this assumption valid, we can no longer directly apply Φ sim to the second image. Let I r be another real scan of the patient taken at a time later than that of the baseline image I b 0 . There might be some morphological changes (e.g. atrophy) in I r compared to I b 0 .

To simulate a new synthetic image with the same atrophy as that of I s 0 but with the intensity resampled from I r , we must first perform a non-rigid registration between I r and I b 0 . If Φ reg is the deformation field obtained from the non-rigid registration between I r and I b 0 , it can be used to get an image Φ reg ⋆ I r which is aligned to I b 0 . In the ideal case, Φ reg ⋆ I r and I b 0 are perfectly aligned with the only differences lying in the intensity characteristics and the noise.

We can now compose the deformation fields Φ sim and Φ reg to generate a new synthetic image as follows:

I s 2 = (Φ sim • Φ reg ) ⋆ I r .
I s 2 has the same atrophy as that of I s 0 but with the intensity characteristics of I r . Figure 4.3 illustrates how we obtain I s 0 , I s 1 and I s 2 . These three simulated images have the volume changes as encoded by Φ sim , but have intensity characteristics coming from three different real images of the same patient.

Figure 4.4 illustrates how the approach described in this section can be used to generate multiple sets of longitudinal simulated sequences having identical morphological evolution but different variations of intensities. The three shaded regions in Figure 4.4 are the sets of longitudinal sequences with identical volume changes but with different variations of intensities.

Simulation Examples with Simul@trophy

This section presents simulation examples of synthetic longitudinal MRIs with prescribed atrophy patterns and realistic intensity variations1 . The real input MRIs used for the simulations presented in this section were obtained from the publicly available OASIS dataset [START_REF] Marcus | [END_REF]]. All these real MRIs had undergone intensity inhomogeneity correction using ANTs -N4BiasFieldCorrection [Avants 2011], and had been transported to a common space using affine registration with FSL -FLIRT [Jenkinson 2001]. More precisely, these images had undergone the Pre-Processing and Position Correction steps of the Longitudinal Log-Demons Framework (LLDF) detailed in [Hadj-Hamou 2016].

Figure 4.5 shows a simulation example where uniform atrophy patterns are prescribed in the hippocampi, the gray matter (GM), and the white matter (WM) regions. The ventricles and sulcal CSF regions are allowed to expand as required to compensate for the volume loss in the brain parenchyma. The figure shows two simulated images whose intensities are resampled from two different images: i) the input baseline image I b ii) another follow-up image of the same subject, I r . The figure also shows intensity histograms of these two simulated images for a selected ROI. The selected ROI is a 2D WM region where the simulated images do not have a distinct morphological changes from I b . Thus, the differences in the intensity histograms of I b and the simulated images for this ROI is mostly due to the variation in intensity characteristics of the different images. We can see from the figure that the intensity characteristics of the simulated image resampled from I b closely matches the intensity characteristics of I b . And resampling the intensity from a different image I r of the same subject allows simulating realistic variation of intensities.

To simulate multiple time-point images, the following approach can be used:

• Get u 0 by solving the system of Eqs. (4.1) using the initial atrophy map a 0 and the initial segmentation image L 0 as input.

• For each time step t = 1 to n:

-Warp a t-1 and L 0 using u t-1 • u t-2 ... • u 0 to get a t and L t respectively.

-Solve for u t using a t and L t as input.

Once all the deformation fields Φ s i corresponding to u i for i = 0, 1, ..., n are obtained, these deformation fields can be used as shown in Figure 4.4 to simulate different sequences of longitudinal images. In Figure 4.6, a simulation example of two longitudinal sequences each having three new time-point images is shown. Both sequences were simulated by prescribing a smoothly varying atrophy pattern. The smoothly varying atrophy pattern prescribed in this example is more complex than the simple pattern used in the previous example. It is adapted from the divergence of a stationary velocity field obtained by performing LCC log-Demons registration [START_REF] Lorenzi | [END_REF]] of the input baseline image with a follow-up image of the same subject. The first sequence consists of all the images whose intensities are resampled from the same input baseline image I b , while the second sequence consists of the images whose intensities are resampled from different real MRIs of the same subject. Thus, as shown in Figure 4.7, the first sequence does not have the realistic variation of intensities while the second sequence has the realistic variation of intensities. With this example, we also illustrated that we can generate multiple sequences of longitudinal images with same atrophy patterns but different variations of intensities.

Figure 4.8 shows a simulation example where we prescribe growth instead of atrophy in the brain tissue. The prescribed atrophy in this case is the negative of the atrophy map prescribed in Figure 4.6. From the segmentation image shown in Figure 4.8, we can see that the ventricles were allowed to adapt the volume changes as required to compensate for the volume changes in the brain parenchyma. From the three simulated time-points, we can see that these ventricles are shrinking and the brain parenchyma regions are expanding. The example shows that Simul@trophy can be used to simulate images of not only future time-points, but also the past time-point images.

In Figure 4.9, we show an example where synthetic sequence of images is simulated by starting from a baseline image of a healthy subject. However, the pre- scribed atrophy is derived from an atrophy estimated from the AD patient used in Figure 4.6. The input baseline images of both the AD patient and the healthy subject were segmented using FreeSurfer [Fischl 2002]. In all the segmented regions including the white matter parcellations of the AD patient, the average values of the smoothly varying atrophy map were computed. These regional average values of the atrophy computed from the AD patient were then transported to the corresponding regions of the healthy subject. Thus, in Figure 4.9, we can see that the prescribed atrophy is region-wise uniform instead of smoothly varying. For comparison, the figure also shows three real time-point images of the healthy subject along with the three simulated time-point images with atrophy derived from the AD patient. It has also been very widely used in a very diverse set of applications that also include the medical field. It is a very powerful library that supports wide range of iterative solvers and preconditioners for large systems of equations. The solvers implemented in PETSc can scale very well to large distributive computer systems. Simul@trophy runs from command lines where the required inputs and optional

Impact of registration on simulated images

In Section 4.2.3, we explained that starting from an input baseline image of a subject, I b , we can generate two synthetic images:

I s 1 = Φ sim ⋆ I f and I s 2 = (Φ sim • Φ reg ) ⋆ I f
where Φ sim is the deformation field obtained from the brain deformation model using I b as the input baseline image, and Φ reg is the deformation field obtained from the non-rigid registration between I b and a real follow-up image I f . Perfect alignment of the two images with a non-rigid registration is possible only in the ideal case scenario.

In such an ideal case, the simulated images I s 1 and I s 2 have identical shapes of the brain structures with the only differences lying in the intensity characteristics. In practice, this is almost never the case, and we present below an example of the impact of registration result on the simulated images. Let us use the following short notations for various images described in this section. The figure also shows two simulated follow-up images SF_in_RB and SF_in_RF that are resampled from (RB) and (RF) respectively. We can see that in most regions of the brain, the two simulated images have almost identical morphological appearances. However, there are also regions such as 2 and 5, where the morphological appearances of the two simulated images are not identical. From the registration results for these regions 2 and 5 in the zoomed patches, we can see that the registration is also not accurate in those regions.

Figure 4.10 illustrates the impact of registration result Φ reg on the simulation results. The figure shows both the registration and simulation results along with zoomed patches of RB, RB_to_RF, SF_in_RB and SF_in_RF. As expected, SF_in_RB and SF_in_RF have different intensity characteristics coming from RB and RF respectively. In the regions where registration is accurate, the two simulated images look almost identical except for the differences in the intensity characteristics. However, in the regions where registration is not accurate enough, SF_in_RB and SF_in_RF do not have identical shapes as expected. Thus, for the proposed method of using deformations obtained by registration for simulation, it might be preferable to use aggressive non-linear registrations with a much bigger weight given to similarity terms than the regularization terms.

Discretization scheme for the divergence computation

In [Khanal 2016b], a standard staggered grid discretization was used for solving the system of Eqs. (4.1). The discretization scheme is shown in Figure 4.11 in 2D for illustration; explanation on 2D extends naturally to 3D. In the figure, we can see that the components of the displacement field variable u lie on cell faces and not at cell centres. However, all the input and output images for the model, including the output displacement field image, are standard images that have their values lying in cell centres or voxels. Our implementation of the solver internally creates the required staggered grid for the given input images. Once u is computed within the solver of system of Eqs. (4.1), its values at cell faces are interpolated to obtain the values at cell centres which are then assembled to send as output displacement field image. Within the solver, the numerical scheme used for the discretization of ∇ • u = -a is:

u i+1/2,j,k -u i-1/2,j,k h x + v i,j+1/2,k -v i,j-1/2,k h y + w i,j,k+1/2 -w i,j,k-1/2 h z = a i,j,k (4.2)
where,

u =   u v w   .
Simul@trophy then provides output displacement field image with the values of u lying at cell centres or voxels by using linear interpolation as follows:

    u i,j,k v i,j,k w i,j,k     =     u i+1/2,j,k + u i-1/2,j,k /2 v i,j+1/2,k + v i,j-1/2,k /2 w i,j,k+1/2 + w i,j,k-1/2 /2     (4.3)
To compare divergence maps of this output field with the ones obtained from tools external of Simul@trophy, the only accessible values are the interpolated ones. ITK is widely used in registration based brain morphometry algorithms, but the default derivative computation of ITK has the following centred difference stencil:

u i+1,j,k -u i-1,j,k 2 * h x + v i,j+1,k -v i,j-1,k 2 * h y + w i,j,k+1 -w i,j,k-1 2 * h z = a i,j,k (4.4)
Replacing the components of u at cell centres from Eq. 4.3, we get,

u i+3/2,j,k + u i+1/2,j,k -(u i-1/2,j,k + u i+3/2,j,k ) 4 * h x + ... = a i,j,k (4.5)
The scheme in Eq. (4.5) does not match the one that was used internally by Simul@trophy shown in Eq. ( 4.2). This results in discrepancy if we compare input prescribed atrophy maps against the externally computed divergence maps ∇ • u. Thus, in this work, we have added an implementation for the scheme in Eq. (4.5) so that users can choose either of the two possible schemes of Eq. ( 4.2) and Eq. (4.5). The latter scheme is consistent with the divergence computed by the default derivative computation options of ITK. At each 3D cell, the scheme in Eq. ( 4.2) involves 6 variables of the displacement field, while the scheme in Eq. (4.5) involves 12 variables. In the rest of the paper, they will be referred to as 6-point and 12-point schemes respectively. Figure 4.12 shows the error in specified vs. obtained atrophy when using the two different numerical schemes. As expected, we can see that when a consistent numerical scheme is used, there is no difference between the specified and obtained atrophy. When the schemes are not consistent, the error is larger on the areas where the prescribed atrophy values change sharply.

If the simulated ground truth images using Simul@trophy are used for the evaluation of atrophy estimation algorithms, one must also be careful about the measure of volume change used in addition to the numerical scheme used. For instance, many TBM based brain morphometry algorithms use Jacobian determinants as a measure of volume change. To compute ground truth volume changes of the simulated images for the evaluation of such algorithms, users should compute Jacobian determinants of the output displacement fields u obtained from Simul@trophy by using the same numerical scheme as used by the atrophy estimation algorithm being evaluated.

Figure 4.12: Error due to non-consistent numerical schemes in Eq. ( 4.2), and Eqs. (4.4 and 4.5). ∇ • u shown in the figure are computed external of Simul@trophy by using the default ITK derivative computation scheme shown in Eq. (4.4). When this divergence computation is consistent with the one used in Simul@trophy, we should obtain zero error with ∇ • u + a = 0. This is indeed the case, as seen on the right, when we use 12-point stencil of Eq. 4.5. We see non-zero errors when using 6-point stencil from Eq. ( 4.2) because this scheme and the default ITK scheme are not consistent. The figure shows that the error gets larger at areas where prescribed atrophy has discontinuous jumps.

Implementation of image warping

When implementing an algorithm to warp an image with a given deformation field, it is more convenient to use the inverse of the deformation field. If Φ s is the output deformation field obtained from the brain deformation model by using I b as the input baseline image, Φ s maps any point x in I b to a point y in the simulated image I s as follows: y = Φ s (x).

However, y is not guaranteed to be a discrete voxel location. Since we do not know the intensity values of I s a priori in the nearby discrete positions, the problem of interpolation is much more complex. Thus, we start from a discrete voxel location y in I s where the value of intensity is to be found. Then, the corresponding position x in I b can be obtained by using the inverse deformation field:

x = Φ -1 s (y).
If the transformed point x is not a discrete point, we can interpolate the intensities of I b from neighbouring discrete locations. Let us denote the interpolation by square brackets. Thus i = I[x] describes a mapping of a point x to an intensity, i, of the MR image I at x. Using this notation, the intensity of the simulated image at any position x is given by: I b Φ -1 s (x) . We could avoid the inversion of the deformation field obtained from the model, if we provide negative values of the actual desired atrophy map. For instance, if we want a 10% volume loss in hippocampus, we can instead prescribe 10% volume expansion and solve the model to obtain a deformation field. Assuming this field as already being inverted and warping the image will simulate an image with 10% volume loss in the hippocampus as desired.

The simulator can be used to choose whichever method the user prefers by using the following option:

--i n v e r t _ f ie l d _ t o _ w a r p # Invert u; default : do not invert

The implementation of the inversion is adapted from a fixed-point scheme implementation available in ITK [Luethi 2010]. By default, the simulator uses B-spline interpolation of order three to warp the input images.

Standalone utility tools and scripts for pre-processing and post-processing

There are some standalone tools and scripts available for various pre-processing and post-processing operations that are detailed in the documentation of the released software. Some of these tools for pre-processing and post-processing operations are C++ executables based on ITK, while others are python scripts. In this work, all the input segmentation of the model were obtained by using FreeSurfer. As explained in [Khanal 2016b], these segmentation maps were processed to obtain in the format required by the model. Although the provided scripts are developed for FreeSurfer segmentation maps, they can be easily modified to adapt to other pre-processing tools. Finally, the registration and simulation deformations were composed using ComposeMultiTransform of Advanced Neuroimaging Tools (ANTs) [START_REF] Avants | [END_REF]].

The core component of Simul@trophy is the implementation of the brain deformation model. Resampling of the intensity is straightforward once the deformations from the model and from registration are available. The simulator is not dependent on any one particular registration algorithm. Although we used LCC-LogDemons for illustrative purposes, this can be replaced with any other non-rigid registration algorihtms. Similarly pre-processing is also independent of Simul@trophy. We used FreeSurfer in the simulation examples shown in this work, but any other skull stripping and segmentation algorithms can be used. Simul@trophy provides some example scripts and some utility scripts, which could be modified when using other tools for the pre-processing step.

Discussion

In [Khanal 2016b], we presented a method to generate a subject-specific atrophy pattern by first measuring the atrophy from the available time-points, and then simulating a new time-point by prescribing the measured atrophy. In [Khanal 2016c], we extended the method to interpolate an unavailable intermediate time-point MRI. In this work, we added realistic variation in the intensity of the synthetic images.

The simulation examples were shown using three types of atrophy patterns: i) very simple uniform volume changes in small number of regions, ii) uniform atrophy in large number of regions, and iii) smoothly varying atrophy patterns. For each subject, we could generate large number of synthetic images by perturbing these atrophy patterns in different ways. Even with the same atrophy pattern, we can generate multiple sets of longitudinal sequences of varying intensity characteristics using the approach illustrated in Figure 4.4. Thus, by changing the atrophy patterns and the image intensities, Simul@trophy could be used to generate a database of very large number of simulated images. Such a database might be useful for training of machine learning algorithms.

Simul@trophy can be used in evaluating atrophy estimation algorithms in similar ways as done by [START_REF] Pieperhoff | Detection of structural changes of the human brain in longitudinally acquired MR images by deformation field morphometry: methodological analysis, validation and application[END_REF][START_REF] Camara | [END_REF], Sharma 2010]. The ability to prescribe atrophy at any time point allows the user to introduce volume changes at different regions of the brain at different times. Thus, another interesting application of the simulator is to train and/or validate disease progression models such as the models proposed in [Chen 2012, Fonteijn 2012, Jedynak 2012, Dukart 2013, Schmidt-Richberg 2016]. Having a database of longitudinal MRIs with known spatio-temporal distribution of atrophy can be useful to validate such algorithms. Furthermore, since the algorithms use a data driven approach, the simulator could be useful to train or fine-tune such models.

Another possible application is in filling up unavailable time-point MRIs of some of the subjects, when performing group-wise longitudinal analysis. In such studies, usually the available time-point images of each subject are used to estimate subjectspecific volume changes. These subject-specific measurements are then used to perform group-wise statistics to check whether there are significant differences amongst different groups in some particular regions of the brain. Databases used in such analyses, might not always have all the required time-point images for all the subjects. This could lead to bias if all the subjects are not aligned properly in the temporal dimension of disease progression. Simulating new time-point images for some subjects and using them in the analysis might allow evaluating the impact of such mis-alignments.

We hope to promote two directions of research in the community with opensource release of Simul@trophy. First, the public availability of Simul@trophy enables researchers to build their own simulated databases as needed. This might also hopefully lead to a large public database of ground truth simulated images, that could be used for benchmarking and evaluation of various image based morphometry tools. Second, we hope that Simul@trophy allows other researchers to build upon the biophysical model we presented in [Khanal 2016b], and investigate further, providing more accurate models of brain atrophy.

Finally, Simul@trophy is general enough to be used for other imaging modalities such as CT scans. It could also be used with images of any other organs, where one requires simulating specified volume changes. In this case, the pre-processing should be changed accordingly to generate a segmentation image and atrophy maps. Thus, once the software is public, other researchers might find it useful in applications that we have not foreseen yet.

Conclusions

We proposed a simulation framework that can generate realistic longitudinal MRIs with specified volume changes. The framework allows generating large number of subject-specific multiple time-point images based on a biophysical model of brain deformation due to atrophy. We developed an open-source software Simul@trophy to implement the proposed framework. The major part of Simul@trophy is the implementation of our brain deformation model presented in [Khanal 2016b]. Simul@trophy is based on widely used state of the art libraries PETSc (for solving large systems of equations) and ITK (for medical image processing). Since the software is publicly available in an open-source repository, we hope that researchers can use it to create databases of ground truth images. The framework could be used to generate a common public database, which in turn could be used to validate and evaluate a large number of available atrophy estimation algorithms. Similarly, these databases could be valuable for data driven disease progression models including machine learning algorithms. Validation and training of the models that study temporal relationships, ordering and co-evolution of atrophy in different structures of the brain could be another interesting application.

National Library of Medicine and available at www.itk.org.

4. The multi-platform configuration tool CMake was used for configuring ITK and facilitating its use from our project. CMake was partially funded by the U.S. National Library of Medicine as part of the Insight Toolkit project.

CMake is an open source system and it is freely available at www.cmake.org.

Appendix

4.7.1 Running Simul@trophy from command lines

Once the pre-processing steps described in Section 4.2.1 are performed and the desired atrophy map is generated, these images can be used as input to the model by providing the following command line arguments:

-atrophyFile # Input atrophy map -maskFile # Input segmentation file -imageFile # Input image file

If the model parameters µ and λ have uniform values in Label1 and Label2, they can be provided as an argument to the option -parameters. On the other hand, if they need to have different values in different parts of the brain, one needs to provide them as images similar to other input images as shown below:

-parameters #µ,λ in Region1 , Region2 . Format : µ1, µ2, λ1, λ2 -muFile # Ignore µ from -parameters , use this image -lambdaFile # Ignore λ from -parameters , use this image --useTensorLambda #λ given as DTI ; default is scalar image Some of the important options available are:

-b ou nd ary _c on dit io n # di ric hl et _at _w al ls or dir ic hl et_ at _sk ul l --div12pt_stencil # Use 12 -point scheme ; default : 6 -point scheme --relax_ic_in_csf # Region1 :

∇ • u + kp = 0; default is ∇ • u = -a -relax_ic_coeff # Value of k -numOfTimeSteps # Number of time -steps to solve for
To solve the system of Eqs. (4.1), the argument to -boundary_condition should be dirichlet_at_skull and -relax_ic_in_csf must be provided. Using dirichlet_at_walls instead of dirichlet_at_skull will consider regions with label0 in the same way as the regions with label2, and sets the Dirichlet boundary conditions only at the image borders.

If -numofTimeSteps is greater than one, the simulator provides an output displacement field obtained by composing output displacement fields of each time-steps. For any time-step n < numOfTimeSteps, it also provides output synthetic image by warping the input image with the displacement field obtained by composing output displacement fields from time-step 1 to n. In addition to these outputs, if desired, some other extra outputs can be generated as shown below:

Appendix

84

-resPath # Result path to store all the results -r e s u l t s F i l e n a m e s P r e f i x # Prefix to be provided to all the images --writePressure # Write p as image to disk . --writeForce # Write (µ + λ)∇a as image to disk . --writeResidual # Write solver residual as image to disk .

Chapter 5

Conclusions and Perspectives

Conclusions

The thesis envisioned, at a very top level, a framework for developing a comprehensive biophysical model that could predict and simulate realistic longitudinal MRIs of patients with Alzheimer's Disease (AD). As shown in Figure 5.1, the framework includes three major building blocks: i) Atrophy generation ii) Brain deformation iii) MRI generation. Within this framework, we developed a biophysical model of brain deformation that represents the Brain deformation block (Chapter 2, [START_REF] Khanal | A Biophysical Model of Shape Changes due to Atrophy in the Brain with Alzheimer's Disease[END_REF], Khanal 2016b]). Similarly, inspired from the work of [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF]], we implemented the MRI generation block that allows generating realistic variation of intensity in synthetic longitudinal images (Chapter 4, [Khanal 2016a]). Finally, we also provided a primitive approach to implement the Atrophy generation block, which allows generating subject-specific atrophy patterns from the available timepoint MRIs (Chapter 2, 3, [Khanal 2016c, Khanal 2016b]).

The simulator software developed during the thesis will be released as opensource soon. Making the software open-source and available publicly, we hope that other researchers can use it:

• To create ground truth images for evaluation and validation of atrophy estimation and disease progression models.

• To create database of large number of realistic longitudinal MRIs to train machine learning algorithms.

• As a starting point to further study various hypotheses about spatio-temporal evolution of atrophy and its impact on brain shape changes.

• As a starting point to develop more realistic models of atrophy generation by combining multi-modal imaging and other biomedical information. The error in predicted follow-up from the actual observed follow-up MRI could also be used to optimize for the parameters of the developed models using a feedback system as shown above.

Perspectives

Creating Databases of Synthetic Longitudinal MRIs for Machine Learning Applications

In [Sharma 2013], the authors created a database of simulated images from 18 MRIs by simulating uniform hippocampal atrophies in the range of 1-14% with a step size of 1%. For each pair of atrophy value and patient image, a number of simulated images were created by degrading the simulated atrophies with independent Gaussian noise. Using the simulated database, the authors also estimated the bias in the hippocampal atrophy measurements and also developed a framework to provide confidence intervals on the atrophy estimation. The nature and magnitude of the bias computed were based on the database containing simulations of the images with atrophy only in a single region. The framework presented by [Sharma 2013] and the related database can be enriched by using Simul@trophy to simulate a large variety of images with different atrophy patterns and intensities.

The framework we developed during the thesis could be used to generate a common public database for the validation and evaluation of the available atrophy estimation algorithms. Such a database could also be valuable for data driven disease progression models including machine learning algorithms. Validation and training of the models that study temporal relationships, ordering and co-evolution of atrophy in different structures of the brain could be another interesting application.

Optimisation of Model Parameters

One interesting future work concerns the optimisation of the atrophy parameters to best explain the observed longitudinal images. This is an inverse problem akin to registration methods where one finds a best explanation of the observed changes based on some implicit models of regularisation. For instance, Schweiger et al. used the inverse problem approach to estimate regional volume changes by using a thermoelastic model of brain atrophy and optimising the unknown regional coefficients of expansion [Schweiger 2005]. There are two major challenges that need to be taken care of in this regard. The first one concerns the very large number of parameters: since the atrophy is prescribed at each voxel, the number of parameters equals the number of voxels in the image, or the number of considered regions of interest. The second issue is the computational time required for solving the system of equations describing the model. For the brain MRIs of about 1 mm resolution, solving the system of equations requires between a few minutes to a few hours in a locally available cluster computing resource using 80 cores (depending on the choice of model parameters and the cluster load). Thus, special efforts will be required to develop the optimisation framework in a computationally feasible manner. One possible direction would be to explore the works in optimal control for: i) parameter estimation, for instance [Gholami 2016] ii) regularisation schemes for deformable registration with a constraint on the divergence of the velocity field [Mang 2015].

Towards an Integrative Multimodal Model

The anisotropic nature of the brain parenchyma due to white matter fibres could have an impact on the way it deforms due to atrophy. Since not much is known about this, the proposed model can be useful as it allows such an exploratory study. Changing the scalar parameters µ and λ to be tensors could allow introducing the anisotropic information, for e.g. from DWI images. For the same atrophy map, the effect of anisotropy on the brain deformation is an interesting question to explore.

Reliably simulating neurodegeneration due to AD and its trajectory in structural MRIs is quite challenging as we need accurate models for all three major blocks shown in Figure 5.1. The most difficult part is to generate accurate patterns of atrophy and its evolution with time. As we have shown in Chapter 2 that differential patterns of atrophy can produce similar images, the atrophy estimation algorithms would benefit from a biologically motivated prior on the assumed model for regularizations. Accurate atrophy generation models require more information from other sources in addition to the structural MRIs. In addition to the research in biology pertaining to AD, perhaps a progress in other imaging modalities could also potentially provide information on the spread of imminent neuronal deaths. For instance Aβ plaques seem to occur very early at the beginning of the atrophic process [Chetelat 2010]. Studies such as brain's structural connection breakdown on AD patients using Diffusion Imaging [START_REF] Stebbins | Diffusion tensor imaging in Alzheimer's disease and mild cognitive impairment[END_REF][START_REF] Daianu | [END_REF]], or functional connectivity breakdown along with the structural connectivity [Filippi 2011] could also provide better insight in the future. Similarly, there is ongoing research in developing good tracers to bind to tau proteins and to image in-vivo the neurofibrillary tangles (NFTs) in AD patients [James 2015]. In the future, we might be able to exploit such data to propose basic hypotheses of spatial atrophy distribution using multi-modal images. This could be valuable in developing suitable models for the Atrophy Generation block.

Lack of complete knowledge about the mechanisms of AD makes it very challenging to develop a comprehensive model able to predict accurately the real evolution of the brain shape at individual patient level. Many promising works in the past three decades have focused on the image analysis to accurately measure atrophy and to discover patterns of structural changes in the brain. This decade has seen an increased interests in developing data driven models and generative models of AD progression, and in understanding their relationships with other biomarkers of AD [De Souza 2010, Chen 2012, Fonteijn 2012, Jedynak 2012, Dukart 2013, Young 2014, Schiratti 2015, Young 2015, Ziegler 2015, Schmidt-Richberg 2016]. We believe that more effort also needs to be concentrated towards developing reliable mechanistic models able to accurately predict the longitudinal structural images from the available imaging and other relevant information. The presented work is a step forward in this direction where we have laid the foundation for a comprehensive modeling and simulation system for AD.

Using det(I + εA) = 1 + εtrA + O(ε 2 ) we can approximate J as below,

J ≈ 1 + ∆ttr (∇v) = 1 + ∆t∇ • v = 1 + ∇ • u
Now substituting J in equation (B.1), we have:

∇ • u = -a (B.2)

B.2 Minimization of Strain Energy

Second Piola-Kirchoff Tensor is said to be the derivative of energy density with respect to Green-Lagrange strain tensor. Using the notations introduced in Appendix A, we have, We want to minimize the energy. To find the gradient of the energy, we will first find the direcitonal derivative and then use following relation to get the gradient:

S = ∂W(U) ∂E(U) = ∂ µtr(E 2 ) + λ 2 (trE) 2 ∂E = µ ∂tr(E 2 ) ∂E + λ 2 ∂(
D w Φ = ∇φ, w n , (B.4)
where, w n is a unit vector in the direction of w.

B.2.1 Directional Derivatives of Some Useful Quantities

In the following, we consider a small variation of U by εw and compute the directional directive along w. Since we want to solve the equation for the displacement U and pressure p, we need to expand ∇ • P, which is the divergence of first Piola-Kirchoff stress tensor.

Let us denote displacement gradient by G, so from Equation A.8, we have: Equating quadratic terms in U to zero for linear approximation, we get,

F = I + G, (B.15) G =    ∂U 1 ∂X 1 ∂U 1 ∂X 2 ∂U 1 ∂X 3 ∂U 2 ∂X 1 ∂U 2 ∂X 2 ∂U 2 ∂X 3 ∂U 3 ∂X 1 ∂U 3 ∂X 2 ∂U 3 ∂X 3    , (
P = µ G + G T + λtr (G) I = [µ (G ij + G ji ) + λG kk δ ij ] e i ⊗

  's Disease (AD) . . . . . . . . . . . . . . . . . . . . 1 1.2 Characteristics of AD . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Biomarkers for AD . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3.1 Senile plaques . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3.2 Tau pathology . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.3.3 Neurodegeneration . . . . . . . . . . . . . . . . . . . . . . . . 7 1.4 Context and motivation . . . . . . . . . . . . . . . . . . . . . 10 1.5 Structure of the Manuscript and Contributions . . . . . . . 11 1.6 List of Publications and Awards . . . . . . . . . . . . . . . . . 13

Figure 1

 1 Figure 1.1: Diagram comparing normal adult brain (left) with a brain having AD (right) [Wikipedia 2016a]. The diseased brain has widespread diffuse cortical atrophy, extreme hippocampal atrophy, and ventricular enlargement.

Figure 1

 1 Figure 1.2: Figure reproduced with permission from [Bird 2008], which shows extracellular senile plaque (also called neuritic plaque) in the lower left corner and NFTs in the upper right corner. They are the most important characteristics of AD along with widespread brain atrophy and progressive cognitive impairment.

Figure 1 . 3 :

 13 Figure 1.3: Figure reproduced from [Herrup 2015], with permission, to illustrate the groups of three main categories of possible causes of AD: (top oval) cellular events; (bottom left oval) molecular events showing the sequences that result in two hallmarks of AD: NFTs and senile plaques; (bottom right oval) genetic events showing the risk factor genes where PSEN2, PSEN1 and APP are emphasized to indicate their involvement in familial AD (fAD). See [Herrup 2015] for detailed description.

  ]. Well established biomarkers for AD are the following: • Cerbrospinal fluid (CSF) based and imaging based biomarkers targeting Aβ pathology • CSF based and imaging based biomarkers (very recent and still under research) targeting taupathy • Flurodeoxyglucose based positron emission tomography (FGD-PET) targeting neurodegeneration • Structural MRI based morphometry targeting neurodegeneration and tissue atrophy • Cognitive impairment Cognitive impairment starts to appear many years after the changes that can be observed using the other biomarkers. The next sections describe briefly the biomarkers other than cognitive impairment.

Figure 1 . 5 :

 15 Figure 1.5: Figure reproduced from [Johnson 2012], with permission, to illustrate the variation of PiB PET images in normal control (NC), mild cognitive impaired (MCI), and AD subjects with a range of Aβ plaque deposition. Starting from the left: NC-are normal controls showing no evidence of Aβ deposition. NC+ are (around 25%) normal controls with some Aβ deposition. MCI-are around 40% -50% of patients with MCI showing no evidence of Aβ pathology, while the remaining show either moderate (MCI+) or severe Aβ deposition (MCI++). AD are the large majority of clinically diagnosed AD patients showing heavy Aβ deposition. The figure is based on several sources and studies that are detailed in [Johnson 2012].

Figure 1 . 6 :

 16 Figure 1.6: Figure reproduced from [Villemagne 2015], with permission, showing MRI and PET imaging of tau and Aβ in healthy controls, MCI and an AD patient. Top two rows show that Aβ deposition could be either present or absent in NCs. In second row tau deposits in NC is seen only in hippocampal region. MCI and AD patients have progressively increased depositions of both tau and Aβ. Although Aβ seems to be deposited almost everywhere in the brain regions, tau deposits are more concentrated in the cortical regions. MMSE=Mini Mental State Examination.

Figure 1 . 7 :

 17 Figure 1.7: Coronal slices of an AD patient aligned for comparison from two structural MRIs acquired two years apart. The ventricular expansions and the hippocampal atrophy are clearly observed on both right and left sides. Thinning of the cortical surfaces along with opening up of sulci can also be seen in the temporal lobes on either sides near the red horizontal line.

. 5 )

 5 Figure2.2 shows a simulation example where a 20% hippocampal atrophy is prescribed in addition to a uniform global 2% atrophy. For this case, we used ROBEX[Iglesias 2011] for skull stripping; FSL FAST[Zhang 2001] to segment the skull stripped image into GM, WM and CSF regions; and FSL FIRST[Patenaude 2011] to segment hippocampi and amygdalae to obtain the region of enhanced atrophy. This is a simple example to illustrate that any desired atrophy can be prescribed at voxel level once the desired segmentation of the input image is obtained. Solving the model using equation (2.5) with the prescribed atrophy map results in a deformation field as shown in the middle of Figure2.2. The deformation field is superimposed on the input baseline brain MRI. A simulated follow-up image is obtained by warping the baseline image using the computed deformation field. The difference between the real baseline and the simulated follow-up image is shown on the right of Figure2.2.

Figure 2 . 2 :

 22 Figure 2.2: Atrophy in cortex, white matter and Hippocampus. No atrophy is prescribed in brain stem and other sub-cortical structures such as Thalami, Putamen etc. Left: Prescribed atrophy map with uniform atrophy of 4%, 5% and 20% in WM, GM and Hippocampus respectively. Middle: A coronal slice of the input MRI superimposed with the obtained displacement field. Right: Difference between the input (baseline) and the simulated (follow-up) image.

Figure 2 . 3 :

 23 Figure 2.3: Prescribing uniform atrophy in 34 different cortical regions and in hippocampi. The regions are obtained from FreeSurfer whole brain segmentation. The atrophy values shown in the top row are the mean atrophy reported for amnestic MCI in[Carmichael 2013]. The computed atrophy map is the negative of the divergence obtained from the solution of the model when solved for the prescribed atrophy. This is shown in the bottom row. The negative values of the computed atrophy map (in blue) corresponds to the expansion of CSF to compensate the volume loss prescribed in the brain tissues. We can see here that the resulting atrophy of the parenchyma is exactly the one prescribed at the voxel level while the expansion of the CSF is not uniform.

Figure 2 . 4 :

 24 Figure 2.4: The figure shows several slices of a pair of real baseline and simulatedfollow-up MRIs with large uniform atrophy of 90% in cortical gray matter and Hippocampus, and of 40% in white matter. These large atrophies were simulated in a single time-step in this case. We can see that the model is able to simulate realistic pattern of widening and opening of sulci, the narrowing of gyri, and at the same time the cortical surface does not move unrealistically farther away from the skull.

Figure 2 . 5 :

 25 Figure 2.5: From left to right: Figure shows the initial prescribed atrophy, simulated follow-up images and the difference between the follow-up and baseline images for time-step 1, 3 and 6 respectively.

Figure 2

 2 Figure 2.6: Setting larger values of λ scales the force term to very high values and results in unrealistically large deformation fields.

Figure 2

 2 Figure2.8: We see that the simulated images are different when the values of µ is changed in the brain stem to be 100 times more (difference in µ overlaid over the baseline image in second column) even when there is no atrophy prescribed in this region. The difference is more pronounced in the brain stem and nearby regions.

Figure 2

 2 Figure 2.9: Impact of varying both µ and k. The second row shows that, just like in the previous figure, making the compressibility of the CSF k much smaller compared to µ makes CSF expand uniformly everywhere. The last row shows that when k is large, the CSF can expand very locally in response to the nearby tissue volume loss.

Figure 2

 2 Figure 2.10: For a given baseline image (first column), we prescribe two different atrophy patterns (second and fourth columns). The simulated images in these cases (third and fifth columns) look very similar. The axes in the images show a reference position to aid comparison of the images. The origins of the shown axes are at the same physical position for all the images in the same row.

Figure 2 .

 2 Figure 2.11: Registration results when using the input image of Figure 2.10 as a fixed image (first column) and the two simulated images of that figure as moving images (not shown here). Second and fourth columns show the results of the alignment while the third and fifth columns show the Jacobian determinants of the transformation obtained in these two cases. We see that these Jacobian determinants inferring the volume changes in the image are fairly similar. However these two moving images were created from the fixed image by prescribing very different underlying volume changes as shown in Figure 2.10.

Figure 2 .

 2 Figure 2.12: Brain edge movement reported by gBSI [Prados 2015] for a real baseline and two different simulated follow-up images: Atrophy prescribed exclusively in Gray Matter (on left) vs atrophy prescribed exclusively in White Matter (on right). In both cases the global atrophy prescribed is 4%. The underlying volume changes prescribed for these two cases are shown in the first and third columns of Figure 2.13. The brain edge movement and the reported atrophy are strikingly similar based on visual inspection.

2. 5 .

 5 Investigating the Relationship Between Image Appearance and Atrophy Patterns 36logDemons[START_REF] Lorenzi | [END_REF]], to estimate the local atrophy pattern for the two different scenarios. In this case we notice that the resulting atrophy patterns are different, and mostly localized in white and gray matter respectively. However, the estimated atrophy patterns are still different from the prescribed ones. This is expected as the registration algorithm is unaware of the underlying model used in simulating the images. Furthermore, we also see that changing the regularization schemes changes the results of the volume changes. The two regularization schemes used were the penalization of harmonic energy and the Gaussian smoothing respectively. The harmonic energy penalization (Reg1 in the figure) results in a sharper Log-Jacobian maps while the Gaussian regularization (Reg2 in the figure) results in the more diffused atrophy. These parameters of the implicit model used in the registration algorithm are difficult to relate to the underlying model that generated the time series images.

Figure 2 .

 2 Figure 2.14: For the same prescribed uniform volume loss only in cortical gray matter, registration can provide different values of the measured volume changes in white matter. In the above figure we can see that in the white matter region on the right of the intersecting lines, the registration estimates volume loss differently when changing the parameters of the model without changing the underlying volume changes.

  Complex Patterns of Patient-specific Atrophy for the Evaluation of Atrophy Measurement AlgorithmsIn section 2.5.2 we presented a qualitative analysis of the relationship between the actual underlying deformation and the volume changes inferred by atrophy estimation algorithms. In order to provide a quantitative assessment of the simulation results, in this section we provide a framework that allows simulating patient-specific atrophy patterns in large number of patients. The framework could be used as a starting point for either calibrating the model parameters or as a framework for benchmarking atrophy estimation algorithms. The general pipeline of the proposed framework is shown in Figure2.15.

Figure 2 .

 2 Figure 2.15: Pipeline illustrating the measurement of atrophy from i) RBRF: real baseline with respect to real follow-up, and ii) RBSF: real baseline with respect to simulated follow-up. The two atrophy measurement tools (AMTs) shown above as AMT1 and AMT2 can either be the same tool or different tools depending on the goals of the experiment.

  Figure2.17: FreeSurfer volume measurement comparison on real follow-up vs simulated follow-up image pairs using Absolute Symmetric Percentage Volume Change (ASPVC). The top end of the bars are the mean ASPVC values while the ticks show the plus/minus standard deviation from the mean. Left: Best twenty structures for which the mean ASPVC was found to be the lowest. Right: Only the brain structures for which the mean ASPVC of FreeSurfer cross-sectional run and longitudinal run were reported in[Reuter 2012]. In the rightmost two bars we present the average and standard deviation of mean ASPVCs over all the available regions in Figure7of[Reuter 2012], which used a different dataset TT-115 consisting of testretest same day repeat scan pairs of 115 healthy controls. Although not directly comparable due to the use of different datasets, we do see that the mean ASPVC computed from our real vs simulated follow-up image pairs are of the same order as that of the one presented in[Reuter 2012].

Figure 3 .

 3 Figure 3.2: Figure shows the experimental setup that is used in this paper to simulate mid-point image from the two available images. We use FreeSurfer as the atrophy estimator shown in Figure 3.1, and the atrophy predictor in Figure 3.1 is taken to be a basic linear interpolator.

Figure 3 . 3 :

 33 Figure 3.3: Boxplot of the atrophy estimates for the real mid-point images in the coritcal regions and hippocampus. These are the regions that were used in [Carmichael 2013]. The data shows the distribution of FreeSurfer atrophy estimates in the AD population of the MIRIAD dataset when considering the first and the mid-point scans.

Figure 3 . 4 :

 34 Figure 3.4: Boxplot of the difference in the FreeSurfer atrophy estimate in the real mid-point image and the interpolated mid-point image for all the 46 AD subjects present in MIRIAD dataset. The regions shown are the same as the one shown in Figure 3.3 and are displayed in the same order. The interpolated mid-point image is obtained by simulation using the pair of extremal time-point images of each subject.

Figure 3 . 5 :

 35 Figure 3.5: Absolute symmetric percentage volume change (ASPVC) error computed for the best 20 regions (on left), and Reuter regions (on right). Reuter regions are the regions for which ASPVC results were shown in[Reuter 2012] by computing FreeSurfer volume measurements on same day test-retest scan pairs of 115 healthy controls. The bar plots in orange, cross-run-mean and long-run-mean, are mean error of the Reuter regions shown on the right. cross-run-mean were computed using FreeSurfer cross sectional stream, while long-run-mean were computed using FreeSurfer longitudinal stream. ASPVC errors for the real vs. simulated mid-point images are in the same order as that of the mean error for real baseline repeat scans. It is worth noting that the dataset used in this work is not the same as the one used in[Reuter 2012], which is not publicly available.

Figure 3 .

 3 Figure 3.6: For each region, the box plots on the left (green) show estimated atrophy from the real mid-point images using FreeSurfer (a m in Figure 3.2), while the one on the right (orange) are from the FreeSurfer atrophy estimates of simulated mid-point images (a s m in Figure 3.2). The brain structures shown are the regions on Reuter-regions shown on the right of Figure 3.5. We see reduced variablity in atrophy estimation of the population when using simulated follow-up images.

Figure 3 .

 3 Figure 3.7: Same as Figure 3.6, but the results are for the brain structures shown on the left of Figure 3.5.

  Part of this work was funded by the European Research Council through the ERC Advanced Grant MedYMA 2011-291080.

Figure 4 . 2 :

 42 Figure 4.2: Examples of two different kinds of atrophy maps. The first row prescribes atrophy map that is uniform in different regions of the brain, while the second row prescribes smoothly varying atrophy. Both of these atrophy maps have same average values in each ROIs. The example also shows that we can prescribe volume changes in ventricles, if desired, by adapting the input segmentation map accordingly. The simulated images, as shown, are different although they have same mean regional atrophy values.

Figure 4 .

 4 Figure 4.2 shows two simulated images from the same input baseline image but with two different atrophy patterns.

Figure 4 . 3 :

 43 Figure 4.3: I b 0 and I b 1 are the repeat scans of a subject taken within a short period of time during which there is no morphological changes in the brain of the subject. I r is taken at a later time when the brain could have undergone some morphological changes. The deformation field Φ reg is obtained by registering I r to I b 0 , while Φ sim is obtained from the brain deformation model using I b 0 as the input image. The three simulated images I s 0 , I s 1 and I s 1 are all same time-point images but have different intensities that come from I b 0 , I b 1 and I r respectively.

Figure 4 . 4 :

 44 Figure 4.4: A general approach to simulate ground truth synthetic longitudinal images with realistic intensity variations; simulated images are shown within the shaded regions. The deformation fields with a prescribed atrophy for three timepoints (Φ sim 1 , Φ sim 2 , and Φ sim 3 ) are obtained from the biophysical model using I b 0 as the input baseline image. Several different sets of longitudinal images can then be simulated by resampling intensities from different combinations of available real images. The topmost shaded region shows a longitudinal sequence with no realistic intensity variations where the synthetic images are all resampled from I b 0 . The remaining two shaded regions have longitudinal sequences with realistic intensity variations where the simulated images are resampled from other available images of the same subject. In the ideal case, the three sets of longitudinal sequences have exactly the same morphological changes but with different variations in intensity characteristics.

Figure 4 . 5 :

 45 Figure 4.5: Two simulated images are shown on the third row where the image on the left is resampled from the input baseline image I b , and the image on the right is resampled from another image I r of the same subject. Both I b and I r had already been corrected for the bias field intensity inhomogeneity. The intensity histograms shown are of a selected ROI (shown on the last row) where there is no significant morphological changes between the images. From the histograms we can see that the simulated image I s 2 ,t 1 has a different intensity characteristics than I b , while the simulated image I s 1 ,t 1 has intensity characteristics that closely matches to that of I b .

Figure 4 .

 4 Figure 4.6: Two sets of synthetic longitudinal images are shown which are simulated by prescribing a smoothly varying atrophy pattern. The first row shows the input prescribed atrophy and the input baseline image I b of a subject, while the remaining rows show the two sequences. The sequence shown on the left have simulated images that are all resampled from I b . On the right, each simulated image is resampled from real MRIs of the same subject but taken at different times. As shown by the intensity histograms of Figure 4.7, the longitudinal synthetic images on the right have more realistic intensity variations than the one left.

Figure 4 .

 4 Figure 4.7: Intensity histograms of selected patches of the images simulated in Figure 4.6. When the simulated images are resampled from the same input baseline image I b , as expected, the histograms of the simulated images closely match with each other. However, when simulated images are resampled from other different images of the same patients, the histograms of these simulated images do not match closely. The longitudinal sequence of simulated images I s 2 ,t 1 , I s 2 ,t 2 and I s 2 ,t 3 has realistic variation in intensities as observed in the real sequences.

4. 4

 4 Simul@trophy: choices available and practical considerations Simul@trophy is available as an open-source repository under git version control. Researchers can use it according to their needs, improve the presented model, and/or add new models of brain atrophy. It is based on two core components: i) The Insight ToolKit (ITK) and ii) PETSc[Balay 2013]. All the input and output images of the brain deformation model shown in Figure4.1 can be in any format that ITK supports. ITK has strongly promoted reproducible science in the medical imaging domain, and has been widely used in computational science applied to medical imaging[McCormick 2014, Avants 2015]. Similarly, implementation of the model solver is based on open-source PETSc, a library based on C programming language.

Figure 4 . 8 :

 48 Figure 4.8: The figure shows an example of simulating a longitudinal sequence with backward time-points. The input baseline image I b is the same one as used in Figure 4.6, and the prescribed atrophy map is the negative of the map used in Figure 4.6. In the figure, we can see the shrinkage of the ventricles and the growth of the brain parenchyma.

Figure 4 . 9 :

 49 Figure 4.9: The figure shows an example of simulating follow-up images of a normal subject with baseline image I b , where the prescribed atrophy pattern is adapted from an AD patient. The prescribed atrophy is adapted from the atrophy estimated for the AD patient shown in Figure 4.6. Average values of the smoothly varying prescribed atrophy shown in Figure4.6 is computed in all the ROIs. The ROIs are obtained from the FreeSurfer segmentation including all the white matter parcellations[Fischl 2002]. The simulated images on the right have bigger shrinkage of the brain parenchyma and bigger expansion of the ventricles than the real images on the left.

•Figure 4 .

 4 Figure 4.10: RB and RF are non-rigidly registered and the transformation obtained from the registration is used to align RB to RF which is shown in the image RB_to_RF.The figure also shows two simulated follow-up images SF_in_RB and SF_in_RF that are resampled from (RB) and (RF) respectively. We can see that in most regions of the brain, the two simulated images have almost identical morphological appearances. However, there are also regions such as 2 and 5, where the morphological appearances of the two simulated images are not identical. From the registration results for these regions 2 and 5 in the zoomed patches, we can see that the registration is also not accurate in those regions.

Figure 4 .

 4 Figure 4.11: Standard staggered grid discretization scheme that is used to solve the system of Eqs. (4.1). Displacement variables are at faces (edges in 2D) of the cells, while pressure and atrophy values are at centres of the cells.

Figure 5 . 1 :

 51 Figure 5.1: High level block diagram for modeling and simulation of longitudinal MRIs in AD patients. Spatial and temporal distribution of neuronal deaths is represented in Atrophy Generation block which causes the brain shape changes represented in Brain Deformation block. This deformation along with the MRI acquisition conditions variability result intensity change in time series structural MRI of AD patients.The error in predicted follow-up from the actual observed follow-up MRI could also be used to optimize for the parameters of the developed models using a feedback system as shown above.

F

  we can express deformation gradient F in terms of displacement as: :D w E(U) = S : D w E = tr S T D w E = tr S T 1 2 D w (F T F) -D w I = tr S T 1 2 (∂w) T F + F T ∂w) = 1 2 tr S T (∂w) T F + S T F T ∂w) = tr(SF T ∂w) since, S T = S = tr(P∂w).Now, tr(P∂w) = tr P ij e i ⊗ e j ∂w k ∂X l e k ⊗ e l = tr P ij ∂w k ∂X l (e j • e k )e i ⊗ e l = tr P ij ∂w k ∂X l (δ jk )e i ⊗ e l = tr P ij ∂w j ∂X l e i ⊗ e l Pw)) = tr (∂(P ij e i ⊗ e j w k ek )) = tr (∂(P ij w k δ jk e i )) = tr (∂(P ij w j e i )) = tr ∂P ij w j ∂X k e i ⊗ e k = ∂P ij w j ∂X i = w j ∂P ij ∂X i + P ij ∂w j ∂X i = (∇ • P) T w + tr(P∂w). Using tr (∂(Pw)) = ∇ • (Pw) in above equation, tr(P∂w) = ∇ • (Pw) -(∇ • P) T w. Using B.8, D w W(U) = ∇ • (Pw) -(∇ • P) T w.(B.9)B.2.2 MinimizationOur functional to minimize is:R(U, p) = W(U) -p (∇ • U + ã). (B.10)So the Euler-Lagrange equations which are the necessary conditions for the minimizing function to satisfy are:∇ U R(U, p) = 0 (B.11) ∇ p R(U, p) = 0, (B.12)where, ∇ U R(U, p) and ∇ p R(U, p) are gradients of the functional R(U, p) with respect to the functions U and p respectively. Now,D w R(U, p) = Ω D w W(U) -Ω pD w (∇ • U + ã) = Ω ∇ • (Pw) -(∇ • P) T w --Ip) w) T n + Ω -(∇ • P) T w + (∇p) T w = -∇ • P + ∇p, w + ∂Ω ((P -Ip) w) T n,where n is a unit normal vector in the boundary ∂Ω, and we have used Divergence theorem to change some of the volume integrals into surface integrals. Assuming ∂Ω ((P -Ip) w) T n to be zero, we have, D w R(U, p) = -∇ • P + ∇p, w .Hence,∇ U R(U, p) = -∇ • P + ∇p.The constraint ∇ • U + ã is true for every small volume element in the domain Ω, hence the necessary conditions for minimum from Euler-Lagrange equations in B.11 are∇ • P -∇p = 0 ∇ • U + ã = 0. (B.14)

  2µE + λtr(E)I] F T = µ G + G T (I + G) T + λtr (G) I (I + G) T = µ G + G T + GG T + G T G T + λtr (G) I + λtr (G) G T .

  Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2 Biophysical Model of Brain Deformation . . . . . . . . . . . 19 2.2.1 Impact of Loss of Volume on Conservation Equation . . . . .

	1.6. List of Publications and Awards	13
	1.6 List of Publications and Awards	Chapter 2
	Journal Articles A Biophysical Model of Brain
	• [Khanal 2016b] Bishesh Khanal, Marco Lorenzi, Nicholas Ayache and Xavier Pennec. A biophysical model of brain deformation to simulate and analyze Deformation to Simulate and
	longitudinal MRIs of patients with Alzheimer's disease. NeuroImage, 2016. Chapter 2. Analyse Longitudinal MRIs of
	• [Khanal 2016a] Simulating Realistic Synthetic Longitudinal Images with known Volume Changes. Under review in NeuroImage. Chapter 4. Patients with Alzheimer's Disease
	Selective Peer-Reviewed Conference Papers
	Contents
	2.1

•

[START_REF] Khanal | A Biophysical Model of Shape Changes due to Atrophy in the Brain with Alzheimer's Disease[END_REF]

] Bishesh Khanal, Marco Lorenzi, Nicholas Ayache and Xavier Pennec. A Biophysical Model of Shape Changes due to Atrophy in the Brain with Alzheimer's Disease. MICCAI, Boston, USA, October 2014. This work was extended to the journal article

[Khanal 2016b

], Chapter 2 with some significant additions.

•

[Khanal 2016c

] Bishesh Khanal, Marco Lorenzi, Nicholas Ayache and Xavier Pennec. Simulating Patient Specific Multiple Time-point MRIs From a Biophysical Model of Brain Deformation in Alzheimer's Disease. Computational Biomechanics for Medicine X, MICCAI Workshop, Munich 2015. Chapter 3. Awards • 1st prize in "Popular Vote Awards" at the MICCAI 2014 Educational Challenge 2 , Boston, USA, October 2014 for the video: "Statistics on Lie groups for Computational Anatomy", Nina Miolane and Bishesh Khanal (out of 16 international submissions ). • "Student Travel Award" at MICCAI conference 2014 in Boston, USA for the paper: "A Biophysical Model of Shape Changes due to Atrophy in the Brain with Alzheimer's Disease." [Khanal 2014]. 20 2.2.2 Constrained Minimization of the Elastic Energy . . . . . . . 20 2.2.3 Modeling CSF Region . . . . . . . . . . . . . . . . . . . . . . 21 2.3 Implementation of the Biomechanical Tissue Loss Model . 23 2.3.1 Skull Stripping and Brain Segmentation . . . . . . . . . . . . 23 2.3.2 Input Prescribed Atrophy Map . . . . . . . . . . . . . . . . . 23 2.3.3 Staggered Grid Discretization and Finite Difference Method . 23 2.4 Using the Model to Explore Different Atrophy Patterns . . 24 2.4.1 Simulating Images with Any Desired Atrophy Maps . . . . . 24 2.4.2 Simulating Realistic Atrophy Patterns . . . . . . . . . . . . . 25 2.4.3 Simulating Large Atrophy With Multiple Time-steps . . . . . 26 2.4.4 Role of Different Model Parameters . . . . . . . . . . . . . . 27 2.5 Investigating the Relationship Between Image Appearance and Atrophy Patterns . . . . . . . . . . . . . . . . . . . . . . . 32 2.5.1 A Synthetic Example with Binary Image . . . . . . . . . . . . 32 2.5.2 Varying the Spatial Distribution of Atrophy in Real Brain Images 34 2.6 Simulating Complex Patterns of Patient-specific Atrophy for the Evaluation of Atrophy Measurement Algorithms . . 37 2.7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 2.8 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

  biologically plausible deformation for any given desired volume changes at the voxel level of the brain MRI. Large number of brain regions are segmented in 46 AD patients and the atrophy rates per year are estimated in these regions from the available two extremal time-point scans. Assuming linear progression of atrophy, the volume changes in scans closest to the half way time period is computed. These atrophy maps are prescribed to the baseline images to simulate the middle time-point images. The volume changes in real middle time-point scans are compared to the ones in simulated middle time-point images. This present framework also allows to introduce desired atrophy patterns at different time-points to simulate non-linear progression of atrophy. This opens a way to use a biophysical model of brain deformation to evaluate methods that study the temporal progression and spatial relationships of atrophy of different regions in the brain with AD.

	Chapter 3
	Simulating Patient Specific
	Multiple Time-point MRIs From a
	Biophysical Model of Brain
	Deformation in Alzheimer's
	Disease

5. We would like to thank neurologist David Bensamoun of the Nice Resource & Research Memory Centre (Nice, France) for providing qualitative assessment of our simulations of atrophy progression in the brain MRIs. Contents 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.2 Biophysical Model of Brain Deformation due to Atrophy . 49 3.2.1 Modeling CSF Region . . . . . . . . . . . . . . . . . . . . . . 50 3.2.2 Staggered Grid Discretization and Finite Difference Method . 51 3.3 Simulating New Time-points from Two Available scans of MRIs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.4 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . 53 3.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 Abstract This paper proposes a framework to simulate patient specific structural MRIs from the available scans of Alzheimer's Disease(AD) subjects. We use a biophysical model of brain deformation due to atrophy that can generate

  Contents 5.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.2 Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.2.1 Creating Databases of Synthetic Longitudinal MRIs for Machine Learning Applications . . . . . . . . . . . . . . . . . . 86 5.2.2 Optimisation of Model Parameters . . . . . . . . . . . . . . . 87 5.2.3 Towards an Integrative Multimodal Model . . . . . . . . . . . 87

  e j . Patrizia Formichi, Carla Battisti, Elena Radi and Antonio Federico. Cerebrospinal fluid tau Aß and phosphorylated tau protein for the diagnosis of Alzheimer's disease. J.Cell. Physiol., vol. 208, no. 1, pages 39-46, jan 2006. (Cited on page 7.) [Freeborough 1997] P.A. Freeborough and N.C. Fox. The boundary shift integral: an accurate and robust measure of cerebral volume changes from registered repeat MRI. IEEE Transactions on Medical Imaging, vol. 16, no. 5, pages 623-629, 1997. (Cited on pages 34 and 60.) [Frisoni 2010] G. B. Frisoni, N. C. Fox, C. R. Jack, P. Scheltens and P. M. Thompson. The clinical use of structural MRI in Alzheimer disease. Nature reviews. Neurology, vol. 6, no. 2, pages 67-77, feb 2010. (Cited on pages 9, 10, 15, 20, 48, 50 and 60.) [Gholami 2016] Amir Gholami, Andreas Mang and George Biros. An inverse problem formulation for parameter estimation of a reaction-diffusion model of low grade gliomas. Journal of Mathematical Biology, vol. 72, no. 1, pages 409-433, 2016. (Cited on page 87.) [Goedert 1993] Michel Goedert. Tau protein and the neurofibrillary pathology of Alzheimer's disease. Trends in Neurosciences, vol. 16, no. 11, pages 460-465, nov 1993. (Cited on page 7.) [Gorgolewski 2015] Krzysztof Jacek Gorgolewski, Gael Varoquaux, Gabriel Rivera, Yannick Schwartz, Satrajit S Ghosh, Camille Maumet, Vanessa V Sochat, Thomas E. Nichols, Russell A. Poldrack, Jean-Baptiste Poline, Tal Yarkoni and Daniel S. Margulies. NeuroVault.org: A web-based repository for collecting and sharing unthresholded statistical maps of the human brain. Frontiers in Neuroinformatics, vol. 9, no. 8, 2015. (Cited on page 68.) [Gudbjartsson 1995] Hákon Gudbjartsson and Samuel Patz. The Rician Distribution of Noisy MRI Data. Magnetic resonance in medicine : official journal of the Society of Magnetic Resonance in Medicine / Society of Magnetic Resonance in Medicine, vol. 34, no. 6, pages 910-914, December 1995. (Cited on page 61.) [Hadj-Hamou 2016] Mehdi Hadj-Hamou, Marco Lorenzi, Nicholas Ayache and Xavier Pennec. Longitudinal Analysis of Image Time Series with Diffeomorphic Deformations: a Computational Framework based on Stationary Velocity Fields. Frontiers in Neuroscience, vol. 10, no. 236, 2016. (Cited on pages 68 and 82.) [Hardy 1992] J. Hardy and G. Higgins. Alzheimer's disease: the amyloid cascade hypothesis. Science, vol. 256, no. 5054, pages 184-185, apr 1992. (Cited on pages 3 and 4.)

	[Formichi 2006]
											(B.18)
	Now we expand the divergence as,
	∇ • P =	∂P ij ∂X j	e i
	= µ	∂G ij ∂X j	+ µ	∂G ji ∂X j	+ λ	∂G kk δ ij ∂X j	e i
	= µ	∂ ∂X j	∂U i ∂X j	+ µ	∂ ∂X j	∂U j ∂X i	+ λ	∂ ∂X i	∂U j ∂X j	e i
	= µ	∂ 2 U i ∂X 2 j	+ (µ + λ)	∂ ∂X i	∂U j ∂X j	e i
	= µ∆U + (µ + λ) ∇ (∇ • U) .	(B.19)

Oligomer is a macromolecular complex consisting of a few macromolecules like proteins or nucleic acid[Wikipedia

2015].

Details at: http://www.miccai.org/edu/videos.html#mec2014winners

The simulation results are made available at http://neurovault.org/collections/ AUKWWYBC/ [Gorgolewski

2015].

Acknowledgements

I would like to first thank my advisors Xavier Pennec and Nicholas Ayache. Thanks Xavier for mentoring me from the very beginning of this work, and playing a key role in laying out the foundation of the work. I really appreciate your dedication to research and caring of your students. Thank you Nicholas for all helpful advice and to let me have two mentors with complementary qualities. Being in Asclepios under your leadership, I have experienced a learning example on how to build a strong team with a friendly environment. Let me also thank Hervé Delingette and Maxime Sermesant for the valuable discussions and completing up an awesome team of research scientists in Asclepios. Acknowledgements 1. Part of this work was funded by the European Research Council through the ERC Advanced Grant MedYMA 2011-291080. Acknowledgments 1. We would like to thank Mehdi Hadj-Hamou for providing us registration results and the associated deformation fields that were used in this paper to resample intensity from different images. The preprocessing steps involved for this registration are explained in [Hadj-Hamou 2016].

2. Data used in the preparation of this article were obtained from the MIRIAD database. The MIRIAD investigators did not participate in analysis or writing of this report. The MIRIAD dataset is made available through the support of the UK Alzheimer's Society (Grant RF116). The original data collection was funded through an unrestricted educational grant from GlaxoSmithKline (Grant 6GKC).

3. This work benefited from the use of the Insight Segmentation and Registration Toolkit (ITK), an open source software developed as an initiative of the U.S.

Appendix

Appendix A

Deformation theory

A.1 Body, Configurations and Motion

Here we describe briefly mathematical form that is widely used to characterize deformation and motion of materials. The materials presented in this appendix have been adapted from the sources in [Naghdi 1994] and [Kelly 2012] unless cited otherwise.

A body is an abstract mathematical entity that models physical material. It consists of continuous form of matter, where small portions of this matter are known as material particles such that every material particle in the body can be put into one to one correspondence with an Euclidean space E 3 .

Let us consider a body B. Now we can define a configuration which is a region R in E 3 such that each point X c in R is mapped one to one from material particles X of the body B. Thus a configuration may also be considered as a mapping:

Motion can be described by a set of configurations obtained depending on time as:

x = Φ( X, t).

(A.1)

Above equation describes a motion of a body in terms of material particles and hence is known as a material description of motion.

It is common to choose a reference configuration so that the motion can be described relative to this configuration. In this work we will use X as a reference position corresponding to a material particle X of the body. Its corresponding positions in other configurations will be denoted by x. Referential description of motion, also known as Lagrangian description is given by: x = Φ(X, t).

(A.2)

A configuration that body occupies at time t is known as current configuration. Spatial description focuses on the positions in current configuration and finds corresponding positions in reference configuration. This is also known as Eulerian description and is expressed as: Velocity is the time derivative of x, that is:

A.2 Deformation and Strain

A deformation gradient F, is defined as

F describes a local deformation of a material particle whose position is at X in referential configuration. It is also a second order tensor that transforms a line element dX in reference configuration to a line element dx in current configuration, that is:

Displacement gradient in referential co-ordinates is given by,

where I is a second order identity tensor.

Left Cauchy-Green strain tensor gives us a measure of how the lengths of line elements and angle between them change between reference and current configurations and is defined as:

(A.9)

Green-Lagrange strain tensor is a relative measure of strain which vanishes for rigid motions and is given by:

Using A.8, we can express this strain in terms of displacement gradients as:

Jacobian determinant or sometimes simply referred to as Jacobian is the determinant of the Jacobian of transformation from X to x, that is:

J is a measure for the change of volume of a material element due to the deformation, and is related to the volume change by, J = dv dV .

(A.13)

A.3. Saint Venant-Kirchoff Model of Hyperelastic Materials

92

A.3 Saint Venant-Kirchoff Model of Hyperelastic Materials

In this model, stress is related to the Green-Lagrange strain tensor as [Wikipedia 2016b]:

where, λ and µ are Lamé's parameters. S is the second Piola-Kirchoff stress tensor. First Piola-Kirchoff stress tensor is related with S as:

Strain-energy density function is given by:

S can also be obtained from the following relation:

Most of the materials in this part of the Appendix are adapted from [Kelly 2012]. For any two vectors u and v, a dyadic product results in a second order tensor and u ⊗ v is known as a dyad. A dyad can be defined such that it transforms any vector w in the following way:

w).

(A.18)

A second order tensor A with components A ij can be expressed in cartesian co-ordinates using a dyadic product and an index notation as:

(A.19)

In the index notation, if the index is present exactly twice in the same term, it implies summation over that index.

Simple contraction of a tensor and vector is given by:

Double contraction of dyads are defined as:

Transpose of a second order tensor can be expressed as: Let us denote directional derivative of W(u) in a direction w by D w W(u) which is given by,

Other useful results:

Derivation of the System of Equations of the Biophysical Model of Brain Deformation

B.1 Derivation of the Conservation Law

The atrophy rate ã(x, t) at any position x at time t for a representative elementary volume of V xt is defined as the negative rate of change of volume per unit volume:

Let us consider a sufficiently small deformation induced in a time interval ∆t. Let φ denote the deformation of the material during this time. The new position of a material particle initially at reference position X is given by: x = φ(X) = X + u = X + ∆tv where u is the displacement of the particle at position X and v is the particle's velocity.

Let V t and V t+∆t denote the elementary volume of a material at time t and t+∆t respectively. By the definition of atrophy rate, we have,

Now, a = ã∆t which is the amount of atrophy during time ∆t is given by,

where J is the Jacobian determinant given by, J = det (∇φ) = det (∇ (X + ∆tv)) = det (I + ∆t∇v)