Interstitial Fluid Flow and Transport in Glioblastoma and Surrounding Parenchyma in Patients

Background: Glioblastoma is the deadliest, yet most common, brain tumor in adults, with poor survival and response to aggressive therapy. Therapeutic failure results from a number of causes inherent to these tumors. Imaging, computational, and drug delivery approaches can aid in the quest to access and kill each tumor cell in patients. One factor, interstitial fluid flow, is a driving force therapeutic delivery. However, convective and diffusive transport mechanisms are under-studied. In this study, we examine the application of a novel image analysis method to measure fluid flow and diffusion in glioblastoma patients with MRI and compare to patient outcomes. Methods: Building on a prior imaging methodology tested and validated in vitro, in silico and in preclinical models of disease, here we apply our analysis method to archival patient data from the Ivy GAP dataset. Results: We characterize interstitial fluid flow and diffusion patterns in patients. We find strong correlations between flow rates measured within tumors and in the surrounding parenchymal space, where we hypothesized that velocities would be higher. Looking at overall magnitudes, there is significant correlation with both age and survival in this patient cohort. Additionally, we find that tumor size nor resection significantly alter the velocity magnitude. Last, we map the flow pathways in patient tumors and find variability in degree of directionality that we hypothesize in future studies may lead to information concerning treatment, invasive spread, and progression. Conclusions: Analysis of standard DCE-MRI in patients with glioblastoma offers more information regarding transport within and around tumor, can be measured post-resection and magnitudes correlate with patient prognosis.


Introduction
Glioblastoma (GBM), the most lethal form of brain tumor, with a median lifespan postdiagnosis of 12-15 months and a 100% recurrence rate often within several centimeters of the resection cavity. The most recent advancement in GBM therapy was the implementation of concurrent radiotherapy and temozolomide reported by Stupp et al. in 2005[1]. This study pioneered the current standard of care for primary GBM as complete resection, if possible, followed by concurrent temozolomide and radiotherapy administration. Unfortunately, the Stupp protocol only led to a slight improvement in outcomes for patients, increasing the overall 5-year survival to only 27.2% for primary GBM. Despite extensive research efforts since, there has been no advancement in the overall survival of GBM. Although preclinical models exhibit efficacy of several treatment options and therapeutic agents, this success does not translate to clinical trials [2]. A major impediment to the translation of treatments from bench to bedside is the inability to deliver therapeutics effectively within and around the tumor [3].
There are two main modes of drug delivery used in the clinic -systemic and local. Although systemic delivery is often less invasive and easier to implement in the clinic, it must pass through several hurdles before possibly resulting in therapeutically effective response. These challenges include a risk of drug degradation and clearance by the rest of the body, intolerable systemic toxicity, and inability to cross the blood-brain barrier (BBB) at the tumor site [4]. The fact that GBM is a vascularized tumor with abnormal, leaky neovasculature with an often disrupted BBB that can be advantageous for systemic drug delivery [3,5,6]. In fact, this leaky vasculature is exploited to visualize GBM using a combination of paramagnetic contrast agents such as Gadolinium and magnetic resonance imaging (MRI). However, despite the presence of areas with a disrupted BBB, the areas of the adjacent brain vasculature with intact BBB are enough to limit the delivery of drugs to invading tumor cells [7]. Additionally, the high pressures that contribute to a blood-tumor barrier (BTB) and reduce transport of drug from the vasculature into the interstitial space, reduce therapeutic delivery [8][9][10]. Thus, the shortcomings of systemic drug delivery prompted the development of local delivery techniques to bypass the BBB and BTB entirely. Types of local drug delivery include perioperative implant delivery into resection cavity, intraventricular or intrathecal delivery, and convection enhanced delivery (CED) which all aim for therapeutic delivery directly into the tumor and surrounding brain parenchyma [11,12].
The understanding of both the natural and peritumoral flow patterns within the brain are critical for the successful use of both intraventricular or intrathecal delivery and CED, yet the connection between IFF and drug distribution within the brain is understudied [13]. CED emerged as a technique to overcome poor penetration of the tumor and increase therapeutic distribution at the lesion site by creating pressure differentials to increase convective fluid flow directly to the tumor and its surrounding area [14]. Despite being a promising proposal, CED has thus far failed to show significant improvements for patients in clinical trials [14]. One limitation of CED is that infusate can escape the targeted area by following natural flow trajectories within the brain [15]. This observation has reinforced the notion that drug delivery is linked to natural and pathological flow patterns of the brain, which are not always predictable. It is generally thought the increased tumoral pressure drives IFF out of the tumor, into the interstitial space [16,17]. Although this flow pattern is definitely observed, there are also areas with inward flow, parallel flow, and no flow seen at the boundary of a single tumor in implanted murine models [18]. Thus, there is accumulating evidence that although flow patterns are undoubtedly affected by the heightened pressure of the tumor bulk, this does not lead to a single, uniform flow pattern neither intra-nor intertumorally. Hence, having a solid understanding of the mass transport mechanisms, including convection and diffusion, at and around the lesion site is critical to develop effective solutions for the longstanding obstacles in drug delivery.
Fittingly, methods to measure and model these parameters have been a source of growth in the past decades with advances in both imaging via MRI and computational approaches. Diffusion-based imaging techniques, such as Diffusion Tensor Imaging, offer insight into the transport of small molecules throughout the central nervous system [19]. Such advancements in imaging techniques have led to the use of MRI as a tool for estimating drug distribution within the brain. This has been done by creating MRI-visible drug delivery systems, or by correlating specific imaging parameters with drug concentrations at known locations [20,21]. In this study, our goal is to focus on IFF imaging to give further data to computational modelers and drug delivery experts to gain a better understanding of what transport looks like within and around human GBM. DCE-MRI, which utilizes a paramagnetic contrast agent such as Gadolinium, is a well-suited imaging modality to analyze IFF because it allows for the quantitative and noninvasive determination of parameters such as tissue diffusivity and transport within brain tissue [22]. Thus, here we will use DCE-MRI to study the transport of flow and therapeutics within and around tumors.

The Cancer Imaging Archive Ivy GAP Database
Patients with GBM were selected from the Ivy Glioblastoma Atlas Project (Ivy GAP) database, and only analyzed if an axial, T1 weighted DCE-MRI, devoid of motion artifacts was available [23]. We were unable to analyze GBM patients without DCE-MRI availability because this type of MRI is required for our method of IFF analysis. Thus, in this study, we analyzed a subset of 14 patients from the Ivy GAP database who are identified as W13, W18, W29, W30, W31, W33, W34, W35, W36, W38, W40, W43, W48, and W50. Eight of these patients W13, W33, W34, W35, W36, W38, W43, and W48, had DCE-MRI pre-and post-resection available. All data are publicly accessible via The Cancer Imaging Archive [24].

Convection and Diffusion
Analysis: The analysis of IFF in the DCE-MRI acquired from the Ivy GAP database was performed using a computational model previously developed by our group [18]. Assuming that the MR signal intensity is proportional to the contrast concentration within the tissue allows the model to evaluate the spatiotemporal evolution of the contrast agent. This model requires an input of an image stack consisting of at least one pre-contrast, necessary for background subtraction from post-contrast images, and at least three post-contrast images of a single slice, which includes the tumor, from the full brain scan. The graphic user interface (GUI) associated with the model is used to draw a polygon around the region of interest (ROI) (i.e. the tumor) on the image, specify the resolution of the image, timing between the slices of the stack, etc. The model uses the image stack and information input in the GUI to calculate the isotropic diffusion coefficient and velocity field of a ROI by solving the diffusion-advection partial differential equation (PDE) below using the forward-time, central-space finite difference method.
In the above equation, the contrast concentration given by φ(x,t), the isotropic diffusion coefficient D(x,t), and the velocity field u(x,t) evolve in space (x=(x,y)) and time (t). The details regarding the solutions of the above PDE and the model can be found in our previous publication [10]. Using estimates of the spatio-temporal evolution of the contrast agent as input, the model allows us to infer the spatially-resolved diffusion coefficient and the vector field of IFF velocity. The mean and median values of the flow parameters of several slices per tumor were averaged to calculate overall parameter values for the entire tumor. These averages and vector fields are used for the various methods of data visualization presented here.

Statistics and graphing and generation of figures:
Statistics were conducted on individual datasets as described in the results. Graphs were generated using GraphPad Prism v9.0, and graphics were generated using Biorender with a license to JMM. The rose plots were generated using a modified version of the Wind Rose code downloaded from the Mathworks File Exchange, created by Daniel Pereira. The heat maps and images with streamlines were generated using a Python script generated by our group.

Interstitial flow and Diffusion Coefficients can be calculated from DCE-MRI
As we demonstrated in mice earlier, we were able to use gadolinium transport to model both interstitial fluid velocity and diffusion coefficient simultaneously from four sequential images. In mice, a specific sequence was required that took four images after gadolinium entry into the interstitial space over 12 minutes (one image every three minutes). In The Cancer Imaging Archive [24], patient data was available for DCE-MRI image acquisitions which took approximately 1200 images over the course of 2-3 minutes. We chose to analyze a set of images spanning the imaging session from these data and were able to successfully execute our analysis similarly to in mice to determine both IF and DCs in and around the tumors. Our overall process is shown in Figure 1 and includes the acquisition of images from the database, followed by use of the Lymph4D analysis tool [18,25] to generate the data in a pixel-wise fashion, and then subsequent data visualization using Matlab and Python. trates the steps for analysis using patient W13 from the IvyGAP database as an example. The tumor is located on the slice of interest, all timepoints of the slice of interest are extracted from all DCE acquisitions, the tumor is delineated, and analyzed using the Lymph4D analysis tool. The component-wise velocity vectors from the analysis are then input into other Matlab R2020a and Python 3.6 scripts to develop images overlayed with streamline and a quiver plot, as well as the rose plots of velocity magnitude and direction.

Interstitial fluid flow magnitude is variable across patients
Interstitial fluid flow is found throughout tissues, including within and around the tumor [26]. Theoretically, the velocity moving from the tumor bulk and into the surrounding tissue should be highest, as the increased interstitial fluid pressure within the tumor will drive fluid into the surrounding tissue. This phenomenon has been shown in implanted preclinical tumors and in some patient tumors outside of the brain [16,27,28]. Here we analyzed both flow within the tumor, and within the surrounding parenchyma with the hypothesis that velocity would be faster in the parenchyma than within the tumor ( Figure  2A). Six MRI slices per patient were analyzed which encompassed the majority of the tumor in each patient. The average velocity magnitude was calculated per slice and then these were averaged to comprise a total mean velocity magnitude on a per patient basis. Generally, there was about a 10-20% range of mean tumor velocities among the six slices, that was not inherently dependent on location within the brain ( Figure S1).
We saw that there was not a significant difference between the velocity magnitude as measured within the tumor as compared to the surrounding parenchymal space ( Figure  2B). We did see that the rate of flow within the tumor significantly and strongly correlated with that of the parenchyma ( Figure 2C) indicating that the patient-by-patient variability may be more important than the macroregional differences within these individual tumors. To compare, we also wanted to see if diffusion coefficient followed significant trends. We saw that there was no significant difference in diffusion coefficient as calculated in the two regions across patients( Figure 2D), but again we did see a significant, though moderate, correlation between the calculated diffusion within the tumor as compared to the surrounding space ( Figure 2E). We did not observe that the size of the tumor correlated with the velocity (Figure S2A). This is similar to a lack of correlation previously observed in mice.

Patient survival correlates positively with mean velocity magnitude
We aimed to examine the effect of patient-specific variables on flow velocity magnitude within the dataset. We found that the correlation between velocity magnitude and patient weight was nonexistent (r=-0.0055, p=0.984) ( Figure S2B). We did not find that there was a significant difference between sexes (p=0.147), MGMT methylation status (p=0.9497), nor EGFR Amplification (p=0.329) though this dataset may be slightly underpowered to conclude that there is no effect ( Figure S2C). Interestingly, we did find that age significantly correlated with a lower IFF velocity magnitude throughout the tumor ( Figure 3A). This may be explained by a host of literature indicating that fluid flow within the brain slows with age as documented by MRI of ventricles, blood vasculature, and drainage pathways. Most importantly to clinical outcomes, we found that the mean velocity magnitude within the tumor significantly correlated with survival, with higher rates of IFF velocity correlating with longer survival times ( Figure 3B). As expected, age correlated negatively and significantly with survival as well ( Figure 3C) [29,30]. However, to firmly conclude the correlation between velocity and survival, without age as a potential confound, a larger dataset with a range restriction on age would be valuable to examine this novel interaction. Contrastingly, we did not find any correlation between diffusion coefficient and survival in this patient cohort (r=0.182, n.s. Figure S2D).

Resection of tumor does not eliminate interstitial fluid flow
Every patient with glioblastoma undergoes total or subtotal resection after diagnosis. The removal of the tumor decreases the interstitial pressure within the cranium and within the tumor location. Thus, it might be expected that upon removal of the tumor bulk, and thus the source of heightened interstitial pressure, we would see a reduction in interstitial velocity around the tumor into the parenchyma [31]. It appeared that the inherent velocity was a patient-specific parameter more than an interpatient parameter based on our tumor and parenchymal analysis. Thus, we aimed to examine the effect of resection on interstitial velocity magnitude. There was a subset of seven patients in the TCIA Ivy GAP database for which pre-and post-resection DCE-MRI were available ( Figure 4A). Analyzing these patients revealed that there was not a significant decrease nor increase in interstitial velocity magnitude pre and post resection across our cohort ( Figure 4B). However, we did see that there was a change in velocity for individual patients that could be physiologically relevant for better understanding treatment post-resection, with some patients showing decreased flow vs increased flow. Potentially due to this variability in patient response post-resection, we Ant. Post.
did not see a significant correlation between pre-vs. post-resection interstitial velocity magnitude ( Figure 4C). Thus, though there is still inherent flow in the parenchymal space post-resection, the effects of changes in this velocity are unknown within this cohort.

Directional flow velocity is patient-specific
Analyzing IFF in numerous patients revealed inherent variability in IFF directionality within GBM. Some patients have relatively uniform IFF, with little preference for a specific direction, whereas others have a strong tendency to flow in a particular direction. The contrast between these flow patterns are made with analyses for patients W43 (Figure 5 A-D) and W48 ( Figure 5 E-H). The post-gadolinium T1-weighted images of patient W43 ( Figure 5A) and patient W48 ( Figure 5E) indicate the location of these patients' tumors. The streamlines for patient W43 are heavily oriented towards the anterior brain ( Figure  5B), whereas no such clear distinction can be made for the streamlines of patient W48 ( Figure 5F). Furthermore, examining the quiver plot in addition to the streamlines seems to indicate that the areas of faster flow and more directional flow correlate for both patients. This can also be visualized by examining the velocity magnitude heat maps for patient W43 ( Figure 5C) and W48 ( Figure 5G) in conjunction with their corresponding streamlines. Conversely, since faster flows are also observed in the nondominant flow directions, the rose plot of velocity magnitude and direction challenges the notion that IFF magnitude and direction are always correlated ( Figure 5D). Finally, the major advantage of the rose plots is that they effectively represent the directionality vs. uniformity of IFF in patient W43 ( Figure 5D) vs. W48 ( Figure 5H). Although the underlying reasons and accompanying effects of tumoral flow patterns may still be up for debate, here we show both intra-and intertumoral heterogeneity in IFF patterns.

Discussion
Here we have described the use of a previously developed technique to examine interstitial fluid flow within and around the patient glioblastoma microenvironment. Interestingly, we see that interstitial fluid velocities are 10-100 times higher than those we observed in implanted patient tumors in mice [18]. Interestingly, the surface area comparison of total body for human to mouse is ~200, which might be in line with this scaling [32]. Similarly, the surface area of the brain human:mouse scales similarly [33,34]. While the application of such a technique is novel and could potentially give us more insight into a patient's unique tumor in terms of transport parameters, there are limitations. First, our analyses are focused on 2D planes, and we do not have the 3D resolution of the tumor nor transport parameters at this point. Thus, vectors in x and y may be obfuscated or confounded by z direction transport and thus, any holistic approach to applying these data for transport modeling is limited. This may be particularly apparent when examining the inter-slice variability whereby in some patients, there is a wider range of flow magnitudes apparent depending on the slice through the tumor.
Overall, we found that interstitial fluid velocity magnitude correlated with survival in our cohort of archival patient data. As convective flow increased we saw increased survival. We did not see a similar trend with changes in diffusion coefficient across the patients. This suggests that convective transport may be more important in GBM in some regard. Increased convection could lead to better transport of nutrients, drugs, or cells within the tumor. If increased convection is a primary aim of therapeutic delivery strategies, such as CED, this may indicate that better accessed tumors overall would be better treated. We did not examine parameters related to transport into the tumor (such as Ktrans, or vessel permeability) so we do not know if these values would potentially correlate with our measurements of velocity [35,36]. Such parameters would indicate if velocity is simply indicative of increased delivery into the tumor, and not necessarily through the tumor. Reduced values of transport into the tumor are expected with higher interstitial pressures that limit transport across the vasculature and efforts to increase pressure driving from vessels into tumors have shown preclinical success [37,38]. We expect that increased interstitial fluid flow also would reflect increased interstitial fluid pressure since that has been the primary driving force in tumors outside of the brain in preclinical models [15,39,40]. Thus, it is not readily apparent from our data that transport into the tumor and transport through the tumor should correlate, yet this new information should be helpful in further drug delivery efforts.
In fact, the predominating opinion is that increased IFF is a symptom of increased interstitial pressure within the tumor that drives fluid out into the surrounding "normal pressure" parenchyma or tumor microenvironment [28,41,42]. Resection of tumors reduces this interstitial fluid pressure by alleviating the source of increased pressure from the brain. Thus, we expected resection to result in reduced IFF velocity, which it did not. We generally saw similar velocities across patients pre and post surgery, with some cases of reduced or heightened velocities on an individual patient basis. Thus, whether the primary tumor is intact or not, interstitial fluid flow is still occurring within the surrounding tissue. However, one caveat is that we do not know what the normal IFF rate may be in a healthy non-tumor bearing brain. Thus, we do not know if patients who had "lower" velocities never showed increased velocity beyond normal levels, vs patients with heightened velocities. Indeed, measurements of cerebrospinal fluid flow, perivascular and vascular flow both indicate a range of patient-to-patient variability in flows, some of which have been linked to age, but also may be related to other disease and physiological characteristics [43][44][45].
The relation of age with interstitial fluid transport is important since we see reduced velocities in older patients. Drug delivery is intrinsically linked to the diffusion and convection within a tissue. By focusing on IFF, we are specifically looking at the advection component of transport within the brain parenchyma. These convective forces are going to dominate when therapeutics are larger (i.e. antibodies, nanoparticles) as compared to smaller (i.e. small molecules, peptides, etc.). For patients with reduced velocities, which seem to be primarily those in advancing age, we purport that therapeutic delivery may be limited within the tumor, not simply as a matter of delivery to the tumor as one potential reason. However, that is a large leap from the data we show here, but indicates the further need for examination of therapeutics in models of disease that span age as well as sex, molecular subtype, and microenvironment to better address this disease preclinically. Similarly, we know that patient-specific knowledge of tumors and their attributes are of utmost importance due to the heterogeneity of cancer and mathematical oncology efforts to integrate this personalized data are showing success [46][47][48]. Coupling of our imaging methods to better map out therapeutic transport with more advanced mathematical modeling approaches, or coupling with growth models for predictive therapeutics responses could offer another degree of insight into personalized medicine by including IFF as a factor.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1, Figure S1: Mean tumor velocity magnitude by MR slice through tumor; Figure S2: Patient parameters and outcomes on IFF and diffusion