MRI analysis to map interstitial flow in the brain tumor microenvironment

The tumor microenvironment (TME) consists of all cells, extracellular matrix, chemical factors, and biophysical forces aside from the tumor cells. Together, these components create the complete cancer tissue that is both affected by the cancer and can in turn affect the tumor cells.^1,2 The TME has been implicated in therapeutic response, invasion, proliferation, and differentiation of tumor cells. In glioblastoma (GBM), a highly aggressive brain cancer, the microenvironment is known to contribute to the invasion of tumor cells into the surrounding healthy brain.³ This invasion is partially responsible for poor survival seen in patients, as invaded tumor cells cannot be reached by the current standard of care therapy targeting the tumor bulk. Thus, the identification and characterization of mediators of tumor cell invasion could aid in the treatment of GBM.

We and others have identified interstitial fluid flow (IFF) as an integral component of the tumor microenvironment.^4–10, In vitro analyses using microfluidic devices and tissue culture chambers have shown that IFF is involved in increasing proliferation, triggering invasion of tumor cells, and altering the surrounding microenvironment to promote cancer progression. Growing tumors are marked by increased interstitial pressure, due to accumulation of proliferating tumor cells, extracellular matrix, and fluid, which is higher than the pressure in the surrounding tissue.¹¹ This pressure differentially yields increased IFF across the invasive edges of tumors where tumor meets healthy tissue. While it is a potent driver of invasion in brain,^4,5,12 skin,¹³ hepatic,⁶ and breast cancer,^7,9,14,15 IFF has been poorly measured and characterized in vivo. It is thereby pertinent to develop accurate means to quantify IFF in the pre-clinical and clinical settings.

Chary and Jain pioneered the use of intravital fluorescence recovery after photobleaching (FRAP) to approximate fluid velocity in the rabbit ear.¹⁶ However, as FRAP requires optical access, measurements are confined to superficial locations and cannot be acquired quickly across the entire tumor and surrounding microenvironment. Butler et al. implanted micropore diffusion chambers downstream of breast tumors to measure total fluid drainage.¹⁷ While micropore chambers provide good measurements of bulk fluid movement, this method does not afford information on interstitial flow velocities and is difficult to implement in most models. Noninvasive attempts to characterize bulk fluid transport in vivo employ magnetic resonance imaging (MRI). These approaches in implanted brain (intradermal/subcutaneous) and breast tumors (orthotopic) have used multi-compartment models to approximate IFF velocities based on the rate of change of the contrast-enhanced ring at the tumor border over time,¹⁸ or identify the fluid drainage volume and pooling rates.¹⁹ Similarly, other dynamic MRI approaches have estimated fluid velocities in implanted tumor models using equations relating signal intensity to a linear attenuation coefficient.^13,20 Our goal was to improve and expand these techniques by developing a novel methodology to noninvasively measure IFF directly in vivo in GBM. Dynamic contrast-enhanced MRI (DCE-MRI) has been used clinically as a standard imaging method to assess the vascularization of tumors by analyzing the influx of T1 contrast agents (i.e., Gadolinium chelates) and tumor permeability. In GBM alone, DCE-MRI has been used for grading tumors,²¹ discriminating between tumor and radiation necrosis regions,²² and predicting survival time of patients.²³ Here, we take advantage of this common contrast-enhanced MRI technique, to develop our computational methodology to measure IFF in and around brain tumors. Unlike other approaches, we aim to evaluate flow velocities on the basis of biological transport principles. We test our method on several in vitro and in silico phantoms and apply it to four human stem cell xenograft models of GBM. We quantify and map interstitial flow in these models, relating the patterns of flow to a common marker of fluid drainage from tumors, Evans blue dye, to verify that this method is a valid approach for determining regions of interstitial fluid flow in the brain tumor microenvironment.

We utilized MRI to observe contrast dynamics within our tumor models [Fig. 1(a)]. After T2-weighted confirmation and localization of xenograft tumor [Fig. 1(b)] and pre-contrast T1-weighted image for background subtraction [Fig. 1(c)], a bolus injection of gadobenate dimeglumine (Gd) (Multihance, Bracco Diagnostics) was given. Initial influx of intravenous Gd into the tumor was detected by dynamic contrast-enhanced MRI (supplementary material, Video 1), followed by spin echo MRI (T1-weighted) [Fig. 1(d), supplementary material, Video 2]. Four T1-weighted images were successively taken following injection of Gd to identify tissue transport of extravasated contrast agent [Fig. 1(d)]. This technique was found to provide high resolution imaging over our desired time course. Difference maps of final-initial contrast agent localization were used to confirm transport of Gd over the duration of imaging [Fig. 1(e)].

Based on our dynamic magnetic resonance (MR) sequences, we wished to determine the velocities of interstitial flow in vivo. Biological transport principles govern interstitial flow, and we proposed to approximate fluid velocity by assuming that it was equal to the velocity computed for solute movement. We described the movement of our solute, Gd, within the tumor and adjacent brain tissue as a convective process that combines microscopic diffusion and macroscopic bulk motion. As described in more detail in the Methods section, we built a mathematical model of the spatiotemporal evolution of solute, Gd [Eq. (1)], which allowed us to computationally solve the inverse problem, i.e., to estimate the flow velocity field that best explains the observed change in MR signal intensity caused by transport of the Gd contrast agent.

To evaluate our parameter estimation approach, we took advantage of our previously described in vitro hydrogel-tissue culture insert system, in which we could directly measure average velocity,²⁴ and numerical simulations in which we could control the velocity magnitude and direction in a control volume.

We have optimized our tissue culture insert system to study pressure-driven flow through a tissue-mimicking hydrogel.^24,25 This system consists of a tissue culture insert with a collagen I gel and varying concentrations of basement membrane extract (BME) to tune the overall volumetric flow [Fig. 2(a)]. By applying a pressure head atop the hydrogel, fluid is driven through the gel and into the lower chamber over the course of 18 h. We imaged this system as the contrast agent moved through the gel over time during the washout of the agent with non-contrast containing media, after addition of Gd to the upper chamber. MR and analysis revealed that Gd flow exhibited velocities between 0 and 4.5 μm/s [Fig. 2(b), top panel] and was primarily downward directed [Fig. 2(b), middle and bottom panel] which was the expected direction. Unexpected variability in estimated flow directions in the upper central part of the gel is likely due to early washout of Gd in that region [Fig. 2(b), middle panel]. As a result of washout, little to no change in MR signal intensity is detectable during subsequent image acquisitions, limiting the ability of our reconstruction algorithm to reliably identify optimal parameters. By selecting a time window within the MR images where the concentration of Gd was changing throughout the course of acquisition, we could determine an average velocity magnitude through the hydrogel. Using gels with increased concentration of BME reduced the model-calculated average velocity magnitude [Fig. 2(c)], in line with previously published studies.⁷ Similarly, average model-estimated diffusivity decreased as BME concentration increased, with diffusivity ranges comparable to those previously reported for collagen gel systems [supplementary material Figs. 1(a) and 1(b)].²⁶ Average model-estimated velocity magnitudes ranged between 1.37 and 2.24 μm/s [Fig. 2(c)] for the different gels and were significantly correlated to manually measured average velocities (Q/A) of 1.17–3.52 μm/s in the same systems [Figs. 2(d) and 2(e)].

To further test our estimation approach, we computationally solved the forward model [Eq. (1)] for specific flow scenarios, which allowed us to compare the estimated to actual flow parameters. We generated computational phantoms mimicking idealized flow situations of bidirectional [Figs. 3(a) and 3(b)] and multidirectional (360°) flow [Figs. 3(f) and 3(g)] of spatially constant magnitude, and unidirectional flow of spatially varying magnitude [Figs. 3(k) and 3(l)]. Evaluation of our reconstruction approach on these in silico phantoms showed that the algorithm performs well on flow fields of spatially constant velocity magnitude [Figs. 3(b) and 3(g)]. The direction of flow was reconstructed correctly, locally throughout the domain [Figs. 3(c) and 3(h)] as well as globally [Figs. 3(e) and 3(j)]. Flow velocity magnitude was consistently reconstructed to approximately 50% of the simulated velocity, except for boundary regions where a larger relative error was observed. The algorithm performed less well on fields with spatially variable velocity magnitude [Fig. 3(l)], where it overestimated the flow velocity magnitude in regions of slower flow [Figs. 3(m) and 3(n)] up to a factor of greater than or equal to 10 at the domain boundary. However, despite differences in the velocity magnitude and local velocity direction [Fig. 3(m)], the overall flow direction in this scenario was reconstructed within ±30° [Fig. 3(o)].

Using our MRI-based flow reconstruction approach, we aimed to characterize and quantify interstitial velocities in and around murine xenografts of glioblastoma stem cell (GSC) lines. Either the tumor alone, or combined tumor and surrounding interstitium, where we saw changes in the contrast agent, was first demarcated on the first post-contrast T1-weighted image to identify the region for analysis [Fig. 4(a)]. We then applied our reconstruction approach in this region to estimate IFF velocity maps from which velocity magnitude [Fig. 4(b)] and direction [Fig. 4(c)] were derived. Both the velocity magnitude and direction were highly variable across each individual tumor, indicating a heterogeneous frequency of interstitial velocity that was not always in the outward direction [Fig. 4(d)]. The distribution of velocity magnitudes within each tumor was positively skewed and, interestingly, was similar between tumor models [Figs. 5(a) and 5(b), supplementary material Figs. 2(a)–2(d)]. This was also true when the region outside of the tumor bulk was included (the interstitial region) [Figs. 5(a) and 5(c), supplementary material Figs. 2(a)–2(d)]. The distribution of diffusivities was normally distributed, but again similar between models [supplementary material Figs. 3(a), 3(b), and 4]. Comparison of the tumor bulk interstitial velocity magnitudes to the tumor bulk + interstitial space indicated that the average flow rates in these interstitial regions were not significantly different from those in the bulk alone [Fig. 5(d)]. Neither the mean [Fig. 5(d)] nor the median [supplementary material Fig. 2(e)] calculated interstitial velocity magnitudes were significantly different between tumor models. Similarly, neither the mean nor the median calculated diffusivities were significantly different between tumor models, except the average bulk diffusivity between G62 and G528 in which G528 was significantly higher [supplementary material Figs. 3(c) and 3(d)]. Since there was a high degree of variability between the tumor-bearing mice, we decided to examine the correlation of outcomes on a mouse-by-mouse basis to identify possible contributing factors to interstitial flow velocity differences. The tumor size was variable at the same timepoint across models, but, we did not observe a significant correlation between the tumor size and the mean interstitial velocity magnitude across tumor models [Fig. 5(e)], nor was there a correlation between the average model-calculated velocity magnitude and average model-calculated [supplementary material Fig. 2(f)]. Thus, it appeared that intratumoral heterogeneity was more apparent in the four tumor models we examined, than intertumoral variability.

Evans blue and other tracer dyes have been used for decades to determine the movement of fluid from tumors.²⁷ We have used it frequently to identify histological regions of interstitial fluid flow within the tumor microenvironment.^4,5 Upon intravenous injection, Evans blue binds to albumin in the bloodstream and only crosses into the brain through compromised tumor-associated vasculature, thus allowing us to uniquely examine transport of this molecule from the tumor into surrounding parenchyma, similar to Gd. To determine how our reconstruction approach matches with this histological method, we wished to correlate Evans blue intensity with our calculated outcomes. We hypothesized that regions with higher Evans blue intensity would have both a net outward direction to the velocity vector fields and higher velocity magnitudes than in regions of lower Evans blue intensity. Multiple user-defined regions of Evans blue positivity were selected [Figs. 6(a)–6(c)] and defined as high positivity if the intensity of the region was greater than 2/3 of the maximum integrated density measured for the tumor, similar to previous approaches we have used.⁵ We also normalized intensity to the maximum intensity measured for the tumor. Then, the average interstitial velocity of the region was computed using our method on the matched MRI section [Figs. 6(b) and 6(c)]. Flow patterns were grouped into either outward flow (±45° from the normal to the tumor border at that point) or non-outward (indicating either inward flow, or parallel flow) with regard to the tumor border. Regions with outward flow [Fig. 6(d)-i, iii, and iv] were characterized by higher Evans blue intensity in histological sections when compared with the integrated density in regions of non-outward flow [Fig. 6(d)-ii] for 10 out of 12 tumors analyzed [Fig. 6(e)]. Calculation of the average velocity magnitude in these regions indicated that regions with high Evans integrated density showed average velocity magnitudes that were higher than those in low Evans blue integrated density [Fig. 6(f)]. Thus, our flow patterns validate the histological outcome of Evans blue drainage within and around the tumor.

The biophysical parameter of IFF velocity has not been widely studied in vivo, though in vitro analyses suggest that it is an important contributor to tumor cell invasion, proliferation, and shaping of the cancer microenvironment.^{4,12,14,15,28} These in vitro systems have been designed to establish IFF velocity magnitudes in line with those reported from a range of different platforms including computational modeling, extrapolation from in vivo data of bulk flow accumulation, and real-time or endpoint measurements of tracer molecules. Each of these methods has their strengths and limitations as discussed earlier, but overall, they lack the ability to examine the heterogeneity within the tumor and to directly measure interstitial fluid velocity using biotransport principles.

We utilized a commonly employed in vitro interstitial flow system to evaluate velocity patterns generated using our computational reconstruction approach and saw flow in the expected downward direction and on the same order of magnitude as the measured average velocity from the volumetric flow rate. These values changed as the matrix was synthesized to be more hydraulically resistant. Computed diffusivity measurements of the in vitro phantom as derived from the MRI method (supplementary material Fig. 3), though not the focus of the present study, were similarly on the same order of magnitude as previously reported values.^25,26 Future work can validate this component of the model using more standard MRI methods such as diffusion-weighted MRI as a basis of comparison.²⁹ Although the in vitro phantom serves well for assessing average values, spatial accuracy of the reconstructed velocity field could not be assessed in this setup. Therefore, we used in silico phantoms of three idealized flow scenarios to estimate the spatial accuracy of our reconstruction algorithm (Fig. 3). These experiments indicate overall good and spatially consistent performance for multidirectional flow fields of constant velocity magnitude, although based on our phantom results, we believe that we may be underestimating this magnitude by a factor of 2 [(Figs. 3(d), 3(e), 3(i), and 3(j)]. However, we caution about the interpretation of reconstructed velocity magnitude when the flow speed varies spatially [Figs. 3(m) and 3(n)]. In the latter scenario, we found the reconstructed velocity field to differ from the actual field not only by a spatially constant factor of magnitude but in a space-dependent manner, thus rendering comparisons between subregions difficult. As the fluid flow between the tumor and interstitium likely features a combination of the flow characteristics examined here, this phenomenon clearly limits the accuracy of our reconstruction algorithm when applied to in vivo data. Despite relatively large but very localized reconstruction errors in the velocity magnitude [Figs. 3(m) and 3(n)], our in vitro and in silico experiments indicate that reconstructed velocity fields are sufficiently accurate to estimate the bulk flow direction to approximately ±30°, even in boundary regions. This allows inward and outward flow to/from tumor from/to interstitium to be clearly distinguished in most situations.

The spatial accuracy at which the flow velocity fields can be reconstructed is intrinsically limited by the time interval between subsequent MR images, in our case approximately 3 minutes. Assuming an average IFF velocity magnitude of 2 μm/s, we cannot expect to resolve characteristics of the flow field below a spatial resolution of approximately 400 μm, corresponding to about 4 pixels. Therefore, spatial accuracy should not be affected by our choice to reconstruct the velocity in each pixel based on the intensity values from its 3 × 3 pixel neighborhood.

Since our reconstruction approach aims to estimate IFF flow indirectly from the change in concentration of an advected compound, its application is limited to situations in which such a change is detectable. In our experiments, parts of the tumor will either be saturated or void of contrast agent during a fraction of dynamic contrast-enhanced image acquisition. The lack of detectable concentration change in those regions results in less accurate reconstruction of the flow velocity field. This effect can be seen in the velocity field reconstructed from the MR phantom [Fig. 2(b)] where the velocity magnitude and direction are inconsistently estimated at the upper edge of the meniscus that remains void of contrast agent during most image acquisitions. Thus, the selection of appropriate concentrations and non-saturated locations is important to the overall accuracy of the reconstruction.

Interstitial flow has been tied to the increase in interstitial pressure within a confined tissue. It has been documented that larger tumors have higher interstitial fluid pressures,^11,30 but we did not see a significant correlation between the tumor size and the average interstitial velocity magnitude. Although we would expect the fluid velocity to increase with the size of tumor, it is possible that in the confined space of the brain, there is an upper limit to the velocity magnitude, and thus, further work examining changes in intratumoral velocity during progression would yield these types of dynamic measurements. There are other factors that we did not examine that may contribute to changes in the average tumor velocity. The correlation between blood vessels and total fluid leakage and perfusion into tumors has been well-documented.^31–34 The use of MRI to measure K_trans can be predictive of their tumor response to anti-angiogenic therapies, such as bevacizumab;³⁵ however, these studies never examined the movement of this fluid once in the tumor. Leaky vasculature is required for entry of Gd contrast into the space, but our methodology specifically examines the movement of this agent once it has extravasated into the brain tumor, and thus while the initial signal is dependent on the vasculature density and permeability, only a change in signal over time is required for the overall measurement.

As stated, high levels of intratumoral heterogeneity contribute to the inability to have definitive correlative findings between IFF and distinctive xenograft tumor models or size. Classically, interstitial fluid flow has been defined through modeling as moving in the outward direction from the tumor into the healthy tissue; however, we see that it is much more complex than that. Interstitial flow varies spatially throughout the tumor, and while this is unsurprising, the inward direction of flow in some border regions was unexpected. Indeed, these regions of inward flow correlated with lower Evans blue intensities [Fig. 6(e)], as did regions of lower average interstitial velocity magnitudes [Fig. 6(f)]. Evans blue has been used as a tracer of bulk drainage for decades both clinically and experimentally. We have used it in previous works to correlate regions of higher and lower flow based on assumptions that interstitial fluid drainage would correspond to the movement of this tracer and thus indicate the general trajectory of interstitial fluid velocities.^4,5 Here, we demonstrate that these assumptions are valid by both depicting the magnitude of flow and to a lesser extent the direction of the flow corresponds to accumulation of this tracer. However, the relationship between this interstitial fluid drainage and bulk fluid drainage has not been specifically probed with our methodology; here, it stands to reason that this interstitial fluid eventually drains towards classical bulk flow structures within the brain.

Natural cerebrospinal fluid flow pathways in the brain may dictate tumor fluid behavior. In healthy tissue, bulk fluid flow pathways were identified almost 20 years ago along the white matter tracts,³⁶ and natural movement of fluid towards these bulk pathways along white matter tracts may result in non-homogenous flow outward from the tumor. Similarly, in the glymphatic system, fluid drains along perivascular spaces,³⁷ and thus tumor fluid flow, may be directed towards and along the nearby blood vessels leading to heterogeneous flow. Finally, tumor cells alter the surrounding extracellular matrix themselves or through the activation of surrounding stromal cells that then alter the matrix.³⁸ These changes surrounding brain tumors include degradation of the native matrix³⁹ and increases in non-native extracellular matrix (ECM) such as vitronectin,⁴⁰ glycoproteins,⁴¹ Tenascin C,⁴² and Collagens.⁴³ These irregular matrices can alter the permeability and hydraulic conductivity of the overall tissue leading to reduced or increased interstitial flow in regions.⁴⁴ Thus, changes to the TME can alter flow and, in turn, the flow can alter the TME, similar to other mechanobiological forces.^45,46 Understanding this feedback loop as it relates to the interstitial flow is important for dissecting out new mechanisms of cancer progression, not only in the brain but also in other organs as well.

Magnetic resonance imaging is routinely used in GBM diagnosis, surgery, and postoperative monitoring and this straightforward imaging technique has been implemented in the clinical and pre-clinical settings for the study of GBM.⁴⁷ Thereby, we have developed and used MRI methodology coupled with biotransport modeling to map interstitial flow in GBM xenografts. As IFF has been linked to increased invasion in GBM and areas of flow have correlated with invasive fronts in mouse models, we believe that our method could eventually act as a predictive indicator of regional invasion. This could allow physicians to aggressively target certain areas with surgery or anti-tumor chemotherapies, or benefit by identifying transport pathways to utilize during drug delivery. Finally, we can use MRI-derived IFF maps to determine if IFF is correlated with microenvironmental changes or activation of effectors that are leading to progression or worsened disease. Identification of such changes could act to identify new therapeutic targets. Assessment of IFF fields gives us the ability to study the tumor microenvironment in the context of an understudied biophysical force: interstitial flow, and thus, potentially open avenues for new understanding of glioblastoma.

All animal experiments were approved by the University of Virginia Institutional Animal Care and Use committee under protocol number 4021. Briefly, Glioblastoma stem cells (GSCs) were cultured and maintained as previously described,⁴⁸ and then steriotactically injected into 8–10 week old male Nonobese diabetic severe combined immunodeficiency (NOD SCID) mice with 15 000 GSCs [G2 (n = 3), G34 (n = 5), G62 (n = 4)] or 400 000 GSCs [G528 (n = 4)]. Injection coordinates were 2 mm lateral and posterior to bregma at a depth of 2.7 mm.

Ten or eleven days post inoculation, mice were anesthetized, and tail vein catheters were inserted. Mice were then imaged with a 7T Clinscan system (Bruker, Ettlingen, Germany) equipped with a 30 mm diameter cylindrical RF coil. The presence of tumor was confirmed using a T2-weighted sequence with the following parameters: repetition time (TR) = 5500 ms, echo time (TE) = 650 ms, field of view (FOV) = 20 mm × 20 mm with a 192 × 192 matrix, slice thickness = 0.5 mm, number of slices = 30, and two averages per phase-encode step requiring a total acquisition time of about 5 min per sequence. Subsequent pre-contrast T1-weighted imaging was performed. A bolus injection of gadobenate dimeglumine was administered at a concentration of 0.3 mmol/kg in sterile saline (MultiHance, Bracco Diagnostics, Princeton, NJ). Following injection, a series of four post-contrast T1-weighted images were taken axially through the head with the following parameters: repetition time (TR) = 500 ms, echo time (TE) = 11 ms, field of view (FOV) = 20 mm × 20 mm with a 192 × 192 matrix, slice thickness = 0.7 mm, number of slices = 22, and two averages per phase-encode step requiring a total acquisition time of about 3 min per sequence. T1-weighted images were acquired for approximately 13 min post-injection.

The solution of (3) is computed separately for each voxel in a selected MR slice after imposing the additional $D≥0$ constraint on the diffusion coefficient, yielding voxel-wise optimal values for each parameter of the transport model. Derived quantities (diffusion coefficient, velocity magnitude, velocity direction) are computed from these optimized parameter maps.

The optimization problem (3) was originally solved for each voxel independently of the behavior of its neighbors occasionally resulting in very sharp transition within the parameter maps. To provide a tighter coupling between the solutions evaluated at neighboring pixels, the optimization process was modified to extend $y$ and $A$ to include values from within a square region around the original voxel under analysis effectively utilizing all temporal values for all voxels within the region as part of the optimization process.

In silico models were created to investigate the accuracy of the proposed reconstruction algorithm and compare estimated to actual model parameters by computationally solving (FENICS) the forward model (1) for specific flow scenarios. We examined three scenarios: standard scenarios: bidirectional flow [Fig. 3(b)] multidirectional (360°) flow of spatially constant magnitude [Fig. 3(g)], which replicates an idealized tumor model with pressure in the outward direction, and unidirectional flow of spatially varying magnitude [Fig. 3(l)].

For each scenario, an in silico phantom was created by defining an initial distribution of solute (Gd) ϕ(x, y, t = 0) [Figs. 3(a), 3(f), and 3(k)] and computationally simulating its time-evolution according to (1) for prescribed velocity fields, u(x, y), and spatially constant diffusivity, D(x, y) = D. From the solution of this forward model, maps of Gd concentration ϕ(x, y, t = t_n) were extracted at multiple evaluation time points t_n, analogous to MRI scenarios with spatiotemporally evolving contrast enhancement. Flow velocity was estimated from these synthetic images by applying our reconstruction algorithm. Finally, we compared the direction and magnitude of prescribed and reconstructed velocity fields for each simulated flow scenario.

The forward model describing the time evolution of the concentration was solved in 2D using the Finite-Element Method. All simulations assumed a spatially constant isotropic diffusion coefficient (1 × 10⁻⁴ cm²/s)²⁶ and were solved in time steps dt = 0.01 s. Solutions were extrapolated onto a 50 × 50 pixel regular grid to mimic the image size of large tumors at typical small-animal MR in-plane resolution.

Following MRI acquisition (ten or eleven days after inoculation), animals were injected with Evans blue dye administered intravenously. The following day, animals were euthanized by intracardial perfusion with saline and brains were cryoembedded and sectioned at 12 μm. Sections at varying depths within the tumor were immunostained for mouse anti-human nuclei (clone 235-1, Millipore) and 4',6-Diamidino-2-Phenylindole, Dihydrochloride (DAPI, Sigma). Whole section scans were taken using an EVOS FL Auto 2.0 and images processed using Photoshop for figures. Raw images were imported into ImageJ, and four to five user-defined regions of Evans blue positivity were selected and intensity was measured using integrated density. Regions were approximately 0.49 mm², and the integrated density in each region was normalized to the tumor maximum. Regions were defined as ‘high flow’ if the integrated density was greater than 2/3 the maximum for the tumor. Using DAPI and anti-human nuclei staining to identify tumor size and location, sections were matched to closest corresponding MRI slice. The velocity was computed in the matched regions, and average velocity and direction were determined. Briefly, flow was considered ‘outward’ if the average direction was ±45 degrees from the normal to the tumor border, and ‘non-outward’ if the direction was parallel. Statistical analysis was performed in GraphPad to compare Evans blue intensity values with velocity direction and magnitude. Wilcoxon matched-pairs signed-rank tests were performed, to correct for non-parametric distributions, to compare outward vs non-outward flow, and low and high flow for each tumor.

See supplementary material for the summary data and histograms of the in vitro and in vivo calculated diffusivities, as well as the histograms of the in vivo calculated velocity magnitude for the individual tumor models.

The authors would like to thank B. J. Purow, MD, for supply of the glioblastoma stem cells, Jack Roy and Dr. Stuart Berr with the University of Virginia Molecular Imaging Core for their assistance with MRI sequence development and optimization, Daniel Logsdon and Alexander Matthew for helpful discussion regarding model development, and Jeremy Gatesman with the University of Virginia Center for Comparative Medicine for his assistance with tail vein injections. We also thank Desiree Floyd for assistance with tumor injections, Rebecca Della Croce for her help with Evans blue intensity analysis, and Melissa Brevard with the University of Virginia Cardiovascular Center histology core for her help with cryosectioning. This work was supported in part by ACS-IRG-81-001-29, NIH 1R01CA222563-01 to J. M. Munson and funding to K. M. Kingsmore by the NSF-GFRP. Additionally, research reported in this publication was supported by the National Cancer Institute of the National Institutes of Health under award number P30CA033572 to City of Hope and funding to D. Abler from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agreement No. 753878. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.