Structural Magnetic Resonance Imaging (MRI) has been widely used for in vivo observation of morphological changes over time in human brain. Atrophy or tissue volume loss measure from structural MRI is an established biomarker for neurodegeneration (Frisoni et al., 2010). There is a large number of brain morphometry algorithms developed in the literature which estimate global or local atrophy from structural MRIs (Wright et al., 1995; Freeborough and Fox, 1997; Ashburner and Friston, 2000; Smith et al., 2002; Hua et al., 2008). Volume/atrophy measurements obtained from such algorithms have been used to test various clinical hypotheses about neurodegenerative diseases (Wright et al., 1995; Sepulcre et al., 2006; Koch et al., 2016). Similarly, comparison of different neurodegenerative diseases have also been performed based on these measurements (Rosen et al., 2002; Whitwell and Jack, 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 and Ridgway, 2015). Such methods contain model biases because TBM results depend on the choices of regularization used during the registration of images (Ashburner, 2013). Likewise, edge-based methods such as BSI, SIENA etc. are sensitive to unmatched image contrasts between scans, poor signal-to-noise ratio, partial volume effects, segmentation errors etc. (Preboske et al., 2006; Prados et al., 2015). 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 et al. (2013) studied the groupings of 34 cortical regions and hippocampi from the per-individual rates of atrophy estimates in these regions. In Fonteijn et al. (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 et al., 2003; Camara et al., 2006; Karaçali and Davatzikos, 2006; Pieperhoff et al., 2008; Sharma et al., 2010; Modat et al., 2014; Radua et al., 2014; Khanal et al., 2016a). Most of 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 (Camara et al., 2008; Pieperhoff et al., 2008; Sharma et al., 2010), to estimate the bias in such algorithms, and also to estimate uncertainty in the measured atrophy (Sharma et al., 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 (Camara et al., 2008; Pieperhoff et al., 2008; Sharma et al., 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 and Patz, 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.

For example, Sled et al. (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, Camara et al. (2008) added noise to the simulated ground truth images with atrophy.

In addition to the Rician noise described above, other noise, and artifacts are also present in MRIs (Simmons et al., 1994). Some of the artifact sources that have been shown to affect the measurements of atrophy estimation algorithms (Camara-Rey et al., 2006; Pieperhoff et al., 2008; Sharma et al., 2010) are:

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

Many other acquisition artifacts may not be simulated because we do not have faithful models. Inability to produce realistic intensity variation and noise in simulated longitudinal images is one of the key limitations in the state-of-the-art atrophy simulators, including our previous work (Khanal et al., 2016a). In this work, we propose a simple but elegant solution to remove the limitation of previous atrophy simulators. First, our biophysical model of brain deformation (Khanal et al., 2016a) is used to obtain a dense deformation field with specified volume changes. Then, to obtain realistic intensity variations, intensities in the simulated images are resampled from baseline repeat scans of the same patient. Although the method is very simple and straightforward, this allows simulating longitudinal images with variation in intensity and noise taken from real scans themselves without explicitly specifying any noise or artifact models. To the best of our knowledge, this idea was not presented before in the literature. When the repeat scans are not available, we use an approach introduced by Prakosa et al. (2013) 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 the repeat scans if available, otherwise from the real images of the same patient taken at different times.

Figure 1 shows a diagram of the complete framework. To implement this framework, we have developed an open-source atrophy simulator software called Simul@trophy¹. To our knowledge, Simul@trophy is the first atrophy simulator to be made open-source. Simul@trophy uses the biophysical model presented in Khanal et al. (2016a) but introduces a new numerical scheme to compute divergence, which removes the numerical inconsistency presented in the previous work. This is further explained in detail in Section 4.2.

Section 2 explains all the blocks of the framework shown in Figure 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 3 shows some simulation results using Simul@trophy, and also illustrates some potential applications of the simulator. In Section 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.

We use the biophysical model presented in Khanal et al. (2014, 2016a) 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 1, are: (i) Pre-processing, (ii) Brain deformation model, (iii) Realistic intensity simulator.

This section presents simulation examples of synthetic longitudinal MRIs with prescribed atrophy patterns and realistic intensity variations². The real input MRIs used for the simulations presented in this section come from the database made available by Hadj-Hamou et al. (2016). The images had already undergone the Pre-Processing and Position Correction steps of the Longitudinal Log-Demons Framework (LLDF) detailed in Hadj-Hamou et al. (2016). Starting from the publicly available OASIS dataset (Marcus et al., 2010), the images in the database had undergone intensity inhomogeneity correction using ANTs–N4BiasFieldCorrection (Avants et al., 2011), and had been transported to a common space using affine registration with FSL–FLIRT (Jenkinson and Smith, 2001).

Since all the simulated images must undergo interpolation of intensities, numerical scheme used in the interpolation will have an impact on the intensity characteristics of the simulated images. In all the simulation examples that follows, intensities were resampled using B-spline interpolation of order 3.

Figure 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:

Once all the deformation fields Φ_si corresponding to u_i for i = 0, 1, …, n are obtained, these deformation fields can be used as shown in Figure 4 to simulate different sequences of longitudinal images. As time step gets larger, the segmentation map is warped with an increasingly bigger displacement field using nearest neighbor interpolation, which could result in numerical instabilities. As the atrophy map is also warped at each time step, the global atrophy rate prescribed in the beginning is not necessarily preserved during the intermediate time-steps.

In Figure 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. In brain parenchyma regions, it is the negative of the divergence of a stationary velocity field obtained by performing LCC log-Demons registration (Lorenzi et al., 2013) 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 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 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 6. From the segmentation image shown in Figure 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 9, we show an example where synthetic sequence of images is simulated by starting from a baseline image of a healthy subject. However, the prescribed atrophy is derived from an atrophy estimated from the AD patient used in Figure 6. The input baseline images of both the AD patient and the healthy subject were segmented using FreeSurfer (Fischl et al., 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 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.

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 et al. (2013). All the input and output images of the brain deformation model shown in Figure 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 et al., 2014; Avants et al., 2015). Similarly, implementation of the model solver is based on open-source PETSc, a library based on C programming language. 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 choices are provided via command line arguments. The available command lines are detailed in Appendix in Supplementary Material. In this section, we illustrate some examples of how certain choices made during the simulation affect output results.

In Khanal et al. (2016a), 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 et al. (2016b), 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. This fills an important gap in the existing literature to simulate atrophy in longitudinal images with realistic intensity variation without explicitly modeling the noise and acquisition artifacts. 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. 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.

In Figure 6, smoothly varying atrophy pattern was prescribed by taking the negative of the divergence of a stationary velocity field obtained by registering the input baseline image with a follow-up image of the same subject. The objective of the experiment was to illustrate the ability of Simul@trophy to simulate smoothly varying patterns of atrophy in addition to the piecewise continuous atrophy maps. Registration was taken just as a means of getting a realistic smoothly varying atrophy maps; it is worth mentioning that simulating the deformation to be close to the deformation obtained from the registration algorithm was not the objective of this experiment. This is because the actual deformation field depends on the regularization used in the registration algorithm which does not necessarily follow the modeling assumptions used by Simul@trophy.

Although, the proposed method of resampling intensity from an image different from the input image provides more realistic variations, there are nevertheless certain issues one needs to be aware of. Since the simulated image has its intensities interpolated from another image, it can slightly reduce the noise variance. A neighborhood with expansion in the simulated image have intensities with slightly different linear combinations of intensities coming from a smaller set of voxels in the input image. Thus, the simulated image would have a smoother autocorrelation in the neighborhood compared to an equivalent real image. The fact that the simulated image has undergone interpolation and draws intensities from a limited set of raw voxels means that it is inherently smoother than the real scans. Finally, the usual spatial patterns of artifacts on real scans might not be exactly reproduced when warping real images. Any application using the simulated sets of images with the proposed approach should be aware of and ideally take into account these issues when interpreting results.

Use of repeat baseline scans to obtain intensity variation in the simulated images provides a very simple approach without using explicit noise and artifact models. One limitation with this is that the repeat baseline scans are not always available. When repeat scans are not available, we have proposed to use images at other time-points of the same subject, which requires performing non-linear registration. However, none of the non-linear registration methods are perfect and therefore the inaccuracies in registration affect the simulation results. This issue was discussed with illustrative examples in Section 4.1.

Simul@trophy can be used in evaluating atrophy estimation algorithms in similar ways as done by Pieperhoff et al. (2008), Camara et al. (2008), and Sharma et al. (2010). Since the proposed approach to simulate images may need deformations estimated from image registration, the use of these simulated images for the evaluation of some registration algorithms can bring an issue of circularity. This limitation adds to another limitation present in all publications related to atrophy simulation that we are aware of: namely, the models used in simulating images could favor certain kinds of registration algorithms over others. Although the ground truth atrophy can be measured from the combined deformation fields, the users must be aware of both limitations when they use Simul@trophy for the evaluation of registration algorithms.

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 et al. (2012), Fonteijn et al. (2012), Jedynak et al. (2012), Dukart et al. (2013), and Schmidt-Richberg et al. (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 subject-specific 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.

Simul@trophy could also be used in studying the role of morphology and intensity on atrophy estimation algorithms, and in machine learning based AD classification algorithms. Simul@trophy enables to perform such studies as it allows creating a large number of images by simulating atrophy patterns commonly observed in AD patients but with intensities taken from normal subjects and vice versa.

We hope to promote two directions of research in the community with open-source 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 et al. (2016a), 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.

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 core part of Simul@trophy is the implementation of our brain deformation model presented in Khanal et al. (2016a). 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.

BK has worked on the design and implementation of Simul@trophy. NA and XP have supervised him in the design of the experiments and in the methods presented in the paper. BK has written the manuscript with multiple iterations of suggestions and corrections from NA and XP.

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