Modelling the distribution of white matter hyperintensities due to ageing on MRI images using Bayesian inference

White matter hyperintensities (WMH), also known as white matter lesions, are localised white matter areas that appear hyperintense on MRI scans. WMH commonly occur in the ageing population, and are often associated with several factors such as cognitive disorders, cardiovascular risk factors, cerebrovascular and neurodegenerative diseases. Despite the fact that some links between lesion location and parametric factors such as age have already been established, the relationship between voxel-wise spatial distribution of lesions and these factors is not yet well understood. Hence, it would be of clinical importance to model the distribution of lesions at the population-level and quantitatively analyse the effect of various factors on the lesion distribution model. In this work we compare various methods, including our proposed method, to generate voxel-wise distributions of WMH within a population with respect to various factors. Our proposed Bayesian spline method models the spatio-temporal distribution of WMH with respect to a parametric factor of interest, in this case age, within a population. Our probabilistic model takes as input the lesion segmentation binary maps of subjects belonging to various age groups and provides a population-level parametric lesion probability map as output. We used a spline representation to ensure a degree of smoothness in space and the dimension associated with the parameter, and formulated our model using a Bayesian framework. We tested our algorithm output on simulated data and compared our results with those obtained using various existing methods with different levels of algorithmic and computational complexity. We then compared the better performing methods on a real dataset, consisting of 1000 subjects of the UK Biobank, divided in two groups based on hypertension diagnosis. Finally, we applied our method on a clinical dataset of patients with vascular disease. On simulated dataset, the results from our algorithm showed a mean square error (MSE) value of 7.27×10−5, which was lower than the MSE value reported in the literature, with the advantage of being robust and computationally efficient. In the UK Biobank data, we found that the lesion probabilities are higher for the hypertension group compared to the non-hypertension group and further verified this finding using a statistical t-test. Finally, when applying our method on patients with vascular disease, we observed that the overall probability of lesions is significantly higher in later age groups, which is in line with the current literature.


Introduction
White matter hyperintensities of presumed vascular origin (WMH), also called white matter lesions (Wardlaw et al., 2013), are common findings in the brain white matter that appear hyperintense on T2-weighted, fluid attenuated inversion recovery (FLAIR), and proton density-weighted brain MRI images. Even though the pathogenesis of WMH has not yet been well understood (Wardlaw et al., 2013), WMH are strongly associated with cerebrovascular disease and vascular risk factors (Li et al., 2013), and they are also frequently found in neurodegenerative diseases such as Alzheimers disease (Debette and Markus, 2010;Prins and Scheltens, 2015). In the general population, they have been associated with increased risk of stroke, dementia and death .
Several visual rating scales for WMH are available and commonly used (Wahlund et al., 2001;Fazekas et al., 1987). Various segmentation algorithms for voxel-wise assessment of WMH are also available (Caligiuri et al., 2015), including our own (Griffanti et al., 2016). They generate a lesion probability map for a single subject indicating the probability for each voxel of being WMH. However, determining the voxel-wise distribution of WMH at the population-level is important in order to investigate the relationship between the spatial distribution of lesions and various factors such as age, vascular risk factors and cognitive function (e.g. measured with the Montreal Cognitive Assessment (MoCA) score). In a population-level parametric lesion probability map, the value at each voxel reflects the probability of finding a lesion in a specific population for a parametric factor of interest. This would aid in studying differences in lesion distribution between normal and pathological ageing (Biesbroek et al., 2013;Duering et al., 2014) and be useful in the clinical context of the studied population. In fact, although the occurrence of WMH is common in the ageing population, clear relationships between amount and distribution of lesions and the cognitive, demographical, and risk factors are not yet well understood.
The simplest method used in the literature for correlating voxel-wise lesion distribution maps with various parametric factors (such as age, cognitive function and so on) includes averaging subject-level lesion probability maps (Rostrup et al., 2012) that are grouped according to those parametric factors. However, this approach does not necessarily provide continuity or smoothness in space or across the parametric dimension for the lesion probabilities. Alternatively, several mass-univariate methods such as voxel-based lesion-symptom mapping (VLSM) (Bates et al., 2003) and voxel-wise linear regression of lesion probabilities against various factors (Charil et al., 2003(Charil et al., , 2007 have been proposed. For example, Charil et al., 2007 studied the correlation between the cortical thickness and various factors including lesion distribution in multiple sclerosis using voxel-wise linear regression. However, these methods are based on models that work independently in each voxel and do not capture the relationship between neighbouring voxels. Hence, they are not ideal for modelling the lesion distribution, since lesions are clustered regions rather than isolated voxels. Moreover, the linear models used in these methods are not optimal for binary data (Charil et al., 2007;Bates et al., 2003).
Spatially varying coefficient processes establish a local spatial relationship between the coefficients of regression models (Gelfand et al., 2003;Gamerman et al., 2003;Ge et al., 2014). For example, Bayesian Spatial Generalized Linear Mixed Model (BSGLMM) proposed by Ge et al. (2014) is based on spatially varying coefficients to determine the relationship between the spatial distribution of lesions and subject specific covariates such as multiple sclerosis (MS) subtype, age, gender, disease duration and disease severity measures. The spatially varying coefficients are modelled jointly using a multivariate pairwise difference prior model, a particular instance of the Multivariate Conditional Autoregressive model. However, the dimension of the correlation matrices involved in spatially varying coefficient processes is very high and their inversion becomes computationally infeasible for very large imaging datasets (Ge et al., 2014). In fact, the computational load of the model proposed by (Ge et al., 2014) requires the use of graphical processing units for parallel computing.
Several solutions have been proposed to overcome this problem, including Gaussian predictive processes (Banerjee et al., 2008), using a truncated Karhunen-Love expansion to estimate spatially varying coefficients (Crainiceanu et al., 2009), functional principal component analysis approach (Reiss and Ogden, 2010) and covariance tapering (Kaufman et al., 2008). However, all of them perform data reduction (unlike Ge et al., 2014) and depend on strict assumptions (Crainiceanu et al., 2009).
In this paper, we propose an alternative approach by developing a probabilistic model to obtain a parametric lesion probability map in order to investigate the relationship between distribution of WMH and various parametric factors. Our method is not computationally intensive, and uses a spline representation to ensure continuity and smoothness between neighbouring voxels in both spatial and parametric dimensions. In our model, the lesion distribution at each voxel is obtained by the linear combination of all the splines on the neighbouring voxels. Moreover, our model adopts a Bayesian framework in order to overcome the discretely sampled nature of the input binary maps obtained from lesion segmentation. Our approach for modelling the spatial distribution of lesions is suitable for very large imaging datasets such as UK Biobank (Miller et al., 2016) and allows the flexibility of observing the effect of various parametric factors on the lesion distribution. Although in principle our model could work with any factor, in this work we analysed the voxel-wise distribution of lesions with respect to age since a strong relationship between the progression of WMH and age has already been established (Simoni et al., 2012).
We first evaluated our model on a simulated dataset and compared it to existing methods. We then validated our model on two real datasets. The first one is a subsample of the UK Biobank, in which we tested the relationship between WMH and age comparing subjects with and without hypertension. We chose hypertension as our grouping variable, since it has been found to be one of the important risk factors of WMH (Dufouil et al., 2001;Gottesman et al., 2010) in addition to age. Also in this case we compared our results against the existing methods. As a final validation, we applied our model to a clinical population (vascular population of subjects who had a transient ischemic attack or minor stroke) and analysed the results with respect to age.

Modelling the distribution of WMH using Bayesian inference
Our proposed Bayesian spline model takes as input the binary lesion map and the age (or other parametric factor of interest) for each individual subject. We model the distribution of lesions in three steps: 1) constructing a representation for the lesion probability distribution using by spline basis functions, 2) formulation of the posterior lesion probability; and 3) maximisation of the posterior probability by constrained optimisation.
Let D s be the 3D binary lesion map of individual subject S. Our aim was to form a 4D spatio-temporal parametric lesion probability map with the 4 th dimension indicating the parametric factor of interest.

Constructing a representation for the lesion probability distribution using by spline basis functions
Since the binary lesion maps contains localised areas, the lesion voxels are sparse and have discrete values (0 or 1). Therefore, to ensure spatial and temporal continuity in the probability distribution of lesions, we approximate the lesion probabilities using cubic b-splines.
Consider cubic b-splines with spline coefficients C i and basis functions B ij , where i denotes the indices of knot points of the basis functions and j ¼ ðx; y; z; tÞ indicate the spatial coordinates in the first 3 dimensions and the parametric factor in the 4 th dimension. The lesion probability θ j at each voxel is related to C i and B ij by eqn. (1): The probability value at each voxel is calculated as the linear combination of spline basis functions. We formulated the splines as a separable outer product of individual 1D splines in 4 dimensions (3 in space and one in the parametric dimension). To implement the above calculation, we generated 1-dimensional cubic b-splines and convolved the binary map with scaled 1D b-splines independently in all 4 dimensions to get spatially continuous lesion probabilities θ as shown in Fig. 1. This operation is more computationally efficient than off-the-shelf spline fitting toolboxes and is well suited to large datasets.
It is worth noting that after each convolution edge artefacts occurred due to convolution with the tails of the splines and zero padding. In order to get the corrected convolution output f norm , we divided the convolution output f in with a normalising image (formed by convolving the same bsplines B with an unity image uðxÞ, where uðxÞ ¼ 1 inside the valid FOVbrain voxels -and 0 outside, thus forming a similar pattern at the edges) as shown in eqn. (2).

Formulation of the posterior lesion probability
The most common lesion probability estimation method is averaging lesion maps across subjects. However, the accuracy of this approach is sensitive to the amount of data, especially considering the fact that there might not be any subject representing specific age groups or a few subjects might not have any WMH. Bayesian methods are well suited to this problem since they allow for the uncertainty associated with the limited amount of data, while the spline model ensures continuity between the neighbouring voxels.
As we consider age as our parametric factor, the 4 th dimension, t, of our data indicates the bin number corresponding to the age groups. We formed these bins t by grouping the subjects into age groups age t (e.g., age 1 ¼ [20, 22], age 2 ¼ [23, 25], etc.). The age groups can be defined to have a duration as short as we require (even in months) and hence binning the ages does not necessarily restrict the values to be overly discretised in the parametric dimension. Let R (shown in Fig. 2(a)) be the 4D spatio-temporal volume formed by summing the binary lesion maps of the subjects in each age group age t . The value at each voxel R j ¼ Rðx; y; z; tÞ denotes the number of subjects having lesions at that voxel in the age group age t . Let N be a 4D volume in which the value at each voxel N j ¼ Nðx; y; z; tÞ denotes the number of subjects in each age group age t . R and N are formed by, For an age group age t , R j indicates number of observations of binary data (lesion occurrence) in the subjects belonging to this age group age t , with the total number of subjects in age t given by N j . Thus the probability of observing R j number of binary outputs given N j and probability of lesion occurrence θ j for an age group age t is given by the binomial likelihood distribution, Rj!ðNjÀRjÞ! . The full likelihood distribution (assuming independence between age groups and voxel locations) for all age groups and at all voxel locations is given by the product of individual likelihoods, Using Bayes theorem, the posterior lesion probability distribution pðθjN; RÞ is given by, pðθjN; RÞ∝pðRjN; θÞ pðθjNÞ (8) We will assume that the prior probability of lesion occurrence is the same for any point in space and time. Hence with uniform prior pðθÞ ¼ 1, Our aim is to maximise pðθjN; RÞ to obtain a parametric model of the lesion probability distribution over a population. To make the calculation of derivatives simpler, we determine the log-posterior function since the logarithm is a monotonically increasing function. Log-posterior LðθjN; RÞ  is given by, Since the lesion probability values θ have been approximated with spline basis functions as explained in the section above, we substitute the values of θ from eqn. (1) in the above log-posterior equation, The derivative of LðθjN; RÞ with respect to the spline coefficients C k is given by, Note that the derivative function requires more spline modelling steps (basically, each summation requires a 4D spline approximation), however, performing 1D convolutions independently in 4 dimensions considerably reduces the running time.

Maximisation of the posterior probability by constrained optimisation
In order to estimate the final lesion probability θ i , we need to determine the spline coefficient values that would maximise the log-posterior. Also, by maintaining the value of C i within the range [0,1] the final θ i values will be constrained to the range [0,1]. Hence, to determine the value of C i corresponding to the maximum posterior estimate, we formulated it as a constrained local optimisation problem as specified below: We used the steepest gradient descent algorithm, which is a first order optimisation algorithm that uses an iterative polynomial line search method for the above optimisation. The step size γ for each iteration is determined by a polynomial line search method satisfying Armijo -Wolfe conditions (Wolfe, 1969). In our case, the function to be minimised is À LðθjN; RÞ. For L : ℝ n → ℝ, the step size for each iteration follows the following condition: where, θ k is the current guess, q k indicates the search direction and γ is the step size. Moreover, in order to speed up the convergence, we redefined C i using a normalising function (with parameter λ i 2 ð À ∞; þ ∞Þ) in order to maintain the range of C i within [0,1] and relax the optimisation constraints as given by, Using the chain rule, the derivative of LðθjN; RÞ (eqn (16)) can be modified as We perform a step of spline modelling (given by eqn. (1)) on the optimal C i that maximises the log-posterior ÀLðθjN; RÞ to get the final estimated parametric lesion probability map b θ. We determined the average of lesion binary maps of subjects within the age group age t to get the age group-wise average 4D lesion map Average ¼ P j ðR j =N j Þ. We provided smoothed Average as the initial estimate for C i during optimisation. Fig. 3 illustrates the initial and final values of C i and λ i involved in the optimisation step for two different age groups (taken from different age bins of the 4D spatio-temporal map). The initial θ values (Fig. 3(d)) and the final b θ values (Fig. 3(e)) are shown with their difference (Fig. 3(f)).
The results show little difference between the initial θ values and final b θ values for the younger age group ( 0.005 in Fig. 3(f, bottom row)) and a higher difference (% 0.01 in Fig. 3(f, top row)) in the elder age group. This could be attributed to the difference in the amount of lesion voxels between the age groups.

Convergence analysis
Implementation of steepest gradient descent optimisation algorithm was done in MATLAB (Bortz and Kelley, 1998). We specified the initial step-size γ to be 1 Â 10 À3 , while step sizes for the subsequent iterations were determined by the polynomial line search algorithm, for which lower and upper bounds were set to 0:1 Â initial γ and 0:5 Â initial γ respectively. We specified the maximum number of iterations as 50. The main parameter for determining γ for each iteration is the 2-norm kG k of the 4D gradient matrix G (calculated from eqn. (20)). We set the final convergence tolerance value for the decrease in kG k as 1 Â 10 À4 (since gradient of the function must be zero at the function minimum).

Simulations
In order to evaluate the accuracy of our model, we initially tested it on a simulated dataset. To this aim, we compared the ground truth distribution (θ true ) with the probability distribution ( b θ) estimated by our algorithm and compared it with other methods.
We simulated volumes of size 91 Â 109 Â 91 arranged in 6 months age bins to form 60 age groups, thus creating a simulated set of subject images with dimensions 91 Â 109 Â 91 Â 60. For this simulation, we assigned a random number of subjects in each age group N sim 2 ½1; 50. For each subject, we generated a lesion map by randomly sampling voxels within a specified ROI (in our case, the MNI brain mask). For sampling, we used the smoothed average of lesion probability maps from a previous study (vascular cohort in Griffanti et al. (2016)) as the ground truth probabilities to obtain a realistic WMH spatial distribution. The probability of sampling each lesion voxel in the data is initially represented by a random normal variable NðμðxÞ;1Þ. In order to avoid isolated lesion voxels and to establish spatial dependencies between the neighbouring voxels (so that sampled voxels would look like plausible lesions), we smoothed the sampled lesion voxels using a Gaussian kernel K with a standard deviation of 0.8. We obtained binary lesion maps by thresholding the smoothed lesion map K*NðμðxÞ; 1Þ above 0. As a result of convolution with K, the ground truth distribution at voxel x is now represented by NðK*μ; K 2 *MÞ, where M is a binary brain mask. We obtained the true lesion probability for the simulation, θ true , by calculating the cumulative distribution function of the normal variate z ¼ K*μ= ffiffiffiffiffiffiffiffiffiffiffiffi ffi K 2 *M p or more specifically, pðz > 0Þ. Repeating the above process N sim times for each age group, we accumulated the sum of binary volumes as and applied our algorithm to R sim and to model the distribution of lesions over the simulated population, to get b θ (algorithm output). An instance of data obtained from the simulation is shown in Fig. 4. We applied our algorithm to the simulated data and compared its performance with respect to some alternative methods with different levels of complexity. The simplest method for estimating the lesion probability with respect to age groups is to determine the group-level average of lesion probability maps as done in Rostrup et al. (2012). We performed the similar average on binary lesion maps and calculated Average ¼ R t =N t for each age group age t . However, the resulting lesion probability map Average is sensitive to the number of subjects in each age bin age t and does not provide a continuous/generalized estimate of lesion distribution. Hence another option is to smooth the average lesion distribution map Average, using a Gaussian kernel with different standard deviations σ (0.5, 1.5 and 3.0) to get a more spatio-temporally continuous lesion probability map Smoothed AverageðσÞ. We also applied the Bayesian Spatial Generalized Linear Mixed Model (BSGLMM) based on spatially varying coefficients proposed by Ge et al. (2014). We evaluated the performance of the above methods by determining the error values (difference between the ground truth probabilities and the outputs of the above methods) and the corresponding mean-squared error (MSE) values.

Application to real data
We applied our algorithm to two different groups from the UK Biobank data and compared the lesion distribution between the two groups using the methods that provided better MSE values than those reported in the literature, as determined from tests on the simulated dataset. We further validated our observations by applying our method on a clinical dataset consisting of a vascular population of subjects who had experienced a transient ischemic attack or minor stroke.
The dataset used for initial validation is a subset of UK Biobank data (Miller et al., 2016). In the UK biobank data, about 10% of the sample has been identified as being diagnosed with hypertension, corresponding to around 7975 subjects. Our goal was to model the lesion distribution with respect to age, comparing the subjects diagnosed with hypertension (HT group) to subjects without hypertension, where the latter acts like a control group (non-HT group). In this case, to reduce the amount of computation, we selected a balanced subset of randomly sampled 1000 subjects, 500 subjects from the HT group (age range 45.5-78.3, mean age 66.3 AE 6.1 years, with female to male ratio, F:M ¼ 202:298), while the remaining were from the non-HT group (age range 45.5-78.4, mean age 62.0 AE 7.6 years, with female to male ratio, F:M ¼ 257:243).
For these subjects we generated binary lesion masks to be used as input for our algorithm. This was performed using BIANCA (Griffanti et al., 2016) on FLAIR images, also using information from T1-weighted images (Alfaro-Almagro et al., 2018). The lesion probability maps obtained were binarised after thresholding at 0.8. We then registered the single subject binary lesion maps to the MNI space using linear and non-linear registration, thresholded them at 0.5 and binarised them to compensate for interpolation. We then applied our algorithm to the two groups separately and compared the results with those from other methods. When applying our method, we arranged the subjects in each group into one year age bins to form a 4D image of dimension 91 Â 109 Â 91 Â 33. We compared the two groups using an unpaired t-test using permutation testing (randomise, Winkler et al., 2014) to analyse the significance of the difference in the lesion probability values between HT and non-HT groups, and calculated the corresponding z-scores. Results were considered significant for p corr < 0:05, corrected for multiple comparisons by using the null distribution of the max (across the image) voxelwise test statistic. For BSGLMM, the t-values were obtained by dividing the posterior mean by posterior standard deviation (Ge et al., 2014).
Finally, we validated our algorithm on a clinical dataset of subjects at risk of vascular cognitive impairment. The dataset consists of MRI data from 474 consecutive eligible participants from Oxford Vascular Study (OXVASC) (Rothwell et al., 2004), who had recently experienced a minor non-disabling stroke or transient ischemic attack (TIA) (age range 20-102 years, mean age 67.4 AE 14.3 years, with female to male ratio, F:M ¼ 240:234) (see Griffanti et al., 2016 for more details). Similar to the UK Biobank dataset, we obtained the lesion binary masks using BIANCA and registered the single-subject masks to the MNI space. We grouped the data in age bins of 3 years to form a 4D image of dimensions 91 Â 109 Â 91 Â 27 and applied our model. We also performed permutation-based statistical analysis (corrected for multiple comparisons) to test the statistical significance of the effect of age on the lesion probabilities. For this analysis, in each iteration, we permuted the data with respect to the age groups and applied our model on the permuted data to generate results for that permutation, that were then used to assess statistical significance, with correction for multiple comparisons (using the maximum statistic across space).   Ge et al., 2014). Moreover, both MSE values and running times increases as the standard deviation σ of the smoothing kernel increases due to increase in the kernel size. Fig. 6 shows the lesion probability with respect to age at two  representative voxels (one more anterior and one more posterior) for the tested methods. The curves obtained from Average and Smoothed Averageðσ ¼ 0:5Þ are very close to the ground truth, but they show spikes due to their dependence on the number and WMH characteristics of the specific subjects in each group and our method give curve that capture the true trend quite well. BSGLMM posterior probabilities tend to increase monotonically with age throughout the white matter irrespective of the location of the lesions. We also tested the effect of change in cubic b-spline knot spacing and the standard deviation of kernel K used in the ground truth estimation on the error values (see supplementary material). We observed that the cubic b-splines with knot spacing of 2 voxels provided the best result on the simulated data and hence we used the same specifications henceforth on the real data experiments.

Results on the real data
We show the results of our algorithm on UK Biobank data to estimate the lesion distribution probabilities with respect to age within the HT and non-HT subjects in Fig. 7.
The lesion probabilities increased with age for both the groups and are higher in the HT group (% 0.4 as shown in Fig. 7(c)) than the non-HT group. We can observe that this difference is higher in the periventricular regions compared to deep regions, and higher in the older age groups (% 0.3 between the columns in row (c)) compared with the younger group.
So not only is there a strong relationship between lesion probabilities and hypertension, but also this relationship is more pronounced as the population gets older.
In Fig. 8 and Fig. 9, we show the results of our algorithm compared with the posterior probabilities obtained from the BSGLMM method and the outputs of smoothed average maps (Smoothed Averageðσ ¼ 0:5Þ, Smoothed Averageðσ ¼ 1:5Þ) respectively. We chose to compare the above methods with our Bayesian spline algorithm since these methods gave MSE values lesser than 12 Â 10 À5 (Ge et al., 2014) on the simulated data (Table 1). We observed that the lesion probabilities in periventricular regions increased with age for both HT and non-HT groups, while the probabilities in deeper regions increased only in the HT group (between the first two rows in (a) in Figs. 8 and 9). The bottom rows of Figs. 8 and 9 show the z-scores for significant regions (p corr < 0:05), obtained from a statistical test to determine the significance of differences between the HT and non-HT groups. We observed increased lesion probabilities in the HT group compared with non-HT group, especially in the deep regions. This effect is particularly strong in our method and Smoothed Averageðσ ¼ 1:5Þ when compared with BSGLMM and Smoothed Averageðσ ¼ 0:5Þ (Figs. 8 and 9 (bottom row)). Fig. 10 (top row) shows the estimated lesion probabilities from our Bayesian spline method on OXVASC data for two representative age groups. Similar to our results on UK Biobank data, it can be observed that the overall lesion probability values are relatively higher in the older age group compared to the younger age group (% 0.2, between second and Fig. 7. Results on UK Biobank data. Output of our Bayesian spline method on HT group (a) and non-HT group (b) shown along with their difference maps (c). First column shows the age group t ¼ 20 (age t ¼ 64.4-65.4 years) and second column shows t ¼ 33 (age t ¼ 77.4-78.4 years). All the results have been shown at the slice z ¼ 45 (MNI template on the left). third columns). On performing the permutation test, we observed that periventricular voxels change significantly with age (p < 0:05) as shown in Fig. 10 (bottom row).

Discussion and conclusion
In this work we compared several methods including our proposed Bayesian spline method to model the distribution of white matter lesions with respect to a parametric factor within different populations. For the experiments in this paper, we considered age as our parameter of interest, but the framework is general and can work with any continuous parameter of interest. For example, we also modelled the distribution of white matter lesions with respect to MoCA scores (refer supplementary material for details).
During the development, we analysed the effect of various parameters of our algorithm on the resulting lesion probabilities. One of the most important parameters in the algorithm is the knot spacing in the spline model. We studied its effect and how it relates to the size of smoothing kernel K used for the generation of the simulation data by varying them independently and evaluated their combined effect on the algorithm result b θ (see supplementary material). The best results were obtained for the knot spacing of 2 with MSE of 7:270 Â 10 À5 . We also evaluated the effect of the initialisation provided for the optimisation function and observed that the algorithm converges quicker when we provide the smoothed average lesion map as the initial value for λ i , whereas it fails to converge when a zero image or unity image is provided as initialisation Fig. 8. Comparison of the results of the Bayesian spline algorithm with BSGLMM on UK Biobank data. Top two rows show the results of the methods for HT and non-HT groups for two age groups. The bottom row shows the z-scores in the significant regions (p corr < 0:05), obtained from the significance test on the differences of the lesion probabilities between HT and non-HT groups obtained by the corresponding methods. show the results of the methods for HT and non-HT groups for two age groups. The bottom row shows the z-scores in the significant regions (p corr < 0:05), obtained from the significance test on the differences of the lesion probabilities between HT and non-HT groups obtained by the corresponding methods.
for λ i . Also, our algorithm converges to a similar output when the average lesion maps are smoothed with kernels of similar standard deviation. We additionally verified on real data that our algorithm is robust with respect to changes in the initialisation due to various levels of random perturbations (up to 30%) in the group-wise average 4D map R i (see supplementary material for details).
When comparing the results of different methods on the simulated dataset, we found that the MSE values are lower than that reported in Ge et al. (2014) for most of the methods, except Average and Smoothed Averageðσ ¼ 3:0Þ. We obtained the highest MSE value and highest speed for the simplest method Average. However, this also gave a resulting image that was not as spatially continuous as the other methods. Moreover, the accuracy of Average (as well as of its Gaussian smoothed results Smoothed AverageðσÞ) depends on the number of subjects in the age groups. While they provide reliable estimates in larger populations, their estimates are not so accurate when there are fewer subjects or when the lesion load is lower, since the data becomes sparse. Also, in the case of smoothed Smoothed AverageðσÞ method, determining the optimum value of standard deviations σ for the Gaussian kernel is dependent on the dataset and the lesion characteristics. However, the smoothed Smoothed AverageðσÞ outputs provide good initial estimates for our method and lead to better convergence, as we discussed above. We observed that, while results of both the Bayesian spline method and BSGLMM show relatively similar MSE values, the spline modelling in our algorithm is spatially more robust when compared to BSGLMM method. For instance, from Fig. 6, while lesion probabilities obtained from BSGLMM method fit well in periventricular areas, they do not fit the ground truth distribution as well in deep areas (Fig. 6). Another important advantage of our method is that it is much faster than BSGLMM (Table 1) and can easily scale up to very large numbers of subjects. Our modelling algorithm takes less than 20 min (84 s per iteration), when run on an iMac with 2.9 GHz Intel Core i5 processor. This is due to the fact that we independently perform the convolution in all 4 dimensions with a 1D spline which reduces the computational load. Also, by calculating the analytical gradient (Eqn. (16)) we avoid the numerical computation of gradient for 4D data in the optimisation step. Even though the maximum number of iterations allowed for the optimisation step is 50, our algorithm typically converges earlier (10 iterations on simulated data and < 10 on real data). On the other hand, the BSGLMM algorithm is computationally demanding since it has to estimate the posterior distribution via Markov Chain Monte Carlo model using Gibbs sampler at each voxel.
Our results on the real data (UK Biobank and OXVASC datasets) show that the lesion distribution probabilities are higher in the later ages, which is in line with the existing literature (Simoni et al., 2012;Alfaro-Almagro et al., 2018). Moreover, the fact that we observed an increase in lesion probabilities in deep white matter regions in the HT group with respect to the non-HT group for all the methods (Figs. 8 and 9), suggests that deep lesions might be associated with hypertension. The same phenomenon has also been reported in Rostrup et al. (2012). We further verified it with our statistical tests and observed higher value of z-scores with greater significance (p < 0:05), particularly in the deep regions between the HT and non-HT groups. Though this is generally observable in all the methods, it is more evident for our Bayesian spline method and smoothed average map Smoothed Averageðσ ¼ 1:5Þ i.e. deep regions are affected significantly in the HT group compared with the non-HT group for these two methods.
The results on the UK Biobank data are similar to those on the simulated data when comparing the Bayesian spline method with BSGLMM and Smoothed Average. Smoothed Averageðσ ¼ 0:5Þ, did not yield a smooth lesion distribution and had isolated voxels in the lesion probability maps. Since this method gave lower MSE on simulated data with respect to our method, and in absence of ground truth, we cannot exclude the possibility that spline interpolation in our method induces spurious smoothness that is not inherent to the data. However, the simulation results also suggest that Smoothed Averageðσ ¼ 0:5Þ method produces less smooth data because it is affected by the sparsity of the observed data and hence is dependent on data/lesion characteristics. Also, results of the statistical test using Smoothed Averageðσ ¼ 0:5Þ show only a few isolated voxels that are significantly different between the HT Fig. 10. Results on OXVASC data. (Top row) From left to right: MNI brain, lesion probability distribution at bin numbers t ¼ 3 (age t ¼ 29-31 years) and t ¼ 13 (age t ¼ 59-61 years). The results are shown at slice z ¼ 45. (Bottom row) Results of permutation analysis. From left to right: Voxels that change significantly (p < 0:05) with age (red) shown on slices z ¼ 39; 45 and 49. group and the non-HT group, while the relationship between WMH and blood pressure is well known (Dufouil et al., 2001;Rostrup et al., 2012) and detected quite strongly by the other methods (including ours), that produce a smoother lesion distributions.
To conclude, we compared different methods to model the distribution of white matter hyperintensities within a population with respect to a parametric factor of interest, such as age. The evaluation of the proposed Bayesian spline method on the simulated data provides an MSE value of 7:27 Â 10 À5 , with the advantage of being computationally more efficient. We compared the results of various methods on a real dataset between hypertension and control groups and found that deep lesions increase significantly with hypertension, which is inline with existing work. We validated our model on real data showing that periventricular lesions increase significantly with age, which is well established in the literature. The results are consistent between the two datasets, showing that our model is robust across datasets. We also verified that this trend of change in lesion probabilities with age remains same irrespective of the initial thresholds use for obtaining subject-level lesion binary maps (see supplementary material for details). As future work, the Bayesian spline method could be used to model the voxel-wise relationship between the distribution of WMH and other population-level parametric factors that can be useful in the clinical context of the studied population. Another interesting future direction could be to extend the scope of our method in modelling the distribution of other types of lesions (e.g. Multiple sclerosis) and test the robustness of the spline method in modelling lesions that have a different spatial pattern from WMH. In addition, the lesion probability distribution map obtained from our model can be used either as an additional feature or as a spatial prior to improve lesion segmentation algorithms (Sundaresan et al., 2018).