Modelling the transport of fluid through heterogeneous, whole tumours in silico

Cancers exhibit spatially heterogeneous, unique vascular architectures across individual samples, cell-lines and patients. This inherently disorganised collection of leaky blood vessels contribute significantly to suboptimal treatment efficacy. Preclinical tools are urgently required which incorporate the inherent variability and heterogeneity of tumours to optimise and engineer anti-cancer therapies. In this study, we present a novel computational framework which incorporates whole, realistic tumours extracted ex vivo to efficiently simulate vascular blood flow and interstitial fluid transport in silico for validation against in vivo biomedical imaging. Our model couples Poiseuille and Darcy descriptions of vascular and interstitial flow, respectively, and incorporates spatially heterogeneous blood vessel lumen and interstitial permeabilities to generate accurate predictions of tumour fluid dynamics. Our platform enables highly-controlled experiments to be performed which provide insight into how tumour vascular heterogeneity contributes to tumour fluid transport. We detail the application of our framework to an orthotopic murine glioma (GL261) and a human colorectal carcinoma (LS147T), and perform sensitivity analysis to gain an understanding of the key biological mechanisms which determine tumour fluid transport. Finally we mimic vascular normalization by modifying parameters, such as vascular and interstitial permeabilities, and show that incorporating realistic vasculatures is key to modelling the contrasting fluid dynamic response between tumour samples. Contrary to literature, we show that reducing tumour interstitial fluid pressure is not essential to increase interstitial perfusion and that therapies should seek to develop an interstitial fluid pressure gradient. We also hypothesise that stabilising vessel diameters and permeabilities are not key responses following vascular normalization and that therapy may alter interstitial hydraulic conductivity. Consequently, we suggest that normalizing the interstitial microenvironment may provide a more effective means to increase interstitial perfusion within tumours.

Introduction Architectural heterogeneities in cancerous tissue limit the delivery of anti-cancer drugs by inhibiting their ability to circumnavigate the entire tumour to all cancerous cells [1]. In solid tumours, drug penetration to the tumour core is hindered by physiological barriers which can limit the delivery of targeted agents, with penetration exacerbated by the size of the agent [1][2][3][4][5]. Consequently, preclinical tools which provide a better understanding of therapy interactions within the tumour microenvironment are urgently required in order to increase treatment efficacy. In silico modelling is one such tool which can meet this need by testing novel therapeutic strategies at a much faster rate and cheaper cost than preclinical experimentation [6].
For a systemically-administered agent to effectively target diseased tissue, it must travel from the site of delivery to the site of disease, whilst minimally interacting with normal tissues and not degrading [7]. This is difficult to achieve in tumours since atypical endothelial proliferation of tumour vasculature leads to spatial variations in vascular density and branching patterns, distorted and enlarged vessels, and a highly convoluted network topology [8][9][10]. Further, vascular permeability is heightened and heterogeneous and so these immature blood vessels are generally leakier than those in normal tissue [3,11].
The irregular microenivronment is typically characterised by hypoxia, acidosis and elevated interstitial fluid pressure (IFP) [12][13][14], which drive both tumour vascular proliferation and resistance to therapy [15]. Here, drug delivery may be hindered by the atypical nature of the tumour interstitium. The extracellular matrix (ECM) consists of a cross-linked dense network of collagen and elastin fibres, far denser than usually seen in normal tissue [16]. A denser matrix can inhibit oxygen and nutrient delivery, as well as providing significant resistance to the advection and diffusion of therapeutic particles [1], since key determinants of intratumoural fluid and mass delivery include tissue hydraulic conductivity and vascular compliance [17]. Several therapeutic interventions have sought to limit the effects of these physical barriers by manipulating the microenvironment to enhance the delivery of macromolecular agents [16,18]. For example, normalising the tumour vasculature to reduce vessel permeability thereby increasing drug penetration [12]; and manipulating the connective tissue, and therefore interstitial hydraulic conductivity, using a platelet-derived growth factor (PDGF) antagonist to reduce tumour IFP [19].
Heterogeneities in the underlying morphology of tumours, such as vessel diameters and lengths, and inter-branch distance, exist across individual tumours and tumour cell-lines [20]. These variations in tumour architecture lead to spatial variability in drug efficacy, which complicate efforts to design effective treatment strategies [7]. Experimental efforts have been made to understand the effects of tumour heterogeneity on fluid interactions across tumours, for example, wick-in-needle has been used to measure IFP across tumours [21][22][23]. However, this method disturbs the local microenvironment and only provides an IFP measurement at individual locations. Non-invasive methods have also been developed to estimate tumour IFP [24,25]. For example, convection-MRI, which measures low-velocity flow in tumours at a resolution of � 250 μm in vivo [26]. However, these methods fail to capture full spatial maps of flow at the micron-scale which are crucial to understanding how the combined intra-and extravascular spatial flow heterogeneities occurring at the scale of blood vessels affects the macro-scale flow dynamics and consequent delivery of drugs within a solid tumour. Biomedical imaging complemented by in silico methods provides scope to provide such detail.
Mathematical models have been used to investigate the tumour microenvironment and have provided detailed insights which may otherwise be unavailable experimentally. Seminal models have indicated that a leaky tumour vasculature induces elevated IFP, reduced fluid penetration into the interstitium [14,27], and a non-uniform distribution of drug delivered to solid tumours [2,3,11]. Further, they have defined conventional IFP profiles in tumours-a uniform pressure at the core, with a large decreasing gradient towards the periphery. However, these models generally average spatially over the tumour vasculature and so fail to capture the micron-scale flow dynamics; and they assign a fixed pressure boundary condition on the periphery of the tumour which may artificial induce these conventional IFP profiles.
Subsequent studies have sought to incorporate the spatially heterogeneous effects of tumour vasculature using computer-generated synthetic networks which retain key features of tumour vascular architecture [28][29][30][31][32][33][34], or by integrating spatial variations in vascular permeability parametrised against in vivo experimentation [35,36]. More recently, a hybrid image-based framework has been developed which combines realistic tumour vascular architectures and in silico modelling to predict tumour vascular blood flow [37] and intravascular tumour oxygenation [38]. However, computational models are urgently required which incorporate these highly detailed data to enable predictions of fluid and mass transport to the surrounding tissue [6].
Recent advances in ex vivo optical imaging of cleared tissue specimens have enabled large samples (up to 2 cm 3 with > 10 5 blood vessels) to be imaged in three-dimensions, at resolutions down to a few microns [39]. We have developed a platform called REANIMATE (REAlistic Numerical Image-based Modelling of biologicAL Tissue substratEs) which combines optical imaging of cleared tissue with mathematical modelling and in vivo imaging, within a unified framework, to generate quantitative, testable predictions regarding tumour transport [40]. The platform uses high-resolution imaging data from large, intact, optically-cleared tissue samples to make in silico predictions of blood flow, vascular exchange and interstitial transport. REANIMATE enables new hypotheses to be generated and tested in a manner that would be challenging or impossible in a conventional experimental setting. We have previously used REANIMATE to explore the impact of vascular network topology on fluid transport and vascular disrupting therapy (Oxi4503) to two colorectal cell-lines (LS147T and SW1222) [40].
We develop here a computational model to efficiently simulate both intra-and extravascular fluid transport across large, discrete microvascular networks [40]. Our model simulates Poiseuille flow through the vasculature using the optimisation scheme of Fry et al. [41], parametrised and validated against in vivo ASL-MRI data [40]. Following a similar Green's function method for oxygen transport [42], the vascular component is coupled, via a discrete set of point sources of flux, to a Darcy model which simulates the effective fluid transport in the porous interstitium. A linear system is formed whereby only the source strengths need to be resolved, making it more computationally efficient compared to finite difference or element methods which require a spatial, numerically-discretised mesh [42].
In this study, we detail the generalised model which allows for spatially heterogeneous hydraulic conductances and conductivities. A comprehensive description is provided of its application to whole tumour vascular networks. We then apply our model to an orthotopic murine glioma (GL261) and a human colorectal carcinoma xenograft (LS147T) and reproduce physiological conditions observed in literature. We perform sensitivity analysis to the model parameters associated with transvascular fluid delivery, such as vascular hydraulic conductance and interstitial hydraulic conductivity, to explore the impact on the tumour IFP and interstitial fluid velocity (IFV) profiles.
Finally, we use our computational framework to explore the biomechanics underpinning vascular normalization. Vascular normalization is a method that applies anti-angiogenic therapy to restore tumour vascular structure and function to physiological levels [9,[43][44][45]. By modifying vascular architecture, therapy aims to normalize tumour perfusion and oxygenation, thereby increasing the efficacy of chemo, radio and immunotherapy [44,[46][47][48].
Preclinical and clinical evidence indicates that anti-VEGF (vascular endothelial growth factor) therapy creates a transient window of vessel normalization which improves tumour oxygenation and the delivery of therapeutic agents [15,49]. However, the extent and window of normalization, and the duration and dosage of an anti-cancer drug varies with tumour type [18]. Further, the gold standard for detecting normalization of tumour blood vessels (such as perfusion, microvessel density, morphology and permeability) in the clinic is via histological staining [50].
Here, we recruit our in silico model in combination with realistic, static vasculatures to replicate fluid dynamics changes observed experimentally during vascular normalization [51][52][53] by modifying parameters such as blood vessel diameters, and vascular and interstitial permeabilities. In doing so we hypothesise which of these biomechanics are responsible for experimental observations following therapy, in addition to how tumour type, and the inherent differences in vascular networks structures, affects tumour IFP and perfusion following normalization therapy.

Acquisition and processing of real-world tumour datasets
Orthotopic murine gliomas and human colorectal carcinoma xenograft from the GL261 and LS147T cell-lines (n = 6 for each), respectively, were grown subcutaneously in 8-10 week old, female mice. Following 10 to 14 days of growth, in vivo arterial spin-labelling MRI (ASL-MRI) was performed on a subset of GL261 and LS147T tumours, from which a mean tumour perfusion of 130 ± 50 and 19 ± 8 ml/min/100g was measured [40], respectively. Following perfusefixation, tumours were harvested, optically cleared and imaged using optical projection tomography (OPT, Bioptonics, MRC Technologies, Edinburgh). All experiments were performed in accordance with the UK Home Office Animals Scientific Procedures Act 1986 and UK National Cancer Research Institute (NCRI) guidelines [54]. Full details of the experimental protocol is provided in d'Esposito et al. [40].
Whole-tumour blood vessel networks were segmented from the OPT data for both tumour types. A combination of three-dimensional Gaussian and Frangi filters were applied to the data to enhance vessel-like structures allowing for the segmentation of the blood vessels from the background (see Fig 1a). Skeletonisation of these thresholded data was performed in Fig 1. (a) An example of an SW1222 tumour vascular network enhanced using Frangi filters and extracted from the tumour image stack generated by OPT. (b) The skeletonised vasculature is then segmented into a series of interconnected nodes and vessel segments with known diameters, d i and lengths, l i , for i = 1, 2, 3 shown here. (c) A schematic of sensitivity analysis performed on the source parameters: 1) updating intravascular pressure p i,b for k iterations where p 0,b is the initial network pressure distribution approximated by the flow estimation algorithm; 2) the spacing δ between sources distribution across a branching vessel; and 3) the size of the source radius, r 0 . (d) A flow diagram of the computational framework. In vivo imaging is performed on vascularised tissue to obtain perfusion data (and literature values of vascular pressures when available) which are used to parameterise and validate the framework. Ex vivo imaging is performed on equivalent tissue samples to obtain data on the vascular architecture, including coordinates, vessel diameters and lengths, which are then used to parameterise the vascular flow model. Boundary conditions are assigned (see Fig 2a and 2b) and network intravascular blood pressure is solved. Sources of fluid flux are distributed across the vasculature and assigned a radius equal to their corresponding vessel radius. Interstitial flow parameters are assigned and the model is coupled to the vascular flow compartment via Starling's Law. Solved source strengths are used to update Starling's law (i 2 N s ). This iterative scheme is terminated once predefined tolerances are reached. Amira (Thermo Fisher Scientific, Hillsboro, OR), which also converted the data into graph format (interconnected network of nodes and segments with associated radii, see Fig 1b). To ensure that vessel structures were accurately represented, three-dimensional networks were visually inspected against two-dimensional imaging slices for an accurate representation of vessel location and thickness. Full details of the validation can be found in the Supplementary Material of d'Esposito et al. [40].
In this study a GL261 and a LS147T tumour network were chosen from the d'Esposito et al. [40] datasets for in silico development and testing. Vessel diameters ranged from 17.9 ± 9.3 and 8.9 ± 2.8 μm, with branching lengths of 68.7 ± 48.3 and 88.8 ± 49.4 μm, respectively (see Table 1). Vessel branching angles, inter-vessel distance, radii and tortuosity measures were consistent with data from previous studies that extracted vascular architectures using different methods [20,40].

Computational model
Our computational framework is compartmentalised into two models. The first predicts blood flow through the tumour vasculature and the second predicts interstitial fluid flow throughout the cancerous tissue through use of non-singular Green's functions. Our method enables application to whole, large vascular networks (> 2 cm 3 with > 10 5 blood vessels), thereby permitting predictions of whole tumour fluid dynamics which incorporate the inherent architectural heterogeneities occurring at the micron-scale.
The intravascular component incorporates the model of Pries et al. [55] to simulate vascular blood flow, where the structural properties of the segmented tumour networks and haemodynamic parameters are used as inputs. Flow or pressure boundary conditions at all terminal nodes in the vascular network are required to predict blood flow throughout the network. These boundary data are very challenging to measure in vivo, so we deploy the flow estimation algorithm of Fry et al. [41] to estimate boundary data based on the assumption that the microcirculation is regulated in response to haemodynamics stimuli relating to flow and shear stresses [56]. The scheme estimates unknown boundary conditions by minimising the squared deviation from specified target network wall shear stresses and pressures derived from independent information about typical network haemodynamic properties. In essence, the algorithm removes the need to define conditions at all boundary nodes, to one where simulation sensitivity is weighted towards the definition of these two target parameters. This enables physiologically realistic blood pressure and flow distributions to be estimated across an entire vascular network and has been applied to breast tumour [37], colorectal tumours [40], cortex [57], glioma [40] and skeletal muscle [58]. The second component to our computational model describes fluid transport through the porous interstitium using a Darcy model. Here, the vascular flow solution is coupled to Darcy flow via Starling's law which describes fluid transport across the endothelium. The vasculature is represented by a discrete set of points sources of flux where the source strengths are defined by the vascular blood flow solution. A similar approach has been applied to simulate O2 transport across various tissues [42,57,59] and to predict capillary flow in the absence of these network structures [60]. Our generalised approach enables us to explore the affect of vascular architecture heterogeneity on fluid transport within the interstitium for large-scale vascular networks with spatially heterogeneous tissue and vessel wall permeabilities.
The following sections present the mathematics behind our model, followed by a description detailing its application to large tumour networks.
Vascular blood flow. The segmented tumour networks consist of a series of vessel segments connected by nodal junctions or, in the case of boundary nodes, one-segment nodes which form a boundary to the microvascular network (see Fig 1b). We define a positive flow direction from the start node to end node of each vessel segment. Under the assumption of Poiseuille flow and conserving flow at blood vessel junctions, the relationship between nodal pressures, p k and the boundary boundary fluxes Q 0i is given by where N is the set of all nodes, I is the set of all interior nodes and B is the set of all boundary nodes with known boundary conditions. For all interior nodes, conservation of flux at vessel junctions dictates that Q 0i = 0, however, if i is a known boundary node, Q 0i is the inflow (or outflow if negative). Following the notation of Fry et al. [41], the matrix K ik represents network conductance where S is the set of all segments, and if k is the start node of segment j; is the matrix of vessel conductances where l j , d j and μ j denote the length, diameter and effective blood viscosity of segment j, respectively. We apply empirical in vivo blood viscosity laws, which prescribe the effective blood viscosity as a function of vessel diameter and haematocrit, to compute μ j and consequently incorporate non-Newtonian effects in each individual microvessel [61]. Network haematocrit heterogeneity plays an important part in network flow resistance. However, in this study, we set network haematocrit to 0.45 as we do not have sufficient data to parameterise a haematocrit model at this scale. With future availability of appropriate data, the model has the flexibility to incorporate haematocrit heterogeneity [62].
In the absence of measured flow and pressure data at network boundaries, further assumptions are required to obtain a unique solution. The method proposed by Fry et al. [41] sought to solve a constrained optimisation problem, formulated in terms of a Lagrangian objective function, defined by Here, p 0k is the target pressure at node k, τ j is the wall shear stress in segment j, τ 0j is the corresponding target shear stress, k p and k τ are weighting factors associated with the pressure and shear deviations from the target values, λ i is the Lagrange multiplier associated with node i and w k is the vessel length associated with node k. Setting dL/dp i = 0 and combing with (1) yields a sparse linear system with unknowns p k and λ i . Assigning a pressure drop to a proportion of boundary nodes forms a well-posed system which can be solved using standard methods [41].
The blood flow estimation model by Fry et al. [41] has been thoroughly tested using mesenteric networks [63] in which blood flow measurements were taken in individual vessels and used to inform parameter estimation [41,55,61,64].
Interstitial fluid transport. Darcy's law has been effectively used to describe the passage of fluids [3,14,[30][31][32][33]65] or solutes [42,57] through tissues. In this study, we use Darcy's law to describe the relationship between the IFP, p, and volume-averaged interstitial fluid velocity (IFV), u, within the porous interstitial domain: where κ is the hydraulic conductivity of the interstitial tissue. Here we assume that interstitial pressure tends towards a constant value, p 1 , in the far-field region, Tumours are leaky due to large pores along a vessel's lumen, and so the vasculature exhibits a strong fluid and oncotic interaction. Following the approach of Baxter and Jain [3] and subsequent studies [43,66], we used Starling's law to describe fluid transport across the endothelium: where J v and L p are the fluid flux across and the hydraulic conductance of the vessel wall, respectively, S is the surface area of the vasculature, σ is the oncotic reflection coefficient and, Δp and ΔP are the fluid and oncotic pressure gradients between the vasculature and surrounding tissue. The tumour vascular architecture is used to spatially parameterise the locations of a discrete set of sources of flux into the interstitial domain. Assuming these sources both supply or drain the interstitium, conservation of mass yields where x j and ϕ j are the spatial coordinates and strength, respectively, of point source j, and represents a point source of fluid flux from the vasculature to the surrounding interstitial domain. Applying the substitution � p ¼ p À p 1 , the Green's function, G(x, x � ), for the adjoint problem for � p is given by For a given source, we distribute the delta function over a sphere of finite radius, r 0 . We assume continuity of G and rG at the interface r = r 0 , which enables the following radially symmetric Green's function to be derived: We distribute a total of N s sources across the vasculature with a maximum spacing of δ, and assign a source radii, r 0i , equal to the corresponding blood vessel radius (see Fig 1c). Using Green's superposition principle for linear operators, the convolution of G provides the corresponding pressure solution for source i 2 N s and so where q s j is the vector of source strengths. Here, G ij is the Green's function associated with (10) and defined by where r ij is the distance between sources, defined as r ij = |x i − x j |, and κ i is the interstitial hydraulic conductivity at the location of source i. From (6), the Green's function, G ij , can be used to calculate the volume-averaged IFV, given by where Starling's law, (8), can be rearranged into the form where p v,i (calculated by (12)) and P v,i are the blood and oncotic plasma pressure at the vessel wall, p b,i (calculated by (5)) and P b,i are the vascular blood pressure at source i, in the absence of diffusive interstitial fluid transfer, and oncotic fluid pressure, J v,i is the rate of fluid flow per unit volume from blood vessel i to the interstitium σ i is the oncotic reflection coefficient for vessel i. The intravascular resistance to fluid transport across blood vessel i, is defined by where S i and L p,i is the surface area and vascular conductance of vessel i, respectively. Integrating over the volume of vessel i and assuming flux at the interface is continuous, we define where r v,i is the vessel radius and l i is the length of vessel i. Eqs (12), (16) and (18) are then combined to give a set of N s equations to be solved for the source fluxes q s j , Prescribing parameter values (see Table 2), the resulting solutions for q s j can be solved and used to update vascular pressures, p i,b , using Starling's law, (8), in the absence of an oncotic pressure gradient (i.e. ΔP = 0), and (18). Here, for iteration k + 1, p i,b is set equal to the IFP at wall of blood vessel i, p i,v , calculated on iteration k. This iterative system is repeated by updating q s j values until tolerances are reached (Oð10 À 6 Þ μl/min), subsequently, tissue IFP and IFV fields can be computed using (12) and (14), respectively.
In effect, our computational framework enables a detailed, quantitative assessment of blood and interstitial flow for tissues, both healthy and pathological, where their entire vascular networks are characterised, by encapsulating the flow of fluid between the vascular and interstitial domains. Similar to a Green's function model for O2 transport [42], our model does not require the imposition of explicit boundary conditions on the outer surface of the tissue domain, with the only unknowns to the system being the strength of the set of fluid sources and sinks. As such, when compared to finite difference or element methods, our approach minimises boundary condition artefacts and saves on computational expense as the solutions to the entire mesh are not required. An outline of the interaction between the biomedical imaging, and vascular and interstitial flow compartments is given in Fig 1d. The computational framework was coded using C++ [67] and run on a Apple Mac Pro, with 2 × 3.06 GHz 6-Core Intel Xeon processor and 64 GB of RAM. The system (19) was used solved using a biconjugate gradient method [68] and implemented using the Armadillo sparse linear algebra library [69].

Implementation of computational model on whole tumour microvascular networks
It remains practically infeasible to measure vascular flows and pressures in individual microvessels in vivo, which necessitates a pragmatic approach to boundary condition assignment. Under the assumption that vessels along the tumour surface are connected to peritumoural vessels [37,70], we developed an optimisation procedure which assigns vascular pressures to tumour surface vessels, based on a target pressure drop, with iterative adjustments to match in vivo measurements of mean perfusion from ASL-MRI (see Fig 2a and 2b). These in vivo data are acquired for the same tumours that were subsequently subjected to OPT analysis. Using this approach, we are able to ensure good agreement between in silico predictions and measured perfusion data [40]. In this study a vascular pressure of 30 or 20 mmHg for the GL261 and 45 or 15 mmHg for the LS147T tumour was randomly assigned to 5% of surface boundary nodes, in order to meet the required tissue perfusion. To ensure randomness, the peritumoural nodes were represented by a list. The elements in the list were rearranged randomly using a uniform random number generator where the system clock was used to seed the random number engine. The nodes located in the top 5% of the list were then randomly assigned a low or high pressure using an equivalent randomised approach.
During preliminary simulations it was found that if high/low pressures were prescribed in close local proximity to each other, unphysiological flows were predicted due to the steep local vascular pressure gradient. In order to prevent this, a subroutine was designed so that values at opposing ends of the pressure spectrum could not be located within a defined vicinity of each other. This "exclusion region", centred on a given boundary node, was defined as an ellipsoid volume with diameters, along its three principal axes, equal to 5% of the tissue dimensions.
A proportion of boundary nodes are contained within the tumour volume as a consequence of either angiogenesis which form blind-end vessels, or artefacts of OPT. Previous studies approximated the fraction of blind-end vessels in sample MCa-IV carcinomas to be 33% [71]. However, blind-end information is not available for either GL261 or LS147T tumours. Therefore, consistent with previous computational studies [37], blind-ends were randomly applied to 33% of remaining boundary nodes (using the previous randomised approach), with the remaining 62% of boundary nodes left as unknown in the flow optimisation scheme of Fry et al. [41].
Application of our boundary assignment method requires us to accurately compute mean perfusion across the tumours for a comparison against equivalent experimental data gathered in vivo using ASL-MRI. This requires an accurate definition of the tumour surface and volume to give an accurate approximation of a tumour's mass and fluid flow into the tumour volume. For example, an overestimation of the tumour mass, assuming a cuboid tissue volume surrounding the tumour, can drastically underestimate tumour perfusion since, in this case, the tumour shapes are approximately ellipsoidal and perfusion is inversely proportional to the tumour mass. Similarly, an overestimation of flow into the tissue would overestimate tissue perfusion. Next we describe: 1) defining the surface of the tumour; 2) computing the IFV vectors across the tumour surface; and 3) approximating the total tissue perfusion.
1) The hull of a tumour is calculated using the Matlab (MathWorks Inc., Natick, MA) 'boundary' function applied to all nodes defined during vascular segmentation (see Fig 2c). The Matlab 'fast loop mesh subdivision' triangulation algorithm is then applied to further discretise the define tumour surface. 2) To approximate tumour perfusion requires us to define a set of normal vectors to the tumour surface to compute pressure gradients. We identify the centre of the tumour and duplicate the hull, which is then expanded to form a 10 μm gap between the two surfaces (see Fig 2d). IFP is then calculated across all nodal points on each surface, each paired by a vector normal to the opposing surface. A pressure gradient is then computed along each normal vector and the corresponding velocities are calculated. 3) A sphere packing algorithm is applied to the nodes on the original tumour hull, whereby no sphere overlaps neighbouring spheres (see Fig 2e). Any inflowing node (defined by the corresponding pressure gradient) is averaged across the great circle of its corresponding sphere and its contribution is summed together to calculate the interstitial component of tissue perfusion. Total tissue perfusion is calculated by summing over the peritumoural vascular inlets and interstitial perfusion values.
Finally, we prescribe baseline parameter values (see Table 2). We assume that the tumours were isolated in subcutaneous tissue in the absence of lymphatics, therefore the far-field pressure, p 1 , was set to 0 mmHg. Due to a lack of experimental data, the vascular conductance and interstitial conductivity, L p and κ, respectively, were given uniform values based on literature (2.8 × 10 −7 cm mmHg −1 s −1 and 1.7 × 10 −7 cm 2 mmHg −1 s −1 , respectively [31,75]). As the transport of blood plasma proteins is not modelled explicitly in our model, we assume a constant oncotic pressure gradient of 5 mmHg between the vasculature and interstitium, and set the oncotic reflection coefficient to a uniform value of 0.82 [72].

Results
In the following section we apply our computational framework to a GL261 orthotopic murine glioma and a LS147T human colorectal carcinoma xenograft to form baseline flow solutions (see Figs 3 and 4). Next, we explore sensitivity to source parameters, which include source distribution, source radii and bilateral communication between the vascular and interstitial compartments. We then perform sensitivity analysis to the interstitial parameters (for example, vascular hydraulic conductance and interstitial hydraulic conductivity) on IFP and IFV profiles in the LS147T tumour. This is followed by in silico experiments to establish the tumour parameters, such as vascular and interstitial permeability, which change tumour IFP and perfusion following vascular normalization therapy, and show how vascular architectural heterogeneities affect treatment response.

Simulating interstitial fluid pressure in real-world cancerous tissue
Our vascular flow simulations are in good agreement with those in computational [37,76] and experimental literature (see Table 3 and Fig 3). To test the variability induced by our stochastic boundary condition implementation, the optimisation procedure (detailed in Fig 2a and 2b) was repeated for a total of n = 12 for each tumour. Interstitial fluid flow was then simulated for each separate vascular flow solution, providing us with a baseline set of interstitial flow solutions. Blood flow across all simulations exhibited similar spatial distributions, with perfused vessels mainly restricted to the outer rim of the tumours [77]. The mean standard deviation of vascular blood pressures across all simulations was low with values of 0.82 and 1.53 mmHg (with maxima of 4.78 and 13.2 mmHg) for GL261 and LS147T, respectively (see Fig 5a and 5e). The elevated standard deviations in vascular pressure were located at the periphery of the tumours, which is to be expected as the high and low vascular pressures were stochastically assigned here. Furthermore, our mean blood velocity and vessel wall shear stresses agree with  [37].
Tissue perfusion (calculated as 83.7 ± 22.3 and 5.4 ± 0.3 ml/min/100g for the GL261 and LS147T tumours, respectively) was further validated by in vivo ASL-MRI measurements of 110 ± 70 and 19 ± 8 ml/min/100g for GL261 and LS147T, respectively. This implies that imposing physiologically realistic pressure boundary conditions generates physiologically realistic perfusion, and consequently accurate drug delivery solutions [40].
No literature IFP values were available for GL261 cell lines, however, our IFP was slightly elevated compared to that previously measured in LS147T in vivo (13.5 ± 11.3 mmHg [74]). Considering the range of IFP both here and in vivo [74] and the good accordance with in vivo perfusion, our results provided us with the confidence that our simulations can produce physiological IFP predictions (see Fig 4).  Examining the tumour radial IFP profiles, the LS147T network exhibited similar configurations as observed both experimentally in LS147T [17], in other cell lines [3,24,25], and in computational studies [31,32,36,75], with elevated IFP at the tumour core (see Fig 4b and  4d). In addition, the LS147T network displayed a typical IFV profile radially, with an increasing IFV range towards the tumour surface due to the steeper pressure gradients at the periphery of the tumour. This indicated that bulk fluid filtration occurs at the high flowing vasculature located at the tumour extremity in this network (ρ = 0.41, p < 0.001, where ρ and p are the Pearson's correlation coefficient and its corresponding p-value between tumour radius and vascular flow, respectively).
In comparison, the GL261 network also exhibited a traditional, yet steeper, IFP profile with a wider variance throughout the tumour (see Fig 4c). The IFV profile exhibited a typical profile increasing from the centre of the tumour to the periphery. However, its IFV peak was reached at � 80% of the tumour radius, with a substantial decline in the latter 20% (see Fig 4c).
High IFP has been associated with low vascular density in A-07-GFP tumours [78]. Similarly here, we observed that GL261 and LS147T have distinct differences in their vascular architecture (see Fig 3 and Table 1) and that, in the case of LS147T, low variability in vascular density was associated with higher IFP (ρ = −0.584, p < 0.001-including in silico predictions from [40]). This may indicate that the inherent vascular architectural heterogeneity across tumour cell lines [20] directly impacts the IFP and IFV distributions, creating an unorthodox interstitial flow profile.

Sensitivity of tumour interstitial fluid transport to model parameters
Here, we perform sensitivity analysis to understand the impact of parameter variance on predictions of tumour fluid dynamics. In doing so, not only do we understand the sensitivity to the model, but gain insight into the biological mechanisms at play.
For convenience we have split these parameters into two groups. The first we define as source parameters, which include the vascular blood pressure, p b,i and the source radius, r 0,i for segment i, and the maximal spacing between each source along a given vessel, δ (see Fig 1c). The second group we call the interstitial parameters, which describe the biomechanical mechanisms which affect tissue transport. These include the oncotic reflection coefficient, σ, the hydraulic conductance of a vessel wall, L p , the interstitial hydraulic conductivity, κ and the farfield IFP, p 1 . In the following, if sensitivity analysis to an interstitial parameter was not being performed, it was set to the corresponding value in Table 2.
Source parameters. Blood pressure distributions across a vascular network depend not only on the vascular pressure at the network boundaries and network geometry, but also transvascular fluid transport from the tumour vasculature and into the surrounding interstitium, and vice versa. To incorporate flow communication between these two domains, we couple the vascular and interstitial models using an iterative scheme in which vascular blood pressure distributions, p b,i for i 2 N s , were updated on each iteration by incorporating the loss of fluid within a vessel due to fluid flux across the vessel wall into the interstitium (see Fig 1d). We found that due to the small volume of fluid leaving the vessels (relative to the volume of fluid flowing through the vessels), vascular pressures did not vary significantly between iterations (see Fig 5b and 5f). In the case of the GL261 network, the algorithm had converged after two iterations (see Fig 5b). In comparison, the LS147T network had not converged after six iterations but the mean nodal pressure error was Oð10 À 8 Þ mmHg (see Fig 5f). This may indicate that a tumours inherent vascular heterogeneity and corresponding interstitial parameters alter the scale of these flow errors.
Next we varied the maximum distance, Δ max , between sources distributed along the vasculature. Values of 10 and 25 μm, and 50 and 100 μm were chosen for the GL261 and LS147T networks, respectively, due to differences in mean branching vessel lengths observed in each tumour (see Table 1). The GL261 network exhibited minimal sensitivity in IFP between the two maximal source lengths (see Fig 5c). Comparatively, the LS147T experienced greater variability (of less than 1 mmHg) in mean IFP, between the two values of Δ, at the core and periphery of the tumour (see Fig 5g). We hypothesise that sensitivity to Δ is related to the vascular density of each tumour type. The GL261 network is an order of magnitude greater in vascular density compared to the LS147T tumour (see Table 1). An elevated vascular density results in an increased density of fluid sources, and so any reduction in sources may have a minimal effect to IFP due to flow being maintained by local neighbouring sources.
We finally sought to understand how the source radii, r 0,i for i 2 N s , affect flow in the interstitial domain. We initially set the source radii to a constant value equal to the minimum vessel radius in a given tumour network. Three cases were then explored where r 0,i was multiplied by a factor of 10 −1 , 10 0 or 10 1 . Our results show that decreasing the value of the source radii decreases mean IFP, with greater sensitivity exhibited again by the LS147T network (see Fig 5d  and 5h), again indicating that vascular architectural heterogeneities can illicit different responses in IFP between tumours. This indicates that greater care is needed compared to the assignment of other source parameters, to ensure the physiological accuracy of the results. Here, we set the value of r 0,i equal to the radius of its corresponding vessel in order to provide a heterogeneous distribution of radii in which vessels with a larger radius are able to influence the IFP to a greater capacity than compared to relatively smaller vessels.
Interstitial parameters. Sensitivity to the interstitial parameters p 1 , σ, L p , and κ were explored using the LS147T network. This network was chosen due to the scale of simulations to be performed and due to the greater sensitivity exhibited by modifying the source parameters. First of all, we investigated the deviation in the far-field interstitial pressure, p 1 . Ideally, the far-field pressure is set to the mean IFP of the adjoining healthy tissue, however these data are frequently unavailable. We varied the far-field IFP from 3 to 18 mmHg in increments of 3 mmHg. Alteration of p 1 did not significantly modify tumour IFP, with the majority of variation occurring at the periphery of the tumour (see Fig 6a). As a consequence, mean IFV across the tumour decreased with increasing p 1 (see Fig 6e).
Next we investigated the oncotic reflection coefficient, σ. Physiologically, the coefficient σ varies between 0 and 1, indicating the likelihood that a molecule approaching a pore in the vessel lumen will be reflected back and thereby retained in the vascular compartment. As such, we ranged σ uniformly from 0 to 1 with increments of 0.25. Modification of σ resulted in a network-wide shift in tumour IFP magnitude, thereby preserving the IFP gradient (see Fig 6b). This is an intuitive result since we held the oncotic pressure gradient constant across the entire vasculature, therefore any increase in σ resulted in systematic reduction in both IFP and IFV (see Fig 6b and 6f).
Vascular hydraulic conductance, L p , defines the leakiness of the lumen to the transport of blood plasma. Our sensitivity analysis varied L p uniformly across the tumour, from 1.33 × 10 −8 to 3.41 × 10 −6 cm/mmHg s, encompassing values provided in literature for normal and tumour tissue [11]. Similarly to σ, a reduction in L p resulted in decreasing IFP across the tumour but with contrasting regional gradients (see Fig 6c). Our analysis showed that with decreasing L p , IFP tended towards the assigned far-field pressure due to decreasing fluid filtration from the vasculature, and so in its limit, the IFP would uniformly be equal to p 1 . This is due to decreasing fluid filtration, indicating that the model produces a physiologically viable response. In comparison, a reduction in L p induced minimal changes to IFV in the core of the tumour but resulted in decreases towards the periphery (see Fig 6g).
Interstitial hydraulic conductivity, κ, encapsulates the structure of the ECM, such as stiffness and density, and determines the permeability of the interstitium to fluid transport. Here, we ranged κ from 1.33 × 10 −11 to 1.04 × 10 −6 cm 2 /mmHg s uniformly through the tumour. In this case, we established that sensitivity to this parameter are non-trivial to elucidate. An inflexion point was observed for values 1.67 × 10 −9 � κ � 8.33 × 10 −9 cm 2 / mmHg s in which the spatial IFP distribution switched from one configuration to another (see Fig 6d). For κ � 1.67 × 10 −9 cm 2 / mmHg s, IFP predictions displayed near identical configurations (mean IFP of 25.8 ± 2.4 mmHg for all four cases) in which IFP displayed an increased gradient compared to greater values of κ along with decreasing IFV (see Fig 6h). In comparison, κ � 8.33 × 10 −9 cm 2 / mmHg s displayed similar IFP profiles with mean IFP decreasing with increasing κ (26.2 ± 1.5 to 22.5 ± 1.6 mmHg).
In the next section we apply our computational model to predict changes in tumour fluid dynamics following vascular normalization therapy.

Normalizing blood vessel lumen minimally reduces IFP
Vascular normalization is a therapeutic strategy used to restore tumour vasculature to a structural and functional state exhibited by healthy blood vessels [13,15,18,49,51,79]. The impact of vascular normalization agents include a reduction in microvessel diameter, pruning of immature vessels, increased vascular maturity and pericyte coverage, and reduced vessel tortuosity [7,51,80]. As a result, normalization has been shown experimentally to lower tumour IFP [51,52], enhance tumour vascular perfusion [52,53] and improve drug delivery [48,49,81,82]. However, the mechanistic links are missing which necessitates use of in silico modelling.
We next use our computational framework to investigate the effectiveness of vascular normalization in both GL261 and LS147T tumours where normalization therapy is administered intraperitoneally [48,82]. We assume an ideal case in which vascular normalization has occurred in an axisymmetric manner whereby the effectiveness of the treatment increases from the tumour core to its periphery. To achieve this, the extent of stabilisation of vessel diameters and normalization of lumen permeability (through a combination of changes in L p and σ), as a result of increase pericyte coverage, ranges linearly from typical tumour values in the core to physiological levels at the periphery. Here, vessel diameters were decreased by a factor of 1.99 [79], and set L p and σ to 0.44 × 10 −7 cm mmHg −1 s −1 and 0.91 [83], respectively.
Similarly to experimental [51,52] and computational [81,84] studies, vascular normalization reduced vascular perfusion (see Fig 7) and mean IFP (see S1 Fig) across both the GL261 and LS147T tumours. However, our predictions showed that IFP remained elevated at the core and that the overall reduction was induced by an increase in the pressure gradient across the tumours. This resulted in a small increase in IFV at the core of both tumours (see S1 Fig), yet changes in interstitial perfusion were minimal.
Our therapeutic predictions of vascular normalization are inconsistent with the magnitude of observed effects, such as reduction in IFP, compared to experimentation [52,53]. This suggests that normalizing the permeability of the vessel lumen and stabilising vessel diameters are not significant therapeutic affects following vascular normalization. We advocate that vascular normalization reduces tumour IFP and increases perfusion via other factors, such as changes in interstitial permeability or a reduction in vascular network tortuosity.

Vascular normalization locally alters interstitial hydraulic conductivity
Vascular normalization fortifies blood vessels through an increase coverage of mural cells and basement membrane [85]. Considering this and our sensitivity analysis for κ (see Fig 6d and 6h), we hypothesise that normalization locally alters the interstitial hydraulic conductivity of a treated blood vessel. We therefore sought to normalize interstitial hydraulic conductivity, in parallel and in the same manner as with vessel diameters, L p and σ, to a value of 8.53 × 10 −9 cm 2 mmHg −1 s −1 [86] at the tumour periphery.
Normalizing interstitial hydraulic conductivity had substantial and contrasting consequences to fluid transport in the interstitium of both GL261 and LS147T tumours. In the case of GL261, IFP was elevated compared to baseline simulations throughout the tumour (see Figs 8a and 9a), with a mean increase of 70%. This, in conjunction with a small increase in the IFP gradient, led to elevated IFV in 80% of the tumour (see Figs 8a and 9a) and consistent patterns in tumour perfusion compared to baseline (see Figs 4a and 9b).
Treating LS147T significantly modified the IFP profile (see Fig 9d), elevating IFP by � 9 mmHg at its core and developing a steeper gradient towards the periphery, in addition to a narrowing of IFP standard deviation (see Fig 8b). As a result, IFV increased throughout the tumour core (see Fig 9d), with an increase of 50% at its centre, and a predicted decreased compared to baseline in the outer 50% of the tumour (see Fig 8b).
To understand the response to treatment we analysed the source density, the ratio between sources and sinks of fluid flux, across each tumour type pre-and post-normalization. Fig 9c shows that GL261 is dominated by sources of flux pre-treatment, with the number of sources elevating at its core post-treatment. In comparison, a parity is observed between sources and sinks in the core of LS147T pre-treatment. Following normalization, the balance between sources and sinks drastically altered with no sinks existing between 0 to 0.3 of the LS147T radius, this is followed by a decline in the ratio to levels lower that baseline at the periphery, indicating an increase of sinks here.
Our results suggest that observed fluid dynamic changes following vascular normalization in literature [51,52,52,53] may be as a result of local modification of the interstitial hydraulic conductivity (see S2 Fig). Further, that changes in perfusion may be a result of an steepening of the IFP gradient across a tumour and not necessarily a uniform reduction in IFP across the tumour. This suggests the importance of incorporating inherent spatial heterogeneities in vascular networks observed across tumour cell-lines [20] (see Table 1

Discussion
Elevated interstitial fluid pressure is frequently associated with solid tumours, where a conventional profile exhibits a uniformly high IFP in the core of a tumour decreases rapidly towards the levels of physiological tissue at the periphery [3,14,65]. This atypical characteristic forms a barrier to transvascular fluid and drug delivery, thereby diminishing therapeutic efficacy of anti-cancer treatments.
The passage of fluid through the interstitium is influenced by both hydrostatic and oncotic pressures in blood vessels and therefore by the heterogeneous architecture of tumour microvessels. However, the procurement of detailed fluid flow data in vivo across whole tumour networks is currently infeasible using conventional imaging and experimental techniques. However, recent advances in ex vivo high-resolution optical imaging techniques [39] allow whole three-dimensional tissue architectures to be extracted and reconstructed to act as inputs for detailed in silico modelling of fluid transport.
No study has explicitly modelled vascular and interstitial fluid flow using discrete, threedimensional structural data from images of whole-tumour, real-world vasculature [6]. However, in a recent study, we presented our novel REANIMATE platform which extracts three-dimensional, whole tumour vasculature ex vivo from optically cleared tissue, which is then used to parameterise an in silico model of fluid transport guided by in vivo imaging data [40]. This has enabled us to perform quantitative, realistic predictions of fluid and drug delivery to tumours which has led to novel insights into a tumour's inherent physical resistance to anticancer therapies [40].
We present a generalised framework to model vascular and interstitial fluid dynamics based on whole, explicit tumour vasculature, in a computationally tractable way. We detail its derivation and application to whole tumour vascular datasets and show how our model allows highly-detailed predictions of fluid flow within the tumour microenvironment by incorporating explicit tumour vasculature and spatially heterogeneous parameters, such as vascular and interstitial permeability. Our model allows flow heterogeneities to be quantified in a computationally efficient manner, when compared to finite-difference and element methods [42], and consequently can be applied to any vascular tissue.
We initially apply our framework to an orthotopic murine glioma and a human colorectal carcinoma xenograft from the GL261 and LS147T cell-lines, respectively, to present realistic, baseline simulations of the tumour microenvironment. We then perform sensitivity analysis to the underlying model parameters. The first are the source parameters, which are specific to our model, and include point source distribution and size. Secondly, the interstitial parameters, such as vascular conductance and interstitial conductivity, which are frequently represented in literature due to the prominent use of Starling's law. In the second case, we perform analysis of how variation of these parameters modifies the IFP and IFV in an LS147T dataset. We finally mimic vascular normalization therapy, via varying model parameters, to generate hypotheses relating to fluid dynamic response observed in literature.
Our computational framework is based upon a Poiseuille model for vascular blood flow [41] which is coupled to a steady-state Green's function solution to interstitial fluid flow. Here, tumour vasculature segmented ex vivo is represented by a discrete set of sources of fluid flux for bi-directional transport between the interstitium. Previous models either assume a homogeneous vascular network [3,11,14,65], incorporate a computer-generated synthetic tumour network [28][29][30][31][32][33][34], or incorporate boundary conditions and spatial variations in tissue permeability to artificially represent vascular heterogeneity [35,36]. However, vascular averaging methods do not fully encapsulate the intrinsic, local interactions between neighbouring blood vessels which contribute to global interstitial flow and synthetic networks are difficult to validate against real tumour architecture. Here we use vascular architecture from real, whole tumour networks, and through use of Green's function methods, our model significantly reduces the computational size of the computational problem, allowing vessel-vessel interactions to be modelled at the micron-scale. Thus, we provide the means to perform in silico studies to hypothesis test the impact of vascular heterogeneity on the tumour microenvironment with relative ease.
Our simulations were performed on ex vivo structural imaging data from a GL261 and a LS147T tumour. As no in vivo flow or pressure data were available for the numerous boundary vessels, we developed a procedure whereby simulated data is optimised based on in vivo tissue perfusion data gathered using ASL-MRI [40]. This approach has produced solutions which are highly consistent with experimental measurements in the same tumours [40]. In this study, baseline vascular flow solutions across all tumour simulations are in good agreement with the perfusion data, alongside mean flows [76], velocities and vessel wall shear stresses [37], and fluid pressure in the interstitium [3,17,24,25]. This provided validation that our model produces physiologically realistic results, providing a platform to investigate the tumour microenvironment.
We performed sensitivity analysis to the source parameters, such as updating the vascular flow solution, source distribution and source radii, to understand their influence on the flow communication between the vascular and interstitial domains. Our results exhibited a minimal sensitivity to IFP distributions by varying these parameters, with the exception of the assignment of source radii. Here, we hypothesise that greater care is required for spatially sparse tumours with a low vascular density, in order to ensure physiologically accurate simulations. Investigating the sensitivity to the interstitial parameters in the LS147T tumour, we found that raising the far-field interstitial pressure did not significantly alter the IFP distribution across the tumour, with increased IFP occurring at the tumour periphery, compared to baseline. Uniform variation in the oncotic pressure contribution, by modifying the oncotic reflection coefficient, only affected the magnitude of the interstitial pressure, with minor changes to the IFP and IFV profiles, similar to those reported in previous studies [2]. Increasing L p and κ, raised the gradient of the interstitial pressure profile, thereby increasing fluid transport through the interstitium. However, our IFP distributions did not reach 0 mmHg at the periphery of the tumours as in previous studies [3,14,84]. This is due to not applying a fixed pressure at the tumour boundary in our model, which results in a smoother transition of IFP to the surrounding tissue, similar to previous computational modelling of tumour vascular heterogeneity [35,36].
We next predicted the flow response following vascular normalization therapy mimicked by varying model parameters. Normalization of tumour vasculature was achieved by modifying vascular hydraulic conductances, oncotic reflection coefficients and stabilising blood vessel diameters in a linear radially varying fashion in both tumour networks. Our results indicate that whilst normalisation of lumen permeability and vessel diameters exhibited similar trends compared to experimentation such as a reduction in vascular perfusion post-treatment [52,53], the magnitude in IFP reduction was minimal. Similarly, therapy which soley reduces vessel leakiness has been shown not to be effective in silico using synthetic tumour vasculatures [34].
As vascular normalization has been shown to stabilise the tumour vasculature by increasing mural cell and basement membrane coverage [85], we hypothesise that these structural changes alter the interstitial hydraulic conductivity. Subsequent simulations observed significant and unique changes in IFP, IFV and perfusion across each tumour cell-line. In GL261, we predicted an increase in IFP across the entire tumour, whereas in LS1247T, IFP exhibited a steep pressure gradient and elevated IFP in its core. Consequently, IFV was elevated in both tumour cores. We noted that this significantly altered the interstitial perfusion map in LS147T, with elevated perfusion at the tumour core, which dissipated towards its periphery.
We have shown that incorporating realistic tumour vasculature is key to accurately predicting the spatial flow heterogeneities induced by blood vessels. Our results suggest that reducing tumour IFP is not necessary to increase tumour interstitial perfusion and that developing a IFP gradient across the tumour is the key response. Further, we hypothesise that vascular normalization alters interstitial hydraulic conductivity within tumours. However, as our tumour vascular networks are static, we cannot discern the relative impact caused by a reduction in vascular tortuosity observed experimentally following therapy [51,80]. Notwithstanding, we show that normalizing the tumour interstitial environment may provide an effective means to increase tumour perfusion.
Our results and proposed computational framework offer significant scope for future expansion. For example, recent in vivo methods provide a step forward in approximating conditions by measuring interstitial fluid velocity at the macro-scale [26]. These data could lead to greater accuracy when assigning boundary conditions specific to a tumour. Further, parameter values such as the vascular conductance and interstitial hydraulic conductivity were assigned using previous literature values since these tissue-specific measurements can be challenging to procure through experimentation. For example, interstitial conductivity values across healthy tissue have been reported to span four orders of magnitude [17]. New methods need to be developed to accurately quantify these values.
There are also opportunities to expand the computational model to incorporate more complex biological phenomena. For example, incorporating tumour compression of vessels due to increasing shear stresses within the tumour [29,87], volumetric tissue growth and tumour angiogenesis which would allow us to develop vascular normalization treatment strategies which determine when to administer therapy during the normalization window [12].
We expect to find a wide utility for REANIMATE in a range of disease areas, particularly given the current interest in optical clearing methods and their widespread use in biomedical research. Our in silico framework is novel and timely and will find extensive use for hypothesis testing, to enable tumour biology and drug delivery to be better understood, which in turn may enable the next generation of cancer therapies.
Supporting information S1 Fig. Simulated normlization of L p , σ and vessel diameters in (a) GL261 and (b)