MRI analysis to map interstitial flow in the brain tumor microenvironment

Glioblastoma (GBM), a highly aggressive form of brain tumor, is a disease marked by extensive invasion into the surrounding brain. Interstitial fluid flow (IFF), or the movement of fluid within the spaces between cells, has been linked to increased invasion of GBM cells. Better characterization of IFF could elucidate underlying mechanisms driving this invasion in vivo. Here, we develop a technique to non-invasively measure interstitial flow velocities in the glioma microenvironment of mice using dynamic contrast-enhanced magnetic resonance imaging (MRI), a common clinical technique. Using our in vitro model as a phantom “tumor” system and in silico models of velocity vector fields, we show we can measure average velocities and accurately reconstruct velocity directions. With our combined MR and analysis method, we show that velocity magnitudes are similar across four human GBM cell line xenograft models and the direction of fluid flow is heterogeneous within and around the tumors, and not always in the outward direction. These values were not linked to the tumor size. Finally, we compare our flow velocity magnitudes and the direction of flow to a classical marker of vessel leakage and bulk fluid drainage, Evans blue. With these data, we validate its use as a marker of high and low IFF rates and IFF in the outward direction from the tumor border in implanted glioma models. These methods show, for the first time, the nature of interstitial fluid flow in models of glioma using a technique that is translatable to clinical and preclinical models currently using contrast-enhanced MRI.


INTRODUCTION
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. 3  We and others have identified interstitial fluid flow (IFF) as an integral component of the tumor microenvironment. [4][5][6][7][8][9][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. 11 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, 13 hepatic, 6 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. 16 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. 17 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, 18 or identify the fluid drainage volume and pooling rates. 19 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, 21 discriminating between tumor and radiation necrosis regions, 22 and predicting survival time of patients. 23 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.

Spin echo-MRI reveals heterogeneous contrast dynamics
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 T1weighted 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)].

Mathematical model of Gd transport allows us to probe flow field
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.
In vitro and in silico models determine the sensitivity of our computational approach 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, 24 and numerical simulations in which we could control the velocity magnitude and direction in a control volume.

In vitro model
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 lm/s [ Fig. 2 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. 7 Similarly, average model-estimated diffusivity decreased as BME concentration increased, with diffusivity ranges comparable to those previously reported for

In silico model
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 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.
Traditional markers of fluid drainage correlate with our flow patterns Evans blue and other tracer dyes have been used for decades to determine the movement of fluid from tumors. 27 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 tumorassociated 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. 5 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

Measurement of interstitial flow
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. 29 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 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 630 , 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 lm/s, we cannot expect to resolve characteristics of the flow field below a spatial resolution of approximately 400 lm, 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.

Known contributors to transport in the brain tumor microenvironment
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][32][33][34] The use of MRI to measure K trans can be predictive of their tumor response to anti-angiogenic therapies, such as bevacizumab; 35 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.
The effects of the brain tumor microenvironment on interstitial fluid flow 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, 36 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, 37 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. 38 These changes surrounding brain tumors include degradation of the native matrix 39 and increases in non-native extracellular matrix (ECM) such as vitronectin, 40 glycoproteins, 41 Tenascin C, 42 and Collagens. 43 These irregular matrices can alter the permeability and hydraulic conductivity of the overall tissue leading to reduced or increased interstitial flow in regions. 44 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.

Implications in disease treatment and clinical outcomes
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. 47 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.

Cell culture and tumor inoculation
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, 48 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.

MRI acquisition
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.

Computational model and equations
Gadolinium movement within the tumor and interstitium for multidimensional position x ¼ (x, y, z) and time t, the model can be described by the following differential equation: where, at spatial location x and time t, we have /ðx; tÞ is the concentration of Gd mmol , and uðx; tÞ is the velocity field L T Â Ã . Assuming that the intensity of the measured MRI signal, S ¼ Sðx; tÞ, is directly proportional to Gd concentration, /ðx; tÞ, the concentration in (1) can be replaced by MRI signal intensity. Due to the difference in spatial sampling within (104.17 lm) and between (1.68 mm, slice thickness of 0.7 mm and 20% spacing between slices, 0.84 mm from the center of one to the center of the other) MR image slices, each axial slice (x-y plane) was considered independently. Therefore, the continuous-domain model of (1) was discretized in two spatial dimensions, using the forward-time, central-space (FTCS) finite difference method 49 resulting in the following discrete-domain iterative model: where S n ð Þ ðx; yÞ is the measured MRI signal intensity in the axial (x-y) plane at discrete time n. Dt is the time interval between consecutive acquisitions of the same region of interest. Ds is the physical spacing between neighboring MRI pixels (in-plane resolution, identical in x and y directions). r 2 is the discrete approximation of the Laplace operator obtained using the fivepoint kernel r 2 ! While (2) describes the expected temporal evolution of the intensity of each individual pixel given the parameters D n ð Þ ðx; yÞ and u n ð Þ ðx; yÞ, we are interested in the inverse problem: We seek to identify the parameters D n ð Þ ðx; yÞ and u n ð Þ ðx; yÞ, for each location ðx; yÞ and each discrete time n from the temporal changes of the MRI signal under the assumption of the convective transport model in (1). As this system is underdetermined, we employ an optimization approach to identify the solution that minimizes the absolute difference between observed and model-predicted changes in MR signal intensity between two consecutive time points. Under the additional assumption that the diffusion coefficient and the velocity field are constant over the duration of the experiment [D n ð Þ ðx; yÞ ¼ D(x,y) and u n ð Þ ðx; yÞ ¼ u(x,y)], (2) represents a system of N-1 equations (where N is the number of observation time points) and the optimization problem becomes In (3), y ¼ ½S ð2Þ À S ð1Þ ; S ð3Þ À S ð2Þ ; …; S ðNÞ À S ðNÀ1Þ T is a vector of differences in MR signal intensity values between two consecutive temporal frames of a transversal slice, the matrix contains the image-derived values for each temporal frame, and b ¼ ½D; u x ; u y T represents a vector of the unknown parameters. 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.

Preparation of tissue culture insert phantoms
Collagen hydrogels composed of 1.8 mg/ml rat tail collagen I (Corning) and reduced growth factor basement membrane extract (BME, Cultrex) were seeded into 12 mm tissue culture inserts (Millipore, Burlington, MA). A flow condition is created by adding media atop the gel creating a pressure head. [ Fig. 2(b)]. 24 Volumetric flow rates and average superficial velocity were experimentally determined over a period of 30 min as defined by the following equation: where Q is the volumetric flow rate as determined by manually collecting and measuring the total volume flowed through the gel within a given time. A is the cross-sectional area of the tissue culture insert.

In silico phantoms
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 timeevolution 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 À4 cm 2 /s) 26 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.
Tissue post-processing 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 lm. 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 2 , 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 645 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.

SUPPLEMENTARY MATERIAL
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.