An information‐based comparison of diffusion attenuation models in normal and inflamed bone marrow

Diffusion‐weighted imaging has received attention as a method for characterizing inflammatory exudates in bone marrow in immune‐mediated inflammatory diseases and reveals an increase in diffusivity in regions of bone marrow oedema. Various models of diffusion attenuation have been investigated but the model providing the best description of tissue pathophysiology in regions of marrow oedema is unknown. Determining the most appropriate model is an important step towards protocol optimization and the development of a robust and clinically useful method.


Funding information
Determining the most appropriate model is an important step towards protocol optimization and the development of a robust and clinically useful method.
We aimed to determine which of three candidate models of diffusion attenuation most accurately describes the acquired signal from normal and inflamed bone marrow. 11 subjects with spondyloarthritis and evidence of active inflammation (ie bone marrow oedema) on MRI and 17 patients with no evidence of active inflammation underwent diffusion-weighted imaging of the sacroiliac joints (b-values 0, 50, 100, 300 and 600 s/mm 2 ). Monoexponential, intravoxel incoherent motion (IVIM) and kurtosis models were fitted to the acquired signal from regions of interest in areas of bone marrow oedema and normal marrow. The three models were compared in terms of sum of squared error and information content (corrected Akaike information criterion). Model parameters were compared between regions of bone marrow oedema and regions of normal marrow. f the three models investigated, the IVIM model provided the best description of the signal over the 0-600 s/mm 2 range across normal and inflamed bone marrow. There was a particular advantage of the IVIM model in normal marrow, where it was best able to capture the pronounced fast diffusion component observed in several cases. However, IVIM and kurtosis effects both became smaller and the signal behaviour became closer to monoexponential in the presence of bone marrow oedema. Our data suggest that increases in D tissue (in the IVIM framework) might account for the reduced deviation from monoexponential behaviour in oedematous bone.

| NTRODUCTION
Diffusion-weighted imaging (DWI) has been widely used for interrogating structure, physiology and pathology in the central nervous system, 1-4 but it has only been used to image the body 5-7 relatively recently. A promising emerging application is the use of DWI for imaging immune-mediated inflammatory diseases, in which dysregulation of the normal immune response leads to inflammation and ultimately end-organ damage. Examples of these diseases include the inflammatory arthritides such as rheumatoid arthritis and spondyloarthritis, vasculitis, myositis, inflammatory bowel disease and various other disorders affecting a variety of systems. In the spondyloarthritidesinflammatory diseases defined by the presence of inflammation and/or new bone formation involving the spine-DWI has been used to characterize areas of bone marrow oedema, where inflammation leads to an increase in the number of inflammatory cells and probably to an increase in extracellular water. [8][9][10][11] In areas of bone marrow oedema there is an increase in the apparent diffusion coefficient (ADC), thought to be secondary to expansion of the extracellular space in the inflammatory exudate and a consequent increase in the freedom of water diffusion. ADC measurements have been shown to be responsive to tumour necrosis factor inhibitor (TNFi) therapy, and could be useful for monitoring the severity of inflammation in the clinic. 8 The monoexponential ADC model 1 is the simplest and most commonly used model for diffusion attenuation, and provides a useful summary metric of tissue characteristics. However, a disadvantage of this model is that it fails to capture the complexity of biologic tissues, which contain water molecules in different environments with different degrees of diffusion hindrance and/or restriction. This has led to investigation of other diffusion models that aim to capture more detail about these compartments. For example, the intravoxel incoherent motion (IVIM) model proposed by Le Bihan et al was designed to capture translational movements occurring within voxels (such as blood flow in capillaries-perfusion) and can potentially enable separate measurement of diffusion and perfusion. [12][13][14] This model has recently been investigated in spondyloarthritis, and initial results suggest that it can differentiate between patients with differing disease activities. 15 Diffusion kurtosis imaging (DKI)-relying on the assessment of kurtosis to capture non-Gaussian diffusion behaviour 16,17 -has also recently been investigated in spondyloarthritis. 18 DKI was described as a general method for capturing deviations of the water displacement probability distribution from the Gaussian form; the parameter K app provides a dimensionless metric of the degree to which the diffusion displacement distribution is non-Gaussian. 16 DKI can therefore be used as a 'model free' approach for capturing a wide range of tissue phenomena, and minimizes assumptions about the underlying biology of the tissue (although an implicit assumption is that higher order terms-in b 3 and above-in the 'true' signal equation are negligible).
Several studies have demonstrated the feasibility of IVIM and kurtosis imaging in spondyloarthritis (with some promising results), 15,18,19 but there has been no rigorous attempt to assess the accuracy with which these models actually explain the observed signal behaviour. Thus, it remains unclear if these models are in fact better suited to imaging the bone marrow than the monoexponential model, and whether they provide any additional biological or clinical information. In this study, we compared these candidate models in terms of information content, using diffusion data acquired both from normal bone marrow and from areas of bone marrow oedema.

| METHODS
This study received ethical approval from the Queen Square Research Ethics Committee, London (Research Ethics Committee Reference 15/LO/1475). All participants gave written informed consent prior to study entry.

| Subjects
Fifty-three subjects aged 12-24 with known or suspected spondyloarthritis were prospectively recruited by rheumatologists working in a specialist tertiary rheumatology clinic, as previously described. 20 From this cohort, 11 patients with definite evidence of active inflammation (ie unequivocal areas of bone marrow oedema) and 17 patients with no evidence of active inflammation were selected for further analysis. The remaining 25 patients were deemed equivocal for the presence of inflammation and thus excluded from further analysis. The decision regarding active inflammation was based on the mean visual inflammation score from two consultant radiologists (MHC and KR), using an established definition of active inflammation, as previously described. 21

| Acquisition
Subjects were scanned on a 1.5 T Siemens Avanto scanner with integrated posterior coils and anterior surface coils using a conventional Stejskal-Tanner sequence with single-shot echo planar imaging readout, with b-values 0, 50, 100, 300 and 600 s/mm 2 , which produced magnitude image data from the geometric mean of the signals in the three gradient directions. Other acquisition parameters included T R /T E 2400/92 ms, six averages (for all b-values), field of view (FOV) 322 × 322 mm 2 , acquisition matrix 192 × 192, slice thickness 5 mm, axial acquisition plane. Fat suppression was performed using spectral attenuated inversion recovery (SPAIR). The phase encode direction was set to anteroposterior. The acquisition was acceleration with generalized autocalibrating partial parallel acquisition (GRAPPA), using an acceleration factor of 2.
All subjects also underwent 'conventional' MRI 22 consisting of T 2 -weighted short inversion time inversion recovery (STIR; typical parameters T E 75 ms, T R 5000 ms, T I 150 ms, FOV 220 × 220 mm 2 , acquisition matrix 192 × 256, slice thickness 5 mm), T 1 -weighted turbo spin echo (TSE; T E 9.2 ms, T R 350 ms, FOV 240 × 240 mm 2 , acquisition matrix 243 × 256, slice thickness 5 mm) and fat-suppressed post-contrast T 1 -weighted TSE sequences (T E 9.2 ms, T R 396 ms, FOV 220 × 220 mm 2 , acquisition matrix 230 × 256, slice thickness 5 mm), all acquired in two planes (axial and coronal to the sacroiliac joint), as described previously. 21 The conventional MRI scan was used to determine whether subjects had evidence of active inflammation.

| Voxel selection
In patients with evidence of active inflammation by conventional MRI criteria, regions of interest (ROIs) were placed on the largest focal area of bone marrow oedema on the vendor-supplied ADC map in each patient, using the STIR image as an anatomical guide (see Figure 1), by a radiology registrar (TJPB) with five years of musculoskeletal MRI experience. The ROIs were then automatically transferred onto each of the corresponding diffusion-weighted images (which had the same matrix size, FOV and anatomical position as the ADC maps) at each b-value. In patients without evidence of active inflammation on conventional MRI, ROIs were placed on the largest possible single area of normal subchondral bone marrow, taking care to avoid artefacts in the images (including areas of distortion and fat ghosts). For each patient, all of the voxels within the selected ROI F I G U R E 1 Placement of ROIs. ROIs were placed on subchondral bone on the vendor-supplied ADC map, using the STIR image as an anatomical guide. In subjects with inflammation, the ROIs were placed on the largest single area of bone marrow oedema. In subjects without inflammation, ROIs were placed on the largest available area of normal bone marrow on a single slice. ROIs were then automatically transferred to the diffusion-weighted images for the various b-values (each of which had the same voxel size, FOV and anatomical position). All voxels from the ROI were included in the dataset for fitting were included in the dataset for fitting (the number of datapoints was the number of voxels in the ROI multiplied by the number of b-values). To enable assessment of interobserver variability, the procedure was repeated for all included subjects by a second radiology registrar (NS) with three years of musculoskeletal MRI experience.

| Model description and fitting
The three models and the corresponding fitted parameters are summarized in Table 1.
Model fitting was performed using nonlinear regression, implemented with the MATLAB fmincon method. The signal intensities at each b-value were first normalized by the b = 0 image, meaning that S 0 was fixed and not fitted in the model. The parameters for the IVIM and kurtosis models were initialized using the estimates from the monoexponential model-specifically D tissue , D * and D app were initialized using the estimate for D from the monoexponential model, whilst f ivim was initialized with 0. This initialization guarantees that the more complex models will never perform worse than the monoexponential model. R 2 measurements were recorded for each fit. Note that the formulation of the IVIM model has differed slightly between previous works; the formulation used here is essentially the same as that used by Le Bihan and Turner in Reference 14 and is also consistent with previous work in spondyloarthritis. 15 However, it should be noted that some works have used D * +D blood (Reference 23 ) or a single D * parameter (Reference 24 ) rather than D * +D tissue (as used in this work) to specify the diffusivity for the IVIM component.
For each fit, the sum of squared errors (SSE) and R 2 (coefficient of determination) were recorded. From the SSE, we calculated the Akaike information criterion (AIC). The AIC provides an estimate of the relative quality of each statistical model in terms of information content: a lower AIC indicates less information loss (or higher information content), whilst a higher AIC indicates more information loss (or lower information content).
AIC was calculated using the expression where p is the number of model parameters and n is the number of measurements with unique device settings (in this case the number of voxels in the ROI multiplied by the number of b-values).
To address potential overfitting based on standard AIC measurements, we calculated the corrected AIC (AICc), using The AICc is essentially the AIC with an extra penalty term for the number of measurements being small compared with the number of parameters.
As n ! ∞, the extra penalty term converges to 0, and thus AICc converges to AIC.

| Effect of heterogeneity within the ROI
To exclude the possibility that heterogeneity in D between voxels in the ROI might give rise to the observed IVIM or kurtosis effects, we performed voxelwise monoexponential fitting within the ROI in two example subjects with pronounced departure from monoexponential decay.
Two experiments were performed. First, monoexponential fitting was performed such that both D and S 0 were estimated for each voxel in the ROI and then 'simulated' signal intensities were generated from the whole ROI based on the average parameters from the fits in each voxel. The T A B L E 1 Summary of models and fitted parameters. Note that the kurtosis and IVIM models were initialized using the monoexponential model estimate, with parameter start points shown in the right-hand column simulated signal intensities were compared with the true signal intensities in each ROI. Second, the IVIM parameter f ivim was estimated in each voxel using an asymptotic fitting procedure, as described by Pekar et al. 25 Specifically, in each voxel we performed nonlinear least squares fitting, with the first two b-values (0 and 50 s/mm 2 ) excluded; f ivim was then calculated from the intercept and S 0 using f ivim = (S 0 − intercept)/S 0 . Histograms for f ivim were generated from all voxels within the ROI and the median value was compared with the value obtained from the whole ROI.

| Statistics
For both normal and inflamed marrow, we recorded the number of cases in which the monoexponential, IVIM and kurtosis models offered the best performance in terms of AICc. Individual parameters were also compared between inflamed and non-inflamed groups using two-tailed unpaired t-tests. Linear regression was performed between D and D app and between D and D tissue , and the obtained slope values were compared with 1 using two-sided t-tests. Interobserver variability was assessed using Bland-Altman 95% limits of agreement.  Figure 3. In normal marrow, information content was highest for IVIM in 7/17 subjects, monoexponential in 6/17 subjects and kurtosis in 4/17 subjects. In inflamed marrow, information content was highest for monoexponential in 7/11 subjects, IVIM in 2/11 subjects and kurtosis in 2/11 subjects.

Examples
Parameter estimates for each of the three models are shown in Figure 4 and Table 2. Estimates of D, D app and D tissue were all significantly higher in the inflamed group than in the non-inflamed group (P < 0.0001 for all three parameters). Estimates of D * and f ivim did not significantly differ between inflamed and non-inflamed groups (P = 0.12 and 0.23 respectively). K app was significantly reduced in the inflamed group (P = 0.015).
Results of the linear regression analysis between D and D app , and between D and D tissue , are shown in Figure 5. D app resulted in higher estimates of tissue diffusivity than ADC (slope = 1.17, P = 0.094), whilst D tissue resulted in lower estimates than D (slope = 0.85, P = 0.029).
Results of the experiments testing the effect of heterogeneity within the ROI are shown in Figure 6. The simulated curves obtained from voxelwise monoexponential fitting (A-F) deviated only marginally from monoexponential behaviour and did not replicate the departure from monoexponential decay seen in the measured data, suggesting that variations in monoexponential fit parameters across the ROI could not account for the observed signal behaviour. Similarly, median values for f ivim from asymptotic fitting (G, H) were very close to those obtained from the dataset from the whole ROI (the values for the two subjects shown were 0.129 and 0.046 for the asymptotic fit compared with 0.133 and 0.043 for the ROI-averaged fit).
Results of the Bland-Altman analysis for interobserver variability are shown in Figure 7 and Table 3. 95% limits of agreement were relatively narrow and of similar widths for D, D app and D tissue . K app also showed relatively narrow limits but the limits were wider for D * and f ivim compared with the range of values obtained.

| DISCUSSION
Recently, there has been increasing interest in the use of DWI to image inflammatory diseases of the skeleton, and in the use of nonmonoexponential models of diffusion attenuation-including IVIM and kurtosis-to capture tissue complexity. However, there have been no   Of the three models investigated, the IVIM model provided the best description of the signal over the 0-600 s/mm 2 range across normal and inflamed bone marrow. There was a particular advantage of the IVIM model in normal marrow, where it was best able to capture a pronounced fast diffusion component observed in a number of cases. However, both IVIM and kurtosis models showed a reduction in information content compared with the monoexponential model in inflamed marrow, suggesting that the signal behaviour becomes closer to monoexponential decay in the presence of inflammation. In the IVIM framework, this monoexponential behaviour may arise because of the observed increase in D tissue , which means that the difference between D tissue and D * is reduced and the diffusivities of the two compartments becomes similar. In the kurtosis framework, this effect is 'captured' by changes in D k and K app .
The fits obtained from normal marrow ( Figure 2) show that the IVIM model is able to capture the fast diffusion component over b-values from 0 to 100 s/mm 2 more accurately than the monoexponential or kurtosis models, probably accounting for its superior performance. Although the kurtosis model can 'approximate' a large IVIM effect through an increase in D app accompanied by an increase in K app , this can produce spurious positive gradients in the fitted curve over the higher b-values (ie, for b-values above a certain value, increasing b-values can lead to increasing signal-see Figure 2B). This positive gradient arises because of the positive term within the exponent in the kurtosis signal equation; here the assumption that the higher order terms-ie b 3 and above-are negligible is invalid. By contrast the IVIM and monoexponential signal equations can never result in a positive gradient.
Interestingly, there was no significant difference in f or D * between inflamed and normal marrow, and the main alteration in the inflamed group was in D tissue . These results might be seen as surprising because areas of bone marrow oedema are known to demonstrate enhancement after administration of gadolinium, 22 and might therefore be expected to show increased perfusion and thus an increased IVIM effect. However, other authors have also failed to detect significant differences in the IVIM fraction or diffusion coefficient in areas of inflamed marrow. 26 This result implies that enhancement with contrast in areas of bone marrow oedema is not due to an increase in perfusion per se (and might instead reflect increased capillary permeability).
Although the kurtosis model has previously been used to detect deviations from monoexponential behaviour at high b-values (greater than those used here), it can also detect deviations from monoexponential decay at lower b-values, particularly when the underlying diffusivity is high (as is the case with inflamed marrow). This point is highlighted in Figure 8, which shows that the departure from monoexponential decay for a given kurtosis is greater when the underlying diffusivity is higher. Other investigators have used DKI in combination with an acquisition using relatively low b-values to capture complex blood flow in the placenta. 27 However, as discussed above, our data suggest that the IVIM model offers better performance in normal and inflamed marrow.
The results of the voxelwise fitting experiments suggest that the deviation from monoexponential behaviour in normal marrow cannot be explained by heterogeneity within the ROI and is observed even with single-voxel fitting, albeit with some uncertainty due to low SNR in normal bone marrow.
Previous studies using non-monoexponential diffusion models in the context of bone marrow imaging have demonstrated feasibility, but there has been little justification for the specific models used, and they have not typically been compared with the simple monoexponential model. For example, Zhao et al demonstrated that the tissue diffusion term estimated by IVIM (referred to as D slow in their work) could differentiate between patients with and without active inflammation, but did not compare D slow estimates with ADC or provide a rationale for the use of F I G U R E 5 Linear regression of D app and D tissue (from the IVIM model) against ADC. D app results in higher estimates of tissue diffusivity than D, whilst D tissue results in lower estimates than D IVIM. 15 Similarly, Wang et al described the use of DKI to evaluate activity in ankylosing spondylitis, but did not compare the ability of this model to describe the data with that of ADC or IVIM. 18 The current study suggests that non-monoexponential models should be used with caution in light of the potential for overfitting. In line with the results of this study, it has recently been shown that a simple monoexponential model can outperform more complex models for the differentiation of disease activity in ankylosing spondylitis. 28 Our work also provides an indication of the likely biases that will be introduced by using alternative models. Figure 5 shows that D tissue results in lower diffusivity estimates than ADC (the ADC measurement effectively overestimates tissue diffusivity by approximately 10-20%), whilst D app results in higher estimates. This result is logical when inspecting the form of the IVIM in kurtosis models: the additional terms in the IVIM signal equation largely account for the early diffusion, whilst in the kurtosis model D app accounts for the early diffusion component and K app introduces a curvature that accounts for the slower decay at b > 100 s/mm 2 .
A limitation of this study is that model comparison is likely to be influenced by the specific range of b-values. For example, it is conceivable that there may be a more significant kurtosis component revealed at higher b-values (we used a maximum of b = 600 s/mm 2 ), although it may be F I G U R E 6 Effect of heterogeneity within the ROI. Data from two healthy subjects (the same as those in Figure 2) with pronounced departure from monoexponential decay are shown, with a column for each subject. A, B, Histograms of D obtained from voxelwise monoexponential fitting. C, D, The true measured signal in the ROI. E, F, 'Simulated' signals across the ROI based on parameter estimates from monoexponential fits from each voxel. G, H, Estimates of f ivim obtained from voxelwise asymptotic fitting; the dotted red lines indicate the median. The y-axes for A-F use a logarithmic scale difficult to distinguish between a true kurtosis and the effect of the noise floor, given the unrestricted diffusion behaviour in inflamed marrow.
However, as discussed above, kurtosis can also manifest at lower b-values when the tissue diffusivity is high (see Figure 8). In this work, the kurtosis model can be seen as a general method for assessing the departure from monoexponential decay-the parameter K app provides information that is complementary to the IVIM parameters.
It should be noted that the specific form of the IVIM model has varied somewhat between publications. 1,12,13,23,29 In particular, more recent publications have suggested a formulation including a known value for D blood within the IVIM term of the signal equation, on the basis that the diffusion coefficient of flowing blood will be at least as high as that of static blood. 23 However, given that the biological basis of this effect is poorly understood in bone marrow, we avoided this assumption in the present work. To improve continuity, we used the same form of the IVIM model as in previous studies in spondyloarthritis. 15 Further studies using a larger number of datapoints might help to refine the specific form of the IVIM model for bone marrow, and to improve our understanding of the biological relevance of the various parameters. Additionally, this study did not address the performance with which the various parameters could distinguish between inflamed and non-inflamed tissue. Performance is likely to be influenced both by the model selected and by the stability of the fit; further work is needed to assess fit stability for the various models in this work.
F I G U R E 7 Bland-Altman limits of agreement for interobserver variability analysis. Plots of difference against mean are shown for each parameter Finally, although active inflammation can be inferred based on the presence of bone marrow oedema, this does not constitute histopathological proof of inflammation and in some cases other diseases can also produce oedema. This limitation is a strong motivation for further developing microstructural imaging approaches for bone marrow.

| CONCLUSION
Of the three models investigated, the IVIM model provided the best description of the signal over the 0-600 s/mm 2 range across normal and inflamed bone marrow. There was a particular advantage of the IVIM model over the monoexponential and kurtosis models in normal marrow, where it was best able to capture the pronounced fast diffusion component observed in a number of cases. However, IVIM and kurtosis effects both became smaller and the signal behaviour became closer to monoexponential in the presence of bone marrow oedema. Our data suggest that an increase in D tissue (in the IVIM framework) might account for the reduced deviation from monoexponential behaviour in oedematous bone.
F I G U R E 8 Plot of kurtosis model functions showing the effect of variations in D k and D tissue on the signal. For a constant K app , the departure from monoexponential decay becomes more pronounced at higher diffusivity values