Abstract
Background and purpose Radiation-induced changes in brain tissue may relate to post-radiotherapy (RT) cognitive decline. Our aim is to investigate changes of the brain microstructural properties after exposure to radiation during clinical protocols of RT using diffusion MRI (dMRI).
Methods and Materials The susceptibility of tissue changes to radiation was investigated in a clinically heterogenic cohort (age, pathology, tumor location, type of surgery) consisting of 121 scans of 18 patients (10 females). The imaging dataset included 18 planning CTs and 103 dMRI scans (range 2-14, median = 6 per patient) assessing pre-operative, post-operative pre-RT and post-RT states. Diffusion tensor imaging (DTI) metrics were estimated from all scans for a region-of-interest based linear relation analysis between mean dose and change in DTI metrics, while partial volume effects were regressed out.
Results The largest regional dose dependency with mean diffusivity appear in the white matter of the frontal pole in the left hemisphere by an increase of 2.61 %/(Gy x year). Full brain-wise, pooled results for white matter show fractional anisotropy to decrease by 0.85 %/(30Gy x year); mean diffusivity increase by 9.17 %/(30Gy x year); axial diffusivity increase by 7.30%/(30Gy x year) and radial diffusivity increases by 10.63%/(30Gy x year).
Conclusions White matter is susceptible to radiation with some regional variability where diffusivity metrics demonstrate the largest relative sensitivity. This suggests that dMRI is a promising tool in assessing microstructural changes after RT, which can help in understanding treatment-induced cognitive decline.
1 Introduction
Radiation therapy (RT) is a common treatment procedure for both primary and metastatic tumors in the brain, often in combination with surgery and chemotherapy. As radiation impact is not selective to tumor cells, the efficacy of RT is hindered by the radiosensitivity of healthy surrounding tissue [1,2]. The effect of radiation on brain tissue is dynamic and involves structures outside the targeted tumor volume, directly or indirectly [3,4]. Despite the substantial advances in RT technology and application [5,6], the regional differences in sensitivity to radiation dosing is not well-documented, which is especially true for white matter (WM) structures. Current clinical protocols include guidelines for maximum dose for brain parenchyma and several organs at risk for which the direct structure-function relationship is well-understood, or are visible on CT scans, e.g.: brain stem, optic chiasm, optic nerves, cochlea and hippocampi [7]. However, regional constraints on WM and other structures are not part of the standard consideration.
Diffusion MRI (dMRI) is the preferred non-invasive imaging technique to study WM structure in the brain [8–10]. The fundamental principle behind dMRI is that the Brownian motion of water molecules is dependent on the surrounding tissue structure [11]. WM is a highly structured tissue type, which introduces diffusion anisotropy along a certain direction and restricts the extent of local diffusion. This is reflected in the diffusion tensor imaging (DTI) metrics, such as fractional anisotropy (FA) and mean, radial, and axial diffusivity (MD, RD, and AD, respectively), making dMRI a quantitative method to examine WM anatomy and its microstructural properties. Neuro-oncological research has placed dMRI as a potential technique for tumor diagnosis [12], surgical planning [13], pre-treatment prediction of tumor response [14], monitoring early efficacy of treatment [15] and early WM damage post-radiation [16], and, more recently, late-delayed effects of RT on WM [17].
Radiation-induced WM damage has been reported to include axonal injury, demyelination, neuro-inflammation, and necrosis among others [18–20]. Importantly, these structural deficits seem to correlate in time with both verbal and non-verbal cognitive impairment, including executive functioning, working memory, visuospatial processing, and decision-making [21]. In this context, radiation-induced cognitive impairment has been divided into multiple phases post-RT: acute (<2 weeks), early-delayed (2 weeks to 3 months) and late-delayed (6+ months). Notably, while acute and early-delayed symptoms and damage are mostly temporary, late-delayed damage is considered permanent [22–24]. The introduced cognitive impairment currently occurs in 50-90% of survivors [25,26], and this population is increasing with RT advancement [27,28]. The progressive decline affects the physical and mental health of long-term survivors and impairs the patients’ quality of life (QoL) [22,23]. Therefore, a better understanding of the mechanism of post-irradiation cognitive decline is necessary. High RT doses (>60Gy) are well-documented to cause permanent damage, but emerging evidence shows that even at lower doses (<20Gy) late-delayed damage can still occur [29]. Information on regional WM sensitivity to RT dosage in the long term is crucial for better RT planning and ensuring optimal QoL after treatment.
In this work, we aim to quantify the long-term or late-delayed effects of RT on brain tissue microstructures using dMRI. Specifically, we are interested in how the susceptibility of different structures to radiation varies across regions and dose levels.
2 Materials and Methods
2.1 Cohort
We retrospectively identified patients who were treated with RT at the Department of Radiation Oncology MD Anderson Cancer Center, University of Texas, Houston, Texas, USA. Patients were eligible for inclusion upon met the following criteria: inclusion of a pre-RT MRI scan, which is defined as the baseline (BL) of a scan within one month of the RT start, and follow-up (FU) scans between 3 and 24 months to exclude the acute-early and very late effects. The selection resulted in 18 patients (10 females). All patients received intensity-modulated RT (IMRT); most of them were treated to 60 Gy in 30 fractions, others were recalculated to a total equivalent dose in standard 2 Gy fractionation (EQD2) using biologically equivalent dose principles using the linear quadratic model with an α/β ratio of 2 Gy [30]. Detailed information of treatment and demographic factors are shown in Table 1. This study was approved by the institutional review board and informed consent was obtained from all subjects.
Demographics and tumor characteristics of the patients included in the cohort.
2.2 Image acquisition
2.2.1 CT acquisition
The planning CT scans were acquired on a Philips Brilliance Big Bore scanner, with a tube potential of 120 kVp, matrix size of 512 × 512 × 87 and 0.98 × 0.98 × 3.0 mm3 voxel size.
2.2.2 MRI acquisition
MR images were acquired at multiple time points for each patient on various GE Medical Systems MRI scanners. The dataset included 103 dMRI scans (range 2-14, median = 6 per patient) acquired using 27 gradient directions with a diffusion weighting of b = 1200 s/mm2 and one b = 0 s/mm2. Other acquisition settings showed large heterogeneity between patients, but within each patient, the images were acquired on the same type of scanner. At 1.5T, examinations were performed with Signa Excite (in 6% of the 103 scans) and Signa HDxt (51%) scanners. For all 1.5T scans, the dMRI acquisition details were: TR = 10 s, TE = 102 ms, voxel size: 0.86 mm × 0.86 mm × 6.5 mm in 58% of all scans. At 3T, the following scanners were employed: Signa Excite <1%, Signa HDx, 2%, Signa HDxt scanner 39%. For all 3T scans, the dMRI acquisition details were: TR = 8.75 s, TE = 91 ms, voxel size: 0.86 mm × 0.86 mm × 3.5 mm in 42% of all scans.
2.3 Image processing
All non-dMRI data were processed with in-house algorithms developed in Matlab (Mathworks, Natick, Massachusetts), whereas all dMRI data were processed with ExploreDTI [31].
First, the CT images were cropped to increase registration performance and to reduce computational load in later stages. This step resulted in a bounding box around the skull, excluding the neck and the shoulders of the patients. The same cropping was applied for the dose, organs-at-risk (OaR) and planning, clinical, and gross target volume (TV) maps per patient.
dMRIs were corrected for subject motion and eddy current distortions using the non-diffusion-weighted image as the reference, whereupon all images were rigidly registered to the cropped CT image [32]. In order to reduce the blurring effect of multiple registration procedures, the two transformations were concatenated into a single resampling step. During motion correction and subsequently during registration to the CT, the DW gradients were adjusted with the rotational element of the registrations [33].
Robust estimation of the diffusion tensor model was performed with REKINDLE [34]. The region-of-interest (ROI) definition for the analysis was performed by non-linearly registering the ‘wmparc’ brain template and the associated atlas to every individual DWI scan [35], which were already motion-corrected and co-registered to the CT scan. The mentioned atlas contains 179 labels covering the whole brain and is freely available in Freesurfer [36]. A graphical overview of the automated image processing pipeline is shown in Fig. 1.
Graphical overview of the image processing pipeline. First, the raw diffusion weighted images (DWIs) registered to the very first non-DW image. Next, all images were further registered to the CT scan followed by a robust diffusion tensor estimation, which leads to the calculation of diffusion tensor imaging scalars as fractional anisotropy (FA), mean-radial-axial diffusivity (MD, RD, and AD, respectively). A standard template from the Montreal Neurological Institute (MNI) and associated labels from Freesurfer were non-linearly registered to the already CT-coregistered DWIs. The registered labels allow an atlas-defined ROI based analyses of DTI metrics.
The processing pipeline generates DTI metrics based on atlas-defined ROIs in the native CT space, yielding high spatial alignment between anatomical and DTI data for each timepoint. Finally, equivalent dose in 2 Gy per fraction (EQD2) dose values and DTI metrics from multiple time points were accessible per ROI per patient. We computed the following DTI metrics per ROI: average of FA, MD, AD, RD; volume in mm3 and mean of the EQD2 dose. Examples of the CT-dMRI alignment are shown in Figs. 2/a-c, while Figs. 2/d-f show changes of the tissue and the tumor cavity after the start of RT with an exemplar 3D rendering of the inferior temporal white matter in the non-affected hemisphere.
Graphical overview image processing results. Panels a-b-c shows the CT and dMRI alignment via fused CT - direction encoded color (DEC) maps in representative patients, where the background image is the planning CT and the colors show the locally dominant direction of diffusion via the DTI model. The intensity of the color was weighted with fractional anisotropy (FA), which suppresses the color in CSF, because FA is close to 0 in homogeneous media. White arrows pinpoint critical regions of alignment at (a) edge of the ventricles and gyrus, (b) interface of the brain-skull, and at (c) ventricle – white matter interface. Red arrows show the tumor cavity, where the absence of strong anisotropic diffusion resulted in no color coding. In the bottom row: (d)-(e)-(f) show the changes of the tumor cavity of patient (b) with mean diffusivity (MD) as the background image. An example of an atlas region is shown in red: the 3D rendered region of the inferior temporal white matter in the non-affected hemisphere.
To avoid data contamination from tumor cavities, remaining cancerous tissue, or surgical scars, all ROIs that overlapped the GTV were excluded from the analysis. Also, ROIs covering the ventricles and cerebrospinal fluid CSF were ignored. Moreover, regions with the atlas label ‘Unsegmented White Matter’ were excluded due to its undefined nature. In total, 16 of the original 179 labels were ignored and a total of 6.85% of the 18*163 regions were not included for further analysis.
2.4 Longitudinal analysis
Analyzing data from unequally time-spaced acquisitions can be challenging, given that standard statistical models of repeated measures are not applicable, like paired t-tests or ANOVAs. One solution is to perform regression analysis for every ROI in every patient with the following details: the timing of the follow-up scans is the independent variable, while the mean FA, MD, AD and RD per ROI are the dependent variables. The pre-RT MRI scan is defined as the baseline (BL) with day = 0. This calculation yields in the rate of change of the DTI metrics. In another pipeline, we included the volume of the ROIs as a covariate of no-interest as Vos et al. showed that partial voluming can act as a hidden covariate in DTI analyses [37]. Calculations described above were performed with FSL utility tool fsl_glm [38].
Next, the estimated coefficients from all four dMRI metrics were correlated with the mean EQD2 dose per ROI using a permutation test with 5000 iterations performed with the permutation analysis of linear models (PALM) toolbox version alpha104, a Matlab based open-source software package [39–41]. We used nonparametric permutations as they proved efficient in eliminating false positive results when compared with parametric methods [42]. Significance of a correlation was determined at pcorr < 0.05 using family-wise error rate (FWER) adjustment to correct for multiple comparisons. Tail approximation was used for faster calculations [43]. Age at the time of the diagnosis and sex of the patients were included as nuisance regressors. An effect size measure, the Cohen’s D, is reported to grade the practical significance of the results [44]. A graphical overview of the longitudinal analysis is shown in Fig. 3.
Graphical overview of the longitudinal analysis. First the motion-corrected dMRI metrics are grouped by atlas labels as well as doses. Regions that overlapped with the GTV were excluded from the analysis. The rate of metric change after RT was calculated for every region in every patient. The individual rates were then pooled to calculate the dose dependency of the metric changes. Repeating the procedure for all regions and metrics leads to whole brain dose-susceptibility maps.
Additionally, variance groups (VGs) were defined based on the number of scans per patient, since permutations are only valid within the same VG. Defining the VGs is needed because the variance of the changes in the metrics is not equal in all subjects and depends on the number of scans used. This introduces heteroscedasticity in the metrics, which we accounted for by limiting permutations via VGs.
3 Results
Fig. 4 shows the whole-brain dose susceptibility map as measured with MD for those regions which were statistically significant using a critical p-value = 0.05 as significance threshold. In the majority of brain regions, MD increases with increasing dose. Also, some regional differences appear, with the largest dose dependencies in the anterior part of the Corpus Callosum by 1.77 % / (Gy x year); in the caudal anterior cingulate cortex (ACC) in the right hemisphere (also known as Brodmann areas 24, 32 and 33) by 1.78 % / (Gy x year); and in the white matter of the frontal pole (also known as Brodmann area 10) in the left hemisphere by 2.61 % / (Gy x year). Supplementary Figs. 1, 2 and 3 show the spatial distribution of dose susceptibility maps of FA, AD, and RD, respectively. Supplementary Fig. 4 shows the overall mean dose distribution, where mostly the deep grey matter (GM) regions received the highest dose on average. Supplementary Fig. 5, 6, 7, 8 and 9 shows the coefficients of change along with the 95% confidence interval for regions in the left hemisphere of the cortex; right hemisphere of the cortex; left hemisphere of the white matter; right hemisphere of the white matter; and the deep GM regions, along with the brainstem, segments of the corpus callosum and the cerebellum; respectively.
Whole-brain dose-susceptibility map measured with MD, after significance thresholding. The color shading represents the effect size in terms of percentage change in a year per unit dose, while the background image is the MNI template (radiological view: left on the image is right in the brain and vice versa).
Table 2. shows pooled results for 30 Gy per cortical GM, deep GM, and WM; with and without considering region volume as a covariate of no-interest. After accounting for volume, all metrics increased, which resulted in elevated diffusivity susceptibilities and a diminished FA-dose relation. Supplementary Table 1 shows the cohort-based mean values cortical GM, deep GM, and WM. In WM, FA shows hardly any response to the applied dose, whereas the increase in diffusivities shows larger susceptibility to dose than FA in all tissue types. The increase of RD in WM is nearly 50% greater than AD, meaning that the diffusivity increase perpendicular to the dominant fiber direction is much larger than the increase along the dominant direction for a given dose. In cortical GM, RD is nearly 15% higher than AD, whereas in deep GM, RD is 10% lower than AD.
Pooled percentage rates of tissue degradation with different DTI metrics.
4 Discussion
In this work, we presented a systematic ROI-based analysis of microstructural properties of the brain in patients after cranial RT. The results show that the changes estimated with DTI are tissue-dependent and are much more pronounced in the diffusivity metrics (MD, AD and RD) than in FA. Studies on tissue-dose interaction, as the current one, aid in understanding radiation-induced cognitive decline, which may be the consequence of, but not restricted to, demyelination and loss of axonal integrity with increasing radiation doses. Research in neuroscience, especially on aging [47,48] and cognitive impairment [49,50] have already established a link between WM properties as determined with DTI and mental health of the participant [45]. Information about regional sensitivity to RT is the first step towards new cognition-sparing radiation treatment planning resulting in lower adverse side effects.
Volume changes in brain structures have been well-studied because of their potential link to cognitive functions, such as the volume-memory relationship in the hippocampus [46–49]. RT-related changes in hippocampal volume seem to be consistent, show dose-dependency, and correlate with follow-up memory assessments, thereby promoting the hippocampal volume as an important biomarker and an organ to protect from radiation [23,45,50–53]. In summary, region volume is affected by RT and for this reason we incorporated volume as a nuisance during statistical inference. This led to increase of all metric change rates and resulted in having the largest effect on FA. It confirms the findings of Vos et al [54] that the volume of a bundle or region can act as a hidden covariate in quantitative analyses. Furthermore, the effects of partial voluming on FA are remarkable; neglecting it would lead to a conceptually different outcome. The partial volume effects (PVE) are reduced in the DTI metric changes, therefore the results represent only the microstructural-related changes as estimated with DTI. Studies neglecting PVEs are prone to unreliable analysis. To the best of our knowledge, previous works on regional tissue sensitivity to RT discuss volumes only within ROIs receiving certain dose levels or aim to avoid PVE by restricting the masking to connected voxels that only include WM structure by excluding voxels with FA lower than 0.2 [17,29,55].
Previous studies investigating tissue properties with dMRI and radiation are in line with our work: an overall decrease of FA and an increase of the various diffusivities [56–59]. The exact underlying mechanisms are ambiguous, since changes of FA are not specific to one well-defined procedure, like demyelination, scarification [60] and axonal loss [61]. However, some relations are possible to recover. The appearance of edema reduces FA while it increases MD [62,63], which is a result of increased presence of free water [64]. With the measurements of AD and RD, it is possible to look into the change of diffusivities in depth. RD is 50% more sensitive to radiation than AD in WM, meaning that diffusion becomes increasingly isotropic with a progressive dose. These finding are consistent with previous studies, confirming that the increase in RD overcomes the increase in AD [19,55]. Similarly, an increase in RD has been suggested to relate to demyelination in human [65] and animal [20] RT studies. It is important to note that, while cortical and deep GM show large relative changes in FA, these changes are negligible in absolute terms since FA in GM is close to 0 due to the absence of anisotropic media at the typical acquisition resolution.
Regional differences appear, but it remains doubtful whether these are genuine or artefactual as previous examples pooled the bilateral regions into a single one, forcing symmetrical results [55]. Clearly, the environment for testing tissue susceptibility to dose is not equal for all regions as dose distribution is non-homogeneous. Much larger clinical cohorts are necessary in the future to selectively balance out this tendency, which originates from the spatial heterogeneity of brain tumors [66–68].
Tissue changes after radiotherapy can be also addressed via analyzing T1-weighted (T1w) images. In this work, however, we excluded the T1w scans from the analysis because of highly anisotropic voxel sizes. The larger slice thickness will cause significant PVEs that will impede a meaningful morphometric analysis. New guidelines [69] have been proposed to set up robust acquisition protocols and processing pipelines to tackle volumetric [70] and surface features (e.g., cortical thickness [71]). T1-based brain analysis has a long history in neuroscience, revealing similar changes to RT as in healthy aging, but in a period of decades as opposed to our maximum 24 months long follow-up period [72]. Accelerated aging as a consequence of RT has been also described by non-MRI studies investigating physiologic frailty [73] and accelerated neurocognitive function decline [74].
For the modeling of tissue changes we assumed a linear relationship between the tissue characteristics and time passed since RT. We further assumed that tissue change relates linearly to the applied dose. While this framework is widely applied by the RT-community [17,29,55,75], these assumptions may only approximate the true underlying relations. Further biophysical modeling is necessary to establish a more accurate representation of the tissue-dose relation and, in turn, to provide a better understanding of the underlying tissue breakdown mechanisms [76].
In this study, we considered the relation of DTI metrics and radiation, which - regardless of the significance of the association - does not necessarily translate into clinically significant effects. Furthermore, it provides no information of the disturbed functionality of the regions, as one may consider that a certain change in, for instance, the hippocampi can have different consequences than the same change in the temporal cortex.
The current study uses DTI as a means of diffusion modeling. However, it has been implied that DTI has some critical limitations, like the inability to resolve crossing fibers [77,78]. Recent developments in diffusion MRI suggest [79] that more advanced models than DTI, like diffusional kurtosis imaging [80–82] or spherical deconvolution [83–85], should be taken into consideration.
5 Conclusion
In this work, we investigated radiation-induced changes in the brain measured with DTI in a clinical population. Our findings suggest that radiation-caused damage of brain tissue follows a similar process as aging, but in a much faster fashion. Future studies should incorporate auxiliary modalities as functional MRI or perfusion mapping [86] as well as behavioral and cognitive tests for a more detailed investigations to reveal the causal relationship between changed tissue characteristics and changes in mental health.
Data Availability
The data generated during the study is available from the corresponding author on reasonable request by a qualified individual or third-party, subject to authorization of relevant regulatory bodies.
6 Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
7 Funding
The research of S. D., H. Y. M. and A.L. is supported by VIDI Grant 639.072.411 from the Netherlands Organization for Scientific Research (NWO).
8 References
- [1].↵
- [2].↵
- [3].↵
- [4].↵
- [5].↵
- [6].↵
- [7].↵
- [8].↵
- [9].
- [10].↵
- [11].↵
- [12].↵
- [13].↵
- [14].↵
- [15].↵
- [16].↵
- [17].↵
- [18].↵
- [19].↵
- [20].↵
- [21].↵
- [22].↵
- [23].↵
- [24].↵
- [25].↵
- [26].↵
- [27].↵
- [28].↵
- [29].↵
- [30].↵
- [31].↵
- [32].↵
- [33].↵
- [34].↵
- [35].↵
- [36].↵
- [37].↵
- [38].↵
- [39].↵
- [40].
- [41].↵
- [42].↵
- [43].↵
- [44].↵
- [45].↵
- [46].↵
- [47].↵
- [48].↵
- [49].↵
- [50].↵
- [51].
- [52].
- [53].↵
- [54].↵
- [55].↵
- [56].↵
- [57].
- [58].
- [59].↵
- [60].↵
- [61].↵
- [62].↵
- [63].↵
- [64].↵
- [65].↵
- [66].↵
- [67].
- [68].↵
- [69].↵
- [70].↵
- [71].↵
- [72].↵
- [73].↵
- [74].↵
- [75].↵
- [76].↵
- [77].↵
- [78].↵
- [79].↵
- [80].↵
- [81].
- [82].↵
- [83].↵
- [84].
- [85].↵
- [86].↵