article: 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.
  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–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–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–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–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–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 cm3 with > 105 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–45]. By modifying vascular architecture, therapy aims to normalize tumour perfusion and oxygenation, thereby increasing the efficacy of chemo, radio and immunotherapy [44,46–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–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.
  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 perfuse-fixation, 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 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]. 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 cm3 with > 105 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. 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∞, was set to 0 mmHg. Due to a lack of experimental data, the vascular conductance and interstitial conductivity, Lp and κ, respectively, were given uniform values based on literature (2. 8 × 10−7 cm mmHg−1 s−1 and 1. 7 × 10−7 cm2 mmHg−1s−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].
  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 similar numerical modelling of vascular blood flow in the MDA-MB-231 breast cancer cell line [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. 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, pb, i and the source radius, r0, 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, Lp, the interstitial hydraulic conductivity, κ and the far-field IFP, p∞. In the following, if sensitivity analysis to an interstitial parameter was not being performed, it was set to the corresponding value in Table 2. 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 Lp 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 Lp 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 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, Lp and σ, to a value of 8. 53 × 10−9 cm2 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 and S3 Fig) into in silico studies of vascular normalization.
  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, three-dimensional 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 anti-cancer 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–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 Lp 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.
summary: The structure of tumours varies widely, with dense and chaotically-formed networks of blood vessels that differ between each individual tumour and even between different regions of the same tumour. This atypical environment can inhibit the delivery of anti-cancer therapies. Computational tools are urgently required which facilitate a deeper understanding of the relationship between blood vessel architectures and therapeutic response. We have developed a computational framework which integrates the complex tumour vascular architecture to predict fluid transport across all lengths scales in whole tumours. We apply our model to two tumour cell-lines and show that differences in their inherent vascular structures influence flow through cancerous tissue. We also use our platform to predict the fluid dynamic response following vascular normalization therapy in realistic, static tumour networks and show that the response is dependent on tumour vascular architecture. We hypothesise that therapy may alter the permeability of interstitial tissue to fluid transport and show that lowering interstitial fluid pressure is not a necessary therapeutic outcome to increase tumour perfusion.
keywords: medicine and health sciences
  body fluids
  classical mechanics
  fluid mechanics
  cardiovascular physiology
  cardiovascular anatomy
  cancer treatment
  basic cancer research
  angiogenesis
  oncology
  developmental biology
  tumor physiology
  tumor angiogenesis
  blood vessels
  vascular permeability
  fluid dynamics
  blood pressure
  continuum mechanics
  fluid flow
  physics
  blood flow
  interstitial fluid
  blood
  anatomy
  physiology
  biology and life sciences
  physical sciences
  vascular medicine
title: Modelling the transport of fluid through heterogeneous, whole tumours in silico