Dose-dependent degeneration of non-cancerous brain tissue in post-radiotherapy patients: A diffusion tensor imaging study

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).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.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).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.


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][9][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. Neurooncological 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, neuroinflammation, and necrosis among others [18][19][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 . CC-BY 4.0 International license It is made available under a author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
is the (which was not peer-reviewed) The copyright holder for this preprint . https://doi.org/10.1101/19005157 doi: medRxiv preprint Dose-dependent degeneration of non-cancerous brain tissue 4 permanent [22][23][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.

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.
. CC-BY 4.0 International license It is made available under a author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
is the (which was not peer-reviewed) The copyright holder for this preprint . https://doi.org/10.1101/19005157 doi: medRxiv preprint

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 mm 3 voxel size.

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/mm 2 and one b = 0 s/mm 2 .
Other acquisition settings showed large heterogeneity between patients, but within each patient, the images were acquired on the same type of scanner.

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-diffusionweighted 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].
. CC-BY 4.0 International license It is made available under a author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
is the (which was not peer-reviewed) The copyright holder for this preprint . https://doi.org/10.1101/19005157 doi: medRxiv preprint Robust estimation of the diffusion tensor model was performed with REKINDLE [34]. The regionof-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.
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. 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.
. CC-BY 4.0 International license It is made available under a author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
is the (which was not peer-reviewed) The copyright holder for this preprint . https://doi.org/10.1101/19005157 doi: medRxiv preprint Dose-dependent degeneration of non-cancerous brain tissue 8 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.
. CC-BY 4.0 International license It is made available under a author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
is the (which was not peer-reviewed) The copyright holder for this preprint . https://doi.org/10.1101/19005157 doi: medRxiv preprint . CC-BY 4.0 International license It is made available under a author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
is the (which was not peer-reviewed) The copyright holder for this preprint . https://doi.org/10.1101/19005157 doi: medRxiv preprint

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][40][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 p corr < 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.
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.
. CC-BY 4.0 International license It is made available under a author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
is the (which was not peer-reviewed) The copyright holder for this preprint . https://doi.org/10.1101/19005157 doi: medRxiv preprint 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.
. CC-BY 4.0 International license It is made available under a author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
is the (which was not peer-reviewed) The copyright holder for this preprint . https://doi.org/10.1101/19005157 doi: medRxiv preprint

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,  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.
. CC-BY 4.0 International license It is made available under a author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
is the (which was not peer-reviewed) The copyright holder for this preprint . https://doi.org/10.1101/19005157 doi: medRxiv preprint Dose-dependent degeneration of non-cancerous brain tiss 1 Fig. 4 Whole-brain dose-susceptibility map measured with MD, after significance thresholding. T color shading represents the effect size in terms of percentage change in a year per unit dose, wh the background image is the MNI template (radiological view: left on the image is right in the bra and vice versa). issue 13 The while brain . CC-BY 4.0 International license It is made available under a author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
is the (which was not peer-reviewed) The copyright holder for this preprint . https://doi.org/10.1101/19005157 doi: medRxiv preprint . CC-BY 4.0 International license It is made available under a author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
is the (which was not peer-reviewed) The copyright holder for this preprint . https://doi.org/10.1101/19005157 doi: medRxiv preprint

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 tissuedependent 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][47][48][49]. RTrelated 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][51][52][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][57][58][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 . CC-BY 4.0 International license It is made available under a author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
is the (which was not peer-reviewed) The copyright holder for this preprint . https://doi.org/10.1101/19005157 doi: medRxiv preprint 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][67][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 . CC-BY 4.0 International license It is made available under a author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
is the (which was not peer-reviewed) The copyright holder for this preprint . https://doi.org/10.1101/19005157 doi: medRxiv preprint 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][81][82] or spherical deconvolution [83][84][85], should be taken into consideration.

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.

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.

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).
. CC-BY 4.0 International license It is made available under a author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
is the (which was not peer-reviewed) The copyright holder for this preprint . https://doi.org/10.1101/19005157 doi: medRxiv preprint