Accurate automatic estimation of total intracranial volume: A nuisance variable with less nuisance

Data collection

We analysed T1-weighted MPRAGE scans of 288 (m/f 130/158) subjects aged 50–85 collected as part of a “real-world” multi-centre clinical trial (Fox et al., 2005, Gilman et al., 2005). Subjects met the NINCDS-ADRDA criteria for probable Alzheimer's disease (McKhann et al., 1984), mean age was 71.8(7.9) years, MMSE (Folstein et al., 1975) at baseline 20.4(3.3). Scans were coronal volumetric acquisitions lasting ≤ 7.5 min, slice thickness 1.5–1.8 mm adjusted to cover the entire brain, slice FOV 25 × 25 cm, and effective matrix size of 256 × 256 × 124. Acquisition parameters varied over the 17 MRI centres, full details can be found in Fox et al. (2005).

Manual estimation of TIV

Baseline scans were manually segmented for TIV by four validated operators according to the protocol described in Whitwell et al. (2001), which we summarise here for convenience. The intracranial volume is defined as the “volume within the cranium, including the brain, meninges, and CSF.” Measurements were conducted with the MIDAS software (Freeborough et al., 1997). Whole-brain volumes were first manually delimited using a 3D morphological method (Freeborough et al., 1997).

The T1-weighted volumes are rigidly registered to the Montreal Neurological Institute MNI305 brain average (Evans et al., 2012, Evans et al., 1993). A threshold of 30% of the mean brain signal intensity was used to outline the outer border of the dura as an aid to manual delineation of the outer edge of the intra-cranial volume, and the inferior limit of segmentation is set as the lowest slice in which cerebellar tissue is present. Every 10th axial section was segmented starting from the inferior limit working to the most superior slice with any brain tissue present. Slice areas are linearly interpolated to estimate the TIV for the intervening slices. Intra-rater and inter-rater variabilities were reported to show coefficients of variation (CV) 0.16% (n = 10) and 0.62% (n = 5).

Automatic estimation of TIV using FreeSurfer

FreeSurfer determines estimated TIV (known there as eTIV or just ICV) using an atlas scaling factor (i.e. the determinant of an affine transformation matrix) derived from registering images to an average template using a full (12-parameter) affine transformation (Buckner et al., 2004; see also http://surfer.nmr.mgh.harvard.edu/fswiki/eTIV). Segmentation is not used. Here, we use FreeSurfer 5.3.0 (the latest stable release as of April 2014), running “recon-all -autorecon1” and obtaining the ICV using “mri_segstats --etiv-only”.

Automatic estimation of TIV using SPM

There are several methods available to compute TIV using SPM's unified segmentation and spatial normalisation procedure. Methods can be broadly categorised into two main approaches:

• 1.

  The spatial normalisation transformation can be used, either inverse-transforming (preserving voxel values rather than volumes) a template-space TIV mask to the individual and determining the volume of the resultant individual-space mask (Keihaninejad et al., 2010) or (equivalently, apart from numerical errors) performing Jacobian integration (Boyes et al., 2006) over a template-space TIV mask.

• 2.

  Probabilistic tissue class images can be integrated (i.e. voxels are summed, accounting for the voxel volume) to give tissue volume estimates, with TIV simply being the sum of grey matter, white matter and CSF volumes. A subtlety here is that SPM can provide various sets of tissue class images: native, rigidly-reoriented (and resliced) to standard space, or non-linearly warped to standard space. With volume-preserving transformations¹ for the latter, all three sets of images should theoretically agree, except for the fact that they can have different fields of view; this can be important, since the amount of e.g. spinal cord contained in the three fields of view can differ, leading to different TIV estimates. The “modulated” non-linearly warped images (with mwc prefixes) should have the most consistent inferior cut-off, which may be the reason for their slightly better performance compared to the “native” subject-space segments in Ridgway et al. (2011).

It is important to note that the tissue prior probability templates used in SPM are based on averaging multiple automatically segmented images in standard space (for example, SPM12's priors come from segmentations (using New Segment) of images from the IXI data-set, http://www.brain-development.org/ (Heckemann et al., 2003)), so there is no guarantee that the sum of grey matter, white matter and CSF classes will be exactly consistent with accepted definitions of TIV, particularly with regard to the inferior cut-off and the inclusion of blood-filled sinuses. For this reason, we used the SPM12 tissue prior maps (and corresponding average T1-weighted, T2-weighted and proton-density weighted images from the same IXI data) to create a manually-corrected TIV mask consistent with the protocol described above (though segmented at each slice). Fig. 1 shows the TIV mask applied to tissue classes. Supplementary Fig. 1 shows a typical illustration of the non-brain classes, which are almost entirely located outside the ICV.

Although Ridgway et al. (2011) found only very small differences between SPM-based estimates related to approach 1 and approach 2 above, one theoretical advantage of the former is that it yields a contiguous TIV mask, less prone to isolated mis-segmentations far from the intracranial boundaries; a theoretical advantage of the latter approach is that the segmentation can potentially better model finer spatial detail (for example in the slightly more convoluted areas around the temporal lobe and cerebellum) than the regularised (smooth) spatial transformation. In an attempt to combine these advantages, we implemented a “Tissue Volumes” Utility in SPM12, which computes the totals of the modulated warped segmentations within the aforementioned manually-corrected TIV mask. This is available as a built-in SPM utility through the batch editor in the recent beta versions of SPM12.

The unified segmentation algorithm itself in SPM12 is similar to that in SPM8's New Segment, but with recomputed tissue priors (using multi-modal data from IXI, as mentioned above). An additional change is that the global rescaling of tissue priors present in SPM8's default segmentation but not New Segment was reintroduced.² In SPM12, each tissue probability map is rescaled by an additional (non-negative) parameter, and then re-normalised so that the priors sum to one at each point in space. The advantage of this more flexible model of SPM12 is that it allows for global decreases or increases in the amount of each tissue type. This is especially important for dealing with the kinds of atrophy seen in studies of ageing or dementia. The models for old and new segment and SPM12 are further detailed in Appendix A.

SPM12 TIVs were computed using the beta version of SPM12, revision 5647. For comparison, we also use SPM8 (revision 5236), simply summing modulated warped segments without a TIV mask. Finally the Supplementary Table 1 includes results for both the SPM12 Jacobian integration over the TIV mask and volume of the ICV mask transformed to subject-space.

Statistical analysis

Results were analysed in STATA 12. To assess suitability of automated TIV as a replacement for manual measurements we calculated squared correlation coefficients (R²) of automatic with manual measures. As the R² coefficient represents the degree to which variation in each variable is explained by the other, the R² between the two measures indicates the worst-case loss of explanatory power replacing one with the other as a correlate in a linear model. A high R² compared to the gold standard therefore indicates a method that can be used as a proxy for this purpose. Confidence intervals (CI, 95%) on R² coefficients and regression coefficients β were estimated using bootstrapping (20,000 samples each test, bias-corrected and accelerated). The same bootstrapping procedure was used to test the paired difference in R² coefficient between the automated TIV methods.

Bland–Altman (B–A) plots (Bland and Altman, 1986) were used to assess the agreement of values from the automated and manual TIV. It is expected that two measures of the same quantity should report the same result, that is: a slope of regression close to 1 (measurement error reduces the measured slope) and differences between measures due only to random error with mean 0 and a standard deviation that is acceptably small. Plotting difference against mean value allows assessment of bias and deviation from parity, using standard t-tests and linear regression. Pitman's test (significance of correlation of difference to mean) was applied to compare variance of the measures.

We do not attempt to compare ICV classification images as two of the methods do not produce them. Our reference, manual segmentation is performed only for every 10th axial section, and the FreeSurfer estimate uses only the atlas scaling factor.

Methods of adjusting for TIV

We have assumed above that adjustment for TIV will be performed by including it as a covariate in a linear model (also known as analysis of covariance, ANCOVA), rather than e.g. dividing regional brain volumes by TIV (known as the proportion method); it is also possible to use a previously fitted regression model (e.g. from a large normative study) to adjust volumes in individuals (known as the residual[isation] method). The relative merits of these methods have been debated in the literature. Arndt et al. (1991) demonstrated problems with the proportion method (reduced joint reliability), and noted that other methods such as regression, “may yield measurements that are more appropriate than ratios”. Mathalon et al. (1993) observed that, “whereas residual scores were uncorrelated with head size (by definition), measures taken as a proportion of head size tended to persist in showing correlations with head size,” which suggests that the former is preferable for adjustment; however, they also note that this property of the latter can be of interest in its own right, in terms of understanding the scaling laws of the brain (see also the discussion of allometry in O'Brien et al. (2006)). Sanfilipo et al. (2004), “found that the residual method generally was less affected by systematic and random errors in ICV and APV [absolute parenchymal volume] values, with the exception of dependent-related APV systematic error”, with the proportion method only being preferable in the latter case. Barnes et al. (2010) regressed logarithmically transformed regional volumes against log(TIV) and found confidence intervals often excluded unity, which further argues against the use of the proportion method; however, the question of whether to log-transform or not arguably remains open.

O'Brien et al. (2011) observe that the ANCOVA and residual approaches have the flexibility to be extended to model a quadratic effect of TIV, allowing for nonlinear relationships between regional volumes and head size. For the case of mass-univariate voxel-wise or vertex-wise analysis, the ANCOVA approach has the advantage that the model can straightforwardly vary over the brain (Barnes et al., 2010, Peelle et al., 2012). Furthermore, the ANCOVA model allows interaction terms to be modelled, e.g. between diagnostic group and TIV (O'Brien et al., 2011, Sanfilipo et al., 2004, Ueda et al., 2010). One could even consider higher order polynomial expansions of TIV – and/or logarithmically or otherwise transformed TIV – interacting with group or other variables. In combination with the above-surveyed advantages, the greater flexibility of the ANCOVA model leads us to recommend it as the first choice to consider. In situations where it is feasible (i.e. no more than a few ROIs) we would also recommend O'Brien et al. (2011) approach of graphically investigating the relationship between regional volume(s) and TIV.

On the use of “nonlinear-only modulation”

Two popular software packages for voxel-based morphometry, VBM8³ and FSL-VBM⁴ recommend a strategy to adjust for head-size by modifying the volume-preserving Jacobian-modulation step such that the affine component is ignored and only the nonlinear volume changes are preserved. This is equivalent to fully preserving the original volume with the usual (affine and nonlinear) Jacobian modulation and subsequently dividing by an estimated TIV obtained from the determinant of the affine transformation. We have not here evaluated the use of the affine determinant from SPM12's unified segmentation to estimate TIV, but there seems no reason to expect SPM12's affine determinant to perform better than FreeSurfer's affine-based estimate.

Since (a) we have shown that SPM12 can provide significantly better estimated TIV than FreeSurfer's affine-based estimate, and (b) we have discussed the limitations of adjustment by division, nonlinear-only modulation may be seen as a convenient but possibly suboptimal procedure.

Outliers and algorithmic failures

It is often the practice when dealing with automated methods to screen for failures and attempt alternative methods to recover from them (either adjustments to parameters or resorting to another technique). Doing so raises questions of how to define and screen for failures, as well as whether any adjustments should be incorporated as improvements to the method used. On inspecting the range of results we decided to regard as suspicious any estimated TIV greater than 3000 ml and that omitting the two clear failures in the case of FreeSurfer was a fairer comparison of “out-of-the-box” performance. Omission of the corresponding values from the SPM comparisons has little effect on results (Supplementary Table 1). Additionally we attempted to fix the two FreeSurfer failures by suppressing the automated registration checking using the “-notal-check” option for those cases and results for these can be found in Inline Supplementary Table 1. The absence of any notable outliers for SPM12 suggests a good degree of algorithmic robustness.

Only the numerical results were screened for outliers, no routine quality control was applied to the individual images.

Limitations and further work

We have not directly investigated the effect of atrophy on TIV estimates (cf. Pengas et al., 2009). Our sample population were all probable-AD trial participants with varying degrees of atrophy at baseline (mean manual brain volume:manual TIV ratio 0.69 ± 0.05, min 0.56 max 0.81); the high agreement with manual measures shown by SPM12 is apparently not affected by this. Both Nordenskjöld et al. (2013) and Ridgway et al. (2011) find only small variability over time for manual TIV measures. However, further longitudinal evaluation (as in Pengas et al. (2009)) of the SPM12 method could provide greater reassurance that tissue loss does not change TIV estimates.

Although a large number of different sites and scanners were included here, without obvious detriment to the results, all scanners were 1.5 T, so we cannot claim to have demonstrated robustness to different field strengths (cf. Keihaninejad et al., 2010).

Whilst we have shown very good agreement between SPM12 and manually-measured TIV, we have not directly evaluated the effect of replacing manual with SPM12 values in investigations of other factors such as hippocampal volume as in Nordenskjöld et al. (2013), though the simulations in Supplementary Fig. 2 shed some light on this.

It would also be of interest to evaluate SPM12's performance on data other than T1-weighted MRI; for example, Pengas et al. (2009) and Nordenskjöld et al. (2013) favour proton-density (PD) weighted imaging, whilst Vuong et al. (2013) shows the advantages of T2-weighted MRI. Multi-spectral segmentation of quantitative multi-parametric maps Weiskopf et al. (2013) would be expected to yield even better results, since the combination of PD with other contrasts (including R2*, related to T2) should enhance the distinctions between brain tissue and blood-filled sinuses, and between the CSF and bone/air. It is plausible that automatic methods, perhaps with further refinements, could actually yield more accurate measurements than the current manual gold standard, especially with the latter's use of only every 10th slice, however, demonstrating this would be challenging, requiring somewhat indirect evaluation.

A.A. Old and new segment in SPM8

The “old segment” (OS) algorithm of SPM8 is essentially the same as the default tissue segmentation procedure in SPM5 (Ashburner and Friston, 2005). The “new segment” (NS) algorithm of SPM8 is based on the same principles, although some changes were made to the original algorithm; these changes are briefly described in Appendix A of Weiskopf et al. (2011).

Both implementations produce probabilistic segmentations of images into C tissue classes, where each class is defined by a tissue probability map that is warped into alignment with the image. The C = 4 tissue probability maps of OS (grey matter, white matter, CSF and other) were extended in NS to include tissue priors for bone, non-brain soft-tissue and air (C = 6), which resulted in greater robustness in terms of registering the tissue priors to the image. The other main extensions were an increased flexibility of the registration part of the model (by using more parameters), and the ability to do multi-spectral segmentation by simultaneously modelling multiple images of the same subject (eg. T2-weighted and PD-weighted).

When modelling a single image, both versions estimated a correction for intensity nonuniformity, parameterised by coefficients β. In what follows, we denote the nonuniformity correction at voxel i by p_i(β). Intensity distributions for each tissue class are modelled by a mixture of Gaussians, such that the means and variances are adjusted to account for the correction field.

Each of the K_c Gaussians of each tissue class is described by its mixing proportion (γ_ck, such that γ_ck ≥ 0 and ∑_k=1^K_C γ_ck = 1), a mean (μ_ck) and a variance (σ_ck²). With this model, the probability of observing a voxel of intensity y_i in the image is given by:

In both cases, warping the tissue maps into alignment involves estimating a vector of coefficients (α) that parameterise displacement fields. The ith voxel of the warped version of the cth tissue prior is denoted by b_ic(α). For NS, the likelihood model that is maximised (with suitable regularisation) is simply:

For OS, there is effectively an additional set of C scaling parameters that are estimated, η, such that η_c ≥ 0 and ∑_c=1^C η_c = 1. This model accounts for situations where there may be globally more or less of some particular tissue type.

A.B. Segmentation in SPM12

The absence of these scaling parameters has impacted the behaviour of NS in SPM8. This issue has been pointed out numerous times on the SPM mailing list, and is the likely cause of some of the findings in Callaert et al. (2014) as well as Peelle et al. (2012). The issue has now been resolved in the beta version of SPM12, as the segmentation algorithm of SPM12 is essentially NS, but with the η parameters included again. Other changes from NS to SPM12 include the new tissue probability maps (described earlier) and a change to the regularisation of the spatial transformation model. The single-parameter bending energy regularisation mentioned in Appendix A of Weiskopf et al. (2011) has been extended to a more sophisticated regularisation model in SPM12, which now has five penalty terms: absolute displacement, membrane energy, bending energy, linear elasticity and divergence; though the bending energy has the largest weight by default. SPM12 also reintroduces a version of the morphological clean-up procedure from OS, and includes a Markov Random Field (MRF) based clean-up that was introduced in later versions of NS in SPM8, but was not documented in Weiskopf et al. (2011). For the present work, the key changes are the tissue priors and the scaling parameters.