CrOssing fiber Modeling in the Peritumoral Area using dMRI (COMPARI)

Visualization of fiber tracts around the tumor is critical for neurosurgical planning and preservation of crucial structural connectivity during tumor resection. Biophysical modeling approaches estimate fiber tract orientations from differential water diffusivity information of diffusion MRI. However, the presence of edema and tumor infiltration presents a challenge to visualize crossing fiber tracts in the peritumoral region. Previous approaches proposed free water modeling to compensate for the effect of water diffusivity in edema, but those methods were limited in estimating complex crossing fiber tracts. We propose a new cascaded multi-compartment model to estimate tissue microstructure in the presence of edema and pathological contaminants in the area surrounding brain tumors. In our model (COMPARI), the isotropic components of diffusion signal, including free water and hindered water, were eliminated, and the fiber orientation distribution (FOD) of the remaining signal was estimated. In simulated data, COMPARI accurately recovered fiber orientations in the presence of extracellular water. In a dataset of 23 patients with highly edematous brain tumors, the amplitudes of FOD and anisotropic index distribution within the peritumoral region were higher with COMPARI than with a recently proposed multi-compartment constrained deconvolution model. In a selected patient with metastatic brain tumor, we demonstrated COMPARI’s ability to effectively model and eliminate water from the peritumoral region. The white matter bundles reconstructed with our model were qualitatively improved compared to those of other models, and allowed the identification of crossing fibers. In conclusion, the removal of isotropic components as proposed with COMPARI improved the bio-physical modeling of dMRI in edema, thus providing information on crossing fibers, thereby enabling improved tractography in a highly edematous brain tumor. This model may improve surgical planning tools to help achieve maximal safe resection of brain tumors.

Introduction 1 subsequent tractography reflects the underlying anatomy, facilitating surgical decisions. 22 Thus, the overarching goal of this paper is to design a multi-compartment model that 23 will enable tracking through crossing fiber regions in the presence of peritumoral edema. 24 In our previous work, FERNET [7], we modeled the peritumoral edema using a 25 bi-tensor model with free-water and anisotropic compartments, fitted to single-shell 26 data. Although FERNET enhanced tractography in the peritumoral area, as a 27 single-tensor model, it was unable to model crossing fibers. There are dMRI models for 28 tractography that fit a fixed number of fiber populations per voxel [8,9] which have not 29 been tested in the peritumoral region. Recently, dMRI models that facilitate tracking 30 through crossing regions and allow for any number of fibers in each voxel using 31 constrained spherical deconvolution (CSD) have been developed, but these also have not 32 been tested in the peritumoral region. In these methods, including multi-shell 33 multi-tissue CSD (MSMT-CSD) [10]], which models "CSF-like" and "GM-like" signals 34 in addition to the WM fiber orientation distribution (FOD), the response function of 35 WM is estimated from the highest fractional anisotropy (FA) voxels in the brain, 36 usually in corpus callosum, and is fixed for the entire image. The assumption that one 37 response function can be applied to all WM, however, may not hold in pathological WM 38 tissue. In [11,12], special dMRI acquisitions were designed to enable modeling of WM 39 microstructure with unprecedented detail, but their model did not estimate fiber 40 orientations and the acquisition was not clinically feasible, and therefore the method 41 cannot be used for tractography. 42 Consequently, to address these shortcomings in current literature, and in 43 acknowledgement of the need to model the peritumoral tissue for better tractography, 44 we have developed a new tissue modeling approach for clinically feasible product components, to better isolate the underlying anisotropic microstructure of peritumoral 50 tissue. We first evaluate our proposed method on simulated data. We then demonstrate 51 enhanced modeling of crossing fibers in the peritumoral area and improved tractography 52 in brain tumor patients.

54
Multi-compartment model (COMPARI -CrOssing fiber 55 Modeling in the Peritumoral Area using clinically feasible dMRI) 56 We propose a cascaded model, comprising an initial isotropic fitting, followed by 57 multi-compartment modeling of the tissue into isotropic and anisotropic compartments. 58 With the fractions of the compartments fixed, we finally estimate an FOD with 59 multi-compartment spherical harmonics fitting [13]. 60 We fit a multi-compartment (MC) consisting of two bundle compartments: an 61 isotropic bundle and an anisotropic bundle. The isotropic bundle consists of two ball 62 compartments corresponding to free water and hindered water. The anisotropic bundle 63 consists of stick and zeppelin compartments, as described in [14], corresponding to 64 restricted intra-axonal and hindered extra-axonal diffusivities, respectively.

65
Mathematically, the four compartments are modeled by: anisotropic-bundle f r E(b, n n n, µ µ µ, λ ∥ ) restricted intra-axonal attenuation of free water with λ csf = 3 × 10 −9 m 2 /s, and E(b, λ iso ) = exp(−bλ iso ) , 68 represents isotropic signal attenuation of hindered water with magnitude λ iso . λ iso is 69 estimated by fitting a model with a single isotropic ball with diffusivity between intra-axonal and hindered extra-axonal diffusion, respectively.

75
The stick is an anisotropic cylinder model with negligible diameter represented by 76 E(b, n n n, µ µ µ, λ ∥ ) = exp(−bλ ∥ (n n n T µ µ µ) 2 ), where λ ∥ is the magnitude of Gaussian diffusion 77 along µ µ µ orientation, and n n n represents the diffusion gradient vector. In our model, the λ ∥ 78 is fixed and equal in both stick and zeppelin, and is estimated by the voxel selection 79 method described in [15], which automatically chooses high FA voxels, forming a 80 reference for WM. In contrast, the hindered extra-axonal diffusion is modeled as a 81 symmetric Gaussian distribution (zeppelin) represented by aligns the zeppelin to the stick orientation µ, λ ∥ >= λ ⊥ , and D D D h diag = diag(λ ∥ , λ ⊥ , λ ⊥ ). 84 We impose a tortuosity constraint between the diffusivity values of the zeppelin and f r , 85 the fraction of the restricted intra-axonal (stick) portion of the bundle, with [16]. Finally, we apply a condition whereby if λ iso <= λ ∥ , the fraction 87 f csf is set to 1 and the isotropic compartment simplifies to one ball representing free 88 water.

89
Next, we estimate coefficients c c c of a spherical harmonic function Y m l (θ, ϕ of degree l 90 and order m, where θ ∈ [0, π] is the inclination angle and ϕ ∈ [0, 2π] is the azimuthal 91 angle in polar coordinates, utilizing CSD as in the multi-compartment spherical 92 harmonics models (MC-SH) described in [10]. The Fiber Orientation Distribution where f * indicates that the fitted parameters of Equation 1 are fixed, and * S 2 describes 97 convolution on the sphere. The number of acquired gradient directions defines the 98 maximum harmonic order l max . However, estimating higher l max is possible by knowing 99 prior information of required spherical harmonic coefficients of l max , assuming the 100 non-existing directions to be zero. This is the super-CSD approach as defined in [17], 101 and we use the regularization parameter (λ = 1) and threshold orientation density (τ = 102 ten percent of the mean initial FOD amplitude) which were shown to be optimal [18].

103
The schema of our proposed method is shown in Fig. 1 104

105
To evaluate the accuracy of our proposed method using data with known ground truth 106 for GM, WM, CSF and edema, we simulated a diffusion weighted imaging (DWI) signal 107 using Phantomas [19] for a range of conditions and parameters, including volume 108 fractions, SNR, fiber separation angles, and emulating single-shell and multi-shell data. 109 Three protocols of single-and multi-shell data were simulated ( . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted May 9, 2023. ; with the pipeline described in [24] which included DeepMedic [25] segmentation of the 144 tumor and peritumoral area. The T1 image and segmentation were finally registered to 145 the multi-shell DTI acquisition with ANTs software [26].

146
To compare FODs in our model, we fitted COMPARI and MSMT-CSD in each 147 . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted May 9, 2023. ; patient. We calculated the maximum peak amplitude of the FOD in every voxel in the 148 peritumoral region. To illustrate the degree of anisotropy in our model, we used the 149 anisotropic index (AI), an alternative to FA for FODs, as described in [27]. In addition, 150 we selected one patient with metastatic tumor and edema affecting the centrum 151 semiovale, a region of the brain known to have WM fiber crossings, to visualize number 152 of peaks and FODs in the peritumoral region, and to visually compare results of 153 COMPARI in edema to the contralateral WM. In the same patient, we fitted COMPARI 154 and MSMT-CSD [10] as well as NODDI [28], and the multi-shell two-compartment 155 model for free-water elimination from Hoy et al [29]. We visualized orientations of FOD 156 peaks or tensor orientations between these models in the peritumoral region. Finally, we 157 performed whole brain probabilistic tractography with both models, using the IFOD2 158 algorithm [30] and ] and the SIFT model to determine seed points dynamically [31] in 159 our selected patient. The 'cutoff' parameter for COMPARI was 0.2, while it was the parameter estimates. Then, the CSD estimation was employed based on the approach 173 proposed by [10]. The λ ∥ initiated using response function described in [10].

179
We assessed our algorithm in a simulation setup with different parameters. Then, the 180 tuned parameters were utilized in in vivo experiments.

182
We first assessed fiber crossings at angles of 90 • , 60 • , and 45 • in the context of simulated 183 peritumoral edema. The peak separation rate and error of FOD separation angles were 184 plotted for P1 and P2 (single-shell) protocols, with SNR=30 and l max = 8 using CSD 185 method and isotropic volume fractions vf20%, vf40%, vf60%, and vf80% (Fig.3). 186 We assessed the peak separation rate and error in crossing fiber angle with 187 l max = 8, 10, 12, 14 and 16 at SNR=30 in P3 protocol (Fig.4). We observed a significant 188 improvement in resolving smaller separation angles with l max > 8 compared to l max = 8. 189 In addition, the peak separation rate and error in crossing angle were improved for 60 • . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted May 9, 2023. ;  Furthermore, we separately assessed P3 protocol (multi-shell), with SNR= 20,30,40 196 and l max = 12 (Fig. 5). As we expected, increasing Rician noise reduced the peak  (Fig.6). Median amplitudes as well as median AI were higher in COMPARI 207 than in MSMT-CSD in all patients of the clinical dataset, with average FOD amplitude 208 of 0.10 in MSMT-CSD compared to 0.54 in COMPARI.

210
COMPARI was applied in the peritumoral region as well as in the contralateral healthy 211 homologous WM in a selected patient with brain metastasis (Fig.7). The number of  . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted May 9, 2023. ; Fig 5. Simulation results of the effect of SNR and edema on the ability to resolve crossing fibers. For all experiments, the multi-shell protocol (P3) was used and l max = 12. Higher SNR is associated with improved peak separation rate and error. Top: the peak separation rate of recovering peaks for ground truth crossing angles of (45 • ,60 • ,90 • ) with various isotropic volume fractions denoting different extent of peritumoral edema at three levels of Rician noise (SNR= 20,30,40). Bottom: the distribution of errors in the crossing fiber angle are shown as boxplots. . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted May 9, 2023. ; Visualization of multi-shell models 217 We visually compared orientation maps of MSMT-CSD with COMPARI in edema 218 region. Additionally, we applied the NODDI [28] and HOY et al [29] approaches to our 219 multi-shell data to compare the fitting of the edema region using multi-compartment 220 and free-water approaches respectively (Fig.8). In NODDI (Fig.8a) and in the Hoy  Tractography comparison between FOD models 227 COMPARI-based tractography was compared to that of MSMT-CSD (Fig.9).

228
Streamlines of the CC partially crossed the CST using the MSMT-CSD (Fig.9a). AF 229 and SLF III fibers located in the edema could not be visualized. In contrast, the 230 crossing of the CST and CC fibers could be visualized with COMPARI, with a better 231 qualitative definition of the CST (Fig.9b). Additionally, the AF and SLF III were 232 visualized crossing each other in the peritumoral region with COMPARI.

234
In this study, we introduced a multi-compartment modeling approach (COMPARI) to 235 resolve crossing fibers in the peritumoral edema of brain tumor patients. Simulated and 236 in vivo experiments show that crossing fibers were enhanced in edematous regions.

237
In simulation results (Fig.3-5), we tested the effects of b-value, number of DW 238 directions, SNR, volume fraction and spherical harmonic order parameter (l max ) on the 239 ability of our model to detect crossing fibers. As shown in Fig. 3, our model could . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted May 9, 2023. ; . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted May 9, 2023. ; including a low b-value to estimate the isotropic diffusivity and a high b-value with a 244 large number of DW directions to resolve angular separation, is a prerequisite for our 245 model. In results on simulated multi-shell data using CSD fitting (Fig.5), as expected, 246 increased SNR was associated with higher peak separation rate in resolving crossing 247 fibers, although this effect was small. We could not resolve fibers crossing at 45 • , in line 248 with [37]. However, with the super-CSD approach (Fig.4), we could resolve crossings as 249 low as 45 • when volume fraction was as low as 20%. In higher volume fractions, as 250 expected, the ability to resolve crossing fibers decreased, but our model could resolve 251 crossings of 90 • with nearly 100% peak separation rate and a median of approximately 252 5 • of error. It could also resolve 60 • crossings, albeit with a lower peak separation rate. 253 For all values of l max tested, the super-CSD fit (l max > 8) outperformed the CSD fit, 254 with not much difference between values of l max . Additional simulated experiments (S1 255  Table) showed that increasing b-values and number of gradient directions improved  In human data, there is no ground truth to assess the quality of COMPARI fit.

259
Metastatic brain tumors are known to have the lowest level of tumor infiltration in the 260 peritumoral area [38] with a high degree of edema. Accordingly, we compared fiber 261 crossings in the edematous peritumoral area near the centrum semiovale of one patient 262 with brain metastasis to the contralateral homologous area. We found that COMPARI 263 model was able to resolve fiber crossings similarly in both hemispheres (Fig. 7). Upon 264 comparison with other multi-shell models (Fig.9), we found that crossing fibers in 265 edema with non-negligible amplitude were only identified using COMPARI. Since the 266 amplitude cutoff is a parameter essential for tractography, corresponding results of 267 COMPARI were considerably improved over MSMT-CSD (Fig. 9). To ensure that this 268 increase in amplitude is a generalizable result, we repeated this in 22 other patients, and 269 found that COMPARI FODs consistently had higher amplitude and higher AI than 270 those of MSMT-CSD (Fig. 6).

271
The diffusion kurtosis effect [39] is known to cause inaccurate estimation of 272 diffusivity using Gaussian diffusion models at b-values higher than b=1500 s/mm 2 . As 273 May 1, 2023 12/19 . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted May 9, 2023. ; a result, our model performed poorly on a human data set with b=3000 s/mm 2 and 128 274 directions (S21 Fig), since it relies on an estimate of isotropic diffusivity. This effect was 275 not seen in simulated data due to a limitation of the simulation model. Thus, from our 276 experiments on simulated and human data, we concluded that best results with 277 COMPARI would be achieved using a multi-shell dataset with a low b-value shell to 278 capture isotropic diffusivity, and one or more high-b shells to resolve fiber directions, fit 279 with l max = 12 with the super-CSD approach, if possible.  There are several limitations of this study. First, our model is based on several 298 constraints, including the tortuosity constraint and equality of the λ ∥ value for stick and 299 zeppelin compartments, that were made with the goal of modeling crossing fibers in 300 clinically feasible data. Therefore, derived values such as volume fractions should not be 301 ascribed an explicit biological meaning, although they may be useful as clinical markers. 302 A second potential limitation of our method is that it could fail in the grey matter 303 (GM) if isotropic diffusivity is higher than λ ∥ . In this case, FODs fitted to the GM had 304 an excessive number of peaks (Fig. 7), which could result in false positives in 305 tractography connecting GM regions. As such, care should be taken to fit COMPARI 306 only to the WM, or to restrict tractography to the WM and peritumoral area, such as 307 with anatomically constrained tractography [40,41]. Third, while simulated results 308 suggest that best results with COMPARI will be achieved with multi-shell datasets,

309
COMPARI has not yet been tested in in vivo multi-shell datasets with a variety of 310 b-values and direction sets.

311
In clinical practice, tractography could be useful for clinicians treating patients with 312 brain tumors, especially in the context of surgical planning. Indeed, tumor resection is 313 the first treatment for most tumor types [42], and the quality of the resection (i.e.    Table. Further simulated data protocols. . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted May 9, 2023. ; https://doi.org/10.1101/2023.05.07.539770 doi: bioRxiv preprint