Moath Alamer，Xiao Yun Xu

1 Department of Chemical Engineering Imperial College London, South Kensington Campus, London, United Kingdom

Abstract Tumour vasculature is known to be aberrant, tortuous and erratic which can have significant implications for fluid flow. Fluid dynamics in tumour tissue plays an important part in tumour growth, metastasis and the delivery of therapeutics. Mathematical models are increasingly employed to elucidate the complex interplay between various aspects of the tumour vasculature and fluid flow. Previous models usually assume a uniformly distributed vasculature without explicitly describing its architecture or incorporate the distribution of vasculature without accounting for real geometric features of the network. In this study, an integrated computational model is developed by resolving fluid flow at the single capillary level across the whole tumour vascular network. It consists of an angiogenesis model and a fluid flow model which resolves flow as a function of the explicit vasculature by coupling intravascular flow and interstitial flow in tumour tissue. The integrated model has been used to examine the influence of microvascular distribution, necrosis and vessel pruning on fluid flow, as well as the effect of heterogeneous vessel permeability. Our results reveal the level of nonuniformity in tumour interstitial fluid pressure (IFP), with large variations in IFP profile between necrotic and non-necrotic tumours. Changes in microscopic features of the vascular network can significantly influence fluid flow in the tumour where removal of vessel blind ends has been found to reduce IFP and promote interstitial fluid flow. Our results demonstrate the importance of incorporating microscopic properties of the tumour vasculature and intravascular flow when predicting fluid flow in tumour tissue.

Keywords Fluid flow, Vasculature, Tumour, Mathematical modelling

INTRODUCTION

Fluid flow in tissue plays an important role in the delivery of oxygen, nutrients and therapeutic agents.The physiology of tumour tissues exhibits several distinct features that influence fluid flow which can potentially promote tumour invasiveness and limit drug delivery. As tumours grow beyond the limit at which preexisting vasculature is able to sufficiently deliver blood, hypoxia develops which can trigger angiogenesis (Carmeliet and Jain 2000; Liao and Johnson 2007). This leads to the formation of aberrant vascular networks with non-uniformly distributed blood vessels that lack hierarchy and differentiation.Excessive branching can be found with loops, blind ends, arteriovenous shunts and erratic changes in diameter which can result in impaired blood perfusion(Less et al. 1991). Cellular abnormalities of blood vessels lead to poor vessel stability and large interendothelial gaps that allow for excessive leakage of plasma (Baluk et al. 2005). In the tumour extravascular space, the extracellular matrix (ECM) is produced at a high rate, generating mechanical stress which combined with the leaky vessels leads to high interstitial fluid pressure (IFP) (Heldin et al. 2004).These properties act together to cause abnormal fluid flow in tumours which can be heterogenous and vary from patient to patient (Swartz and Lund 2012). Most therapeutics are injected intravenously and need to reach the tumour site, distribute within the vascular network, permeate through the vessel wall and travel across the interstitial space to reach cancer cells (Zhan et al. 2018). Hence, abnormal blood flow, reduced extravasation and pressure gradients caused by high IFP can result in poor accumulation and distribution of therapeutic macromolecules. Additionally, high IFP and poor perfusion have been associated with increased metastatic potential and poor prognosis (Hompland et al. 2014; Munson and Shieh 2014). Therefore, gaining deep insight into fluid dynamics in tumour tissue is necessary as it plays a key role in tumour growth,metastasis and the delivery of therapeutics (Dewhirst and Secomb 2017).

Experimental methods have been applied to examine fluid flow in tumour tissues, however, they are limited by the scale and resolution at which they are able to resolve fluid dynamics (Boucher et al. 1990; Gulliksrud et al. 2009). Mathematical modelling offers a costeffective approach to gaining quantitative understanding of fluid flow in solid tumours at multiple scales. Several mathematical models have been developed to simulate flow in solid tumours (e.g.,Baxter and Jain 1989; Soltani and Chen 2011, 2012;Zhan et al. 2017). These models treated the tumour vasculature as a uniformly distributed source term without accounting for spatial heterogeneities in vascular distribution. However, tumour vasculature is known to be heterogenous which can lead to heterogenous distribution of molecules and liposomes in tumour tissue (Bhandari et al. 2017a; Vavra et al.2004). The influence of tumour vascular heterogeneity on interstitial and transvascular flow has been investigated by several authors (Mohammadi and Chen 2015; Soltani and Chen 2013; Welter and Rieger 2013;Wu et al. 2008). In some of these studies, macroscopic parameters for the whole tumour, such as the surface vascular density term, were predefined based on data in the literature rather than considering the morphology and architecture of the network generated.Other studies applied homogenization methods where properties of the vasculature were averaged over discretized cells. By applying Starlings law using these homogenization techniques, only an approximation of the effect of vascular distribution on fluid flow could be investigated. Furthermore, these models treated intravascular flow and interstitial flow in a decoupled manner, so that the effect of transvascular leakage on blood flow was not captured. Although several recent studies addressed the effect of heterogeneous vasculature by making use of contrast-enhanced magnetic resonance imaging (Bhandari et al. 2017b,2018; Zhan et al. 2014), these models did not explicitly describe the microvasculature and intravascular flow.The high permeability of tumour vasculature results in excessive fluid leakage from the vessels, hence a strong coupling between vascular and interstitial flow is needed.

Several studies have addressed the coupling between vascular and interstitial flow at the single capillary level(Baish et al. 1997; Netti et al. 1996). Pozikridis and Farrow developed a Green’s function based mathematical model where sources and dipoles are distributed along a cylindrical vessel segment describing the distribution of vascular and interstitial pressure over the inner and outer surface of the vessels(Pozrikidis and Farrow 2003). This formulation reduces errors introduced by homogenization methods and allows for the interstitial pressure to be described as a function of the architecture and orientation of the vasculature. The model captures the strong coupling of flows by accounting for the influence of extravasation and IFP on intravascular flow. Pozrikidis and Farrow demonstrated that while neglecting the influence of IFP on intravascular flow might be acceptable when vascular hydraulic conductivity is low as in normal tissue, it could result in significant overestimation of transvascular leakage in tumour tissue where vascular hydraulic conductivity is high. So far, Pozrikidis’s fluid flow model has been applied to a single capillary level or a network of ordered vessels (Pozrikidis 2010).Recently, a similar Green’s function based fluid flow model was applied to a vascular network obtained from optical images of cleared tumour tissues (Sweeney et al.2019), where blood flow was assumed to be conserved at each junction and was not influenced by interstitial flow. In all the models reviewed so far, the permeability of the vasculature was assumed to be uniform across the vascular network which failed to capture the heterogenous pore size distribution in tumour vasculature as it can range from 7–1200 nm (Hobbs et al.1998). In terms of evaluating fluid flow in tumours,there is a need to develop models that can predict flow in an entire tumour with explicit descriptions of individual vessel features, such as radius, length, orientation and permeability, whilst capturing the strong coupling between intravascular and interstitial flow.

In this work a mathematical model is developed by coupling Anderson and Chaplain’s tumour induced angiogenesis model (Anderson and Chaplain 1998) with Pozrikidis’s fluid flow model (Pozrikidis and Farrow 2003). The angiogenesis model is extended to 3D to accommodate unique features of the tumour vascular network including excessive branching, looping and high tortuosity. The model describes the response of sprouting vessel tips to tumour angiogenic factors (TAF)secreted by tumours cells. TAF distribution is initialised to describe different tumour geometries with varying degrees of necrosis centered in the domain surrounded by normal tissue. A discrete method is used to track the movement of tip endothelial cells in response to random motion, TAF gradients and fibronectin distribution in the tissue. This forms a network which is regulated and discretized where vessels are divided into short cylindrical segments with a specific length and radius. A fluid flow model is then coupled to the vascular network where vascular, transvascular and interstitial flow are described by Poiseuille’s, Starlings and Darcy’s law respectively. These equations are solved using Pozrikidis’s method to capture the strong coupling between vascular and interstitial flow. Applying Pozrikidis’s model, the interstitial pressure is integrated over the vessel surface, hence geometric features of the vessel, including length, radius and orientation are considered when calculating the IFP profile in the tumour. We have extended the model from a single capillary scale to a complex large scale vascular network,and applied the model to different tumour cases to understand the effects of various properties such as vascular architecture, distribution, tumour necrosis and microscopic details of the tumour network.Furthermore, we have modified the model to allow for variations in vascular permeability as a function of vessel maturity and flow conditions, pushing the model a step closer towards capturing the complex characteristics of tumour vasculature.

RESULTS

Tumour vascular network

The vascular networks generated using the angiogenesis model were analysed and their morphological parameters were compared with the corresponding properties found in real tumour vasculature. In this work, three models with a normalized radius of 0.25 were simulated with varying degrees of necrosis: a non-necrotic tumour network, a case with a necrotic core accounting for 3% of the total tumour volume and another case with a 26% necrotic core volume. Table 1 shows the values for the morphological parameters (defined in the Methods section) calculated for the three tumour models. The evaluated morphological parameters are within the range reported in the literature, demonstrating that the generated models are representative of real tumourvascular networks. Comparison of the three models suggests that the average vascular density values are similar; this is because as the degree of central necrosis increases, the tips would move away from the core and create more vessels in the peripheral tumour region.

Table 1 Morphological and hemodynamic parameters for different tumour networks with varying degrees of necrosis

After evaluating the vascular networks generated for the different tumour models, the fluid flow model was applied and the calculated hemodynamic parameters were averaged over the entire vascular network for comparison with relevant data in the literature. As shown in Table 1, the mean flow, velocity and shear stress obtained with the fluid flow model are within the reported range, showing the ability of the fluid flow model to capture the intravascular fluid dynamics in tumours.

It is known that abnormal structural and architectural features, including self-loops, vessel compression and blind ends, are often found in tumour vasculature. The presence of blind ends can cause very low or no flow in some vessels. Another feature is the presence of arteriovenous (AV) shunts where short, low-resistance and high-flow pathways form in close proximity to the feeding vessels. A consequence of these AV shunts is that blood preferentially passes through them due to a large pressure gradient, hence blood can completely bypass the entire capillary network. In the vascular networks described in Table 1, some abnormalities of the tumour vasculature can be captured by modifying the non-necrotic tumour vascular network model (Fig. 1).For example, the default vasculature generated in the non-necrotic model does not contain AV shunts. To simulate this feature, vessels in close proximity to the arterial end are set to branch off where one end is connected to the tumour capillary network and the other end of this branch is assigned a venous pressure boundary condition. This allows for tumour vascular networks with and without AV shunts to be generated and the influence of this vascular abnormality to be examined.

Fig. 1 Tumour vascular networks containing abnormal features. A Self-loops. B Vessel compression. C Blind ends. D Arteriovenous shunts

Fluid flow in tumour tissue

Effect of vascular distribution and necrosis

Intravascular and interstitial pressure distributions across the vascular and interstitial space were calculated as described in the Methods section. Initially,cases of tumours with varying degrees of necrosis as described in Table 1 were evaluated to understand the effect of vascular distribution on fluid flow properties.Figure 2 shows various fluid flow properties calculated for the different tumour models. As shown in Fig. 2A,interstitial pressure at the outer surface of vessels is higher close to the tumour center than in the peripheral region. For the non-necrotic tumour case, the hydrostatic pressure reaches a maximum of approximately 12 mmHg at the surface of some tumour vessels with an average surface pressure of 8.6 ±1.41 mmHg. Pressures at the surface of vessels in the normal tissue region are much lower (3–4 mmHg),which can be attributed to the low hydraulic conductivities prescribed to these vessels. The predicted interstitial pressures in the normal tissue region are within the range for normal tissues (−8 to 6 mmHg) (Kurbel et al. 2001). In the 3% necrotic tumour case, the average interstitial pressure on the surface of the vessels is lower compared to the nonnecrotic case, whilst in the tumour with a higher degree of necrosis (26%) the average interstitial pressure on the vessel surface drops even further. The predicted interstitial pressures in the tumour region are within the range of measured values (4–50 mmHg) reported in the literature (Less et al. 1992). Figure 2B shows the transvascular flux distribution across the tumour network where a net outward fluid flux is observed in intratumoural vessels, which is significantly higher than for vessels in the normal tissue region where the net fluid flux is close to zero. It also shows that transvascular flux is lower in the central tumour region than in the periphery.

Fig. 2 A Pressure distribution on exterior wall of vascular network. B Transvascular flux distribution. C Interstitial pressure distribution in tumour tissue at mid-plane (z = 3000 μm)

Figure 2C shows the IFP distribution at a mid-height(z= 3000 μm) in the tumour tissue. IFP at this location is highest in the central tumour region reaching a peak of 11 mmHg which then drops rapidly moving away from the tumour center where a pressure of around 2.5 mmHg is calculated at the tumour surface. This IFP profile corresponds with experimental and modelling based studies in the literature and further justifies the validity of the model in capturing the flow dynamics in tumour tissue. It is clear that the IFP profile is sensitive to the degree of necrosis, with a pronounced reduction in IFP in the avascular core in the 26% necrotic tumour.This is not so apparent in the tumour with 3% necrosis where high pressures are still seen in the tumour core,but with a slightly lower magnitude compared to the highest pressure found towards the periphery of the tumour. In the 26% necrotic tumour, IFP at the tumour center is significantly lower than the highest pressure distributed at the tumour periphery. The necrotic tumours exhibit pressure profiles that differ from the non-necrotic tumour.

Effect of vessel pruning and blind end removal

The vascular networks examined so far feature a number of blind ends that may affect the pressure drop and flow in the vascular network. To quantify this effect, the tumour vascular geometries were modified by removing several blind ends to open up the pathways for flow. To achieve this, the blind ends were identified within the network, which were then either removed or assigned as venous end vessels depending on their order within the network. Likewise, an AV shunt was also placed in a region near the arterial vessels. The macroscopic properties of the tumour microvasculature are similar as seen in Table 2. To quantify the conductivity of the network for blood flow,parameters including mean vascular segment transit times (VSTT) and percentage of hypo-perfused vessels were calculated and compared for the original and pruned networks. As shown in Table 2, removing blind ends leads to a significant reduction in VSTT and the number of hypo-perfused vessels; the latter is reduced by approximately 34%.

Table 2 Morphological and hemodynamic parameters for tumour vasculature

Figure 3A shows the predicted intravascular pressure in the non-necrotic tumour before and after pruning. It is clear that pressure drop along vessels in the pruned vessel network is increased which can promote flow.Figure 3B shows the interstitial pressure distribution at the mid-plane (z= 3000 μm), which displays a similar profile but much lower IFP in the pruned network.These findings suggest that although the vascular distribution and macroscopic parameters of the vasculature in the tumour tissue are almost similar between the original and pruned networks, microscopic changes such as the removal of blinds and addition of AV shunts can significantly alter the flow dynamics in the vasculature and consequently in tissue space.

Effect of vascular remodelling

To examine the effect of vessel adaption on fluid flow in tumours, the vessel remodelling framework was applied. This was achieved by discretizing the time over which the vascular network grew and calculating the corresponding fluid flow at each time point. Shear stress and time were used as input data to simulate changes in radius, pore size and wall thickness. The evolution of the vascular network with time is shown in Fig. 4.

Figure 5 shows the transvascular flux, vessel surface pressure and interstitial pressure distribution in the adaptive tumour vascular network. The interstitial pressure at the surface of the vessel is seen to besignificantly higher compared to the non-adaptive network (Fig. 3), reaching a maximum of around 28 mmHg, which is still within the interstitial pressure values measured in tumours. The transvascular flux is seen to be heterogenous with most vessels exhibiting a flux in the range of 0–0.2 μm/s. Several less mature vessels at tumour center can be seen with an inward flux of <−0.1 μm/s. This is because more recently formed vessels have larger pores and hence the hydraulic conductivity of these vessels is higher. As a result, the IFP is significantly increased in the adaptive network with a maximum pressure of 26 mmHg. In Fig. 5D, IFP in the high range (≥20 mmHg) is displayed to provide a clearer image of the heterogeneity. When a uniform hydraulic conductivity was assumed as shown in Fig. 3,high pressures are mostly concentrated in the central region. The inclusion of adaptive vessel hydraulic conductivity based on vessel maturity and flow dynamics in the tumour leads to a more heterogeneous IFP distribution even in the case of a well vascularized tumour as seen in Fig. 5D.

Fig. 3 A Intravascular pressure in the original and pruned vascular networks. B Interstitial pressure in tumour tissue at mid-plane(3000 μm) for the two models

Fig. 4 Simulated growth of vascular network with time

DISCUSSION

In this work, a mathematical model was developed where the tumour vascular network is described at a microscopic level on the whole tumour scale and is integrated with Pozrikidis’s fluid flow model to investigate the flow dynamics in tumour tissue. The angiogenesis model was extended to 3D which made it possible to capture the explicit distribution and geometry of the tumour vasculature and its abnormal features. The tumour model was coupled to a fluid flow model which has a distinct advantage of enabling vascular and interstitial flow to be strongly coupled through integration of flow in the vessels and interstitial space over the vessel surface, to provide approximations for fluid flow and transvascular flux that are grounded in physical reality. By applying fluid flow-based boundary conditions on the vessel wall and integrating the pressure over the vessel surface, explicit morphological features of the vasculature including each vessel’s orientation, radius and length can be incorporated. The model was integrated with a vascular remodelling framework that allows for vessel permeability to vary depending on their maturity and local shear stress. The integrated angiogenesis and fluid flow model was evaluated by comparing the obtained morphological and hemodynamic parameters with data in the literature, demonstrating a good agreement.

Using this model, we have been able to demonstrate the strong coupling and interplay between vascular and interstitial flow. Whilst previous studies treating the vasculature as a uniformly distributed source term found IFP to be uniformly elevated across the tumour,in this work we showed that heterogeneous vascular distribution can result in non-uniformly elevated IFP,leading to heterogenous transvascular pressure gradients. Our simulations of fluid flow in tumour vascular geometries with varying degrees of necrosis showed that increasing the necrotic core size resulted in a lower maximum IFP. This finding agrees with a macroscopic study by Soltani and Chen who found that for a small tumour of 1 mm in radius, the presence of a necrotic core as small as 10% of the tumour radius caused a reduction in the maximum IFP (Soltani and Chen 2011). In our case the tumour radius was 1.5 mm.Our results further revealed a non-uniform IFP distribution in necrotic tumours where a dip was seen at the tumour centre which became more apparent with increasing size of the necrotic core. This finding has an important implication for the transport of chemotherapeutics. A dip in pressure at the core would allow for pressure gradients to form from high levels in the tumour periphery to low levels in the necrotic core.For large drugs and liposomes that diffuse slowly, such a pressure gradient would drive convective transport of drugs towards the necrotic region where cells are hypoxic and more resistant or are dead, thus limiting the ability of the drug to target actively proliferating cells. In addition, the large pressure drop at tumour boundary with normal tissue can result in drug clearance from the tumour into the normal tissue,reducing the residence time of the drug in the active tumour regions. This proposition is supported by observation from an experimental study by Vavra et al.who investigated the distribution of molecules delivered intravenously in subcutaneous RG-2 tumours(Vavra et al. 2004). They found that within 2 h following intravenous administration, the molecules accumulated mainly in the tumour’s boundary and necrotic regions with significantly lower concentrations in the viable tumour region. On the other hand, the transport of intravenously administered drug in nonnecrotic tumours might be limited by low extravasation from the blood vessels as evidenced by the lower transvascular flux in the tumour center where IFP is the highest (Fig. 2B). Skliarenko et al. demonstrated that higher IFP correlated with enhanced cell survival following treatment with a vascular disrupting agent(Skliarenko et al. 2006).

Fig. 5 Results for the adaptive tumour network. A Vessel surface interstitial pressure. B Transvascular flux. C Interstitial pressure distribution at mid-plane (z = 3000 μm) showing full pressure range. D Interstitial pressure in the high range (≥20 mmHg) to show heterogeneity

In order to investigate the effect of microscopic features of the vasculature, the flow was analysed in two vascular networks with an almost identical vascular distribution and macroscopic parameters that differed only in the removal of blind ends and addition of AV shunts. These subtle changes in the vascular network opened flow pathways and resulted in a pronounced increase in intravascular pressure drop across the tumour network and a reduction in IFP. This highlights that although vascular networks may appear almost identical from a macroscopic point of view as determined by their vascular density and microvascular distribution, small differences on a microscopic scale can result in dramatically altered flow properties. The ability of normalizing vessel networks has been suggested as a target to enhance treatment (Cantelmo et al. 2017; Martin et al. 2019). Chauhan et al. simulated vessel normalization by reducing the permeability and heterogeneity in vessel permeability, showing that IFP was reduced significantly (Chauhan et al. 2012). IFP reduction was also reported by Sweeney et al. who investigated vessel normalization by reducing the vessel diameters and regulating vascular hydraulic conductivity values (Sweeney et al. 2019). Experimental studies have shown that the treatment of tumours with anti-angiogenic drugs to normalize vasculature resulted in a significant reduction in IFP (Tong et al. 2004). Blind ends and self-loops are abnormal features of the tumour vasculature; they cause flow stagnation and reduce perfusion to nearby regions of the tissue. In our work,vascular normalization was simulated by reducing the abnormalities in the architecture of the tumour vasculature by removing blind ends and self-loops. The resultant increase in intravascular pressure drop can promote perfusion and allow intravenously injected drugs to distribute more effectively within the vascular network. In addition, the reduced IFP would allow the drug to permeate through the vessel wall and limit efflux of the drug to the surrounding normal tissue thereby increasing drug residence time.

Finally, the influence of vascular remodelling was examined, which was incorporated by allowing the radius, porosity and wall thickness of each vessel to vary depending on the age and flow properties. Whilst most of the previous studies investigated fluid flow in tumours using a uniform vessel hydraulic conductivity,a vascular remodelling framework was adopted in this study that explicitly defines hydraulic conductivity values for each vessel based on its maturity and perfusion. Our results showed that accounting for vascular remodelling caused a dramatic increase in IFP.The non-uniform vessel permeability led to heterogenous transvascular flux across most of the tumour network. Chauhan et al. incorporated a heterogeneous pore size distribution in their model(Chauhan et al. 2012), they assigned each vessel a pore size assuming a unimodal distribution throughout the vasculature and hence did not explicitly incorporate the effect of fluid flow or vessel maturity on the permeability of the vessel. Their work showed similar findings where IFP was heterogenous, but it could not be determined from their presented data whether this was due to the heterogeneity of the vessel pore size or the vascular distribution. In the model developed by Vavourakis et al., IFP was calculated incorporating heterogenous vessel permeability, however, their model applied the homogenization technique whereby vascular properties were averaged within each element and IFP was evaluated as a function of averaged vascular property parameters within the element (Vavourakis et al. 2017). Hence, the morphological effect of the vascular network was not explicitly incorporated when determining IFP. In their work, the tumour IFP values were found to be elevated, however the radial pressure profile was relatively uniform across the tumour before dropping rapidly at the boundary. In our study, by applying a fluid flow model that takes into account the explicit morphology of the vasculature, we found IFP to be more heterogeneously distributed with steeper pressure gradients. The predicted IFP profile in the tissue could also pose an obstacle to the transport of therapeutics with steeper gradients at the tumour periphery and increased interstitial fluid velocity potentially promoting efflux of therapeutics from the tumour tissue. Non-uniform gradients were found across the tumour tissue with IFP peaks being distributed heterogeneously within the tumour. As described previously, high IFP can influence the proliferation of cancer cells and their chemosensitivity,therefore, heterogenous IFP values with steeper gradients could lead to a non-uniform response to chemotherapy in the interstitial space. The non-uniform transvascular flux can result in intravenously administered drugs to extravasate from the vessels at greater concentrations in specific tissue regions and allow other regions to escape treatment.

The angiogenesis model used to describe the tumour vasculature was able to capture the heterogeneous properties of the vasculature and incorporate various aspects of the angiogenesis process. In the model, the tumour was assumed to be of a fixed size and the effect of tumour growth during angiogenesis was neglected.This can be justified as the time scale on which tumour growth occurs is much larger than that of vascular growth. However, changes in oxygen concentration caused by vascular remodelling were not incorporated which can influence the distribution of chemical species secreted which would in turn influence the angiogenic process and vascular network generated. Mechanical forces such as solid stress induced by the growing tumour can affect the movement of endothelial cells which has not been incorporated in the angiogenesis model (Welter and Rieger 2010). The stress caused by proliferating cancer cells can compress and shut down the function of some vessels. These factors could influence the spatial arrangement and orientation of the vascular networks produced (Jain et al. 2014). These limitations were not addressed in the current study so as to keep the model computationally tractable. In the fluid flow model, blood was modelled as a continuum Newtonian fluid with a constant viscosity. Blood is in fact a non-Newtonian fluid composed of multiple cellular components, which can have an important influence in small capillaries. Spatial and time variations in these components may result in variable viscosity within the vessel network. Additionally, the response of vessel diameter to changes in transmural pressure and metabolic stimuli has been neglected. Existing models that incorporate these changes were developed using data obtained from vessel networks in normal tissues(Pries et al. 2001, 2009). Tumours exhibit a different functional behaviour and would be expected to respond differently. When considering interstitial flow, the tissue was assumed to be a homogenous porous medium. In reality, tumour tissue is heterogeneous with varying porosity and interstitial hydraulic conductivity which consequently affect fluid flow. Further research can be built on this model to incorporate tissue heterogeneity and more complex angiogenesis models.

CONCLUSIONS

Previous models that assume a uniform vascular distribution or treat the vasculature as a uniform source term do not capture the effect of neighbouring vessels on each other, which could potentially influence flow on the macroscopic scale. By using the boundary integral method to distribute sources and dipoles along the vessel surface, our model allows for the topology and geometry of the vasculature to be taken into account, so that the effect of each vessel segment on IFP at any point in space is considered. The model also accounts for the influence of IFP on intravascular flow,in an attempt to provide a more realistic prediction of fluid flow under tumour physiological conditions. Using this method, the work presented in this paper has demonstrated the importance of considering not only the explicit nature of the tumour vasculature but also the intravascular flow properties when predicting and evaluating fluid flow in tumour tissue. Our results shed further light on the intricate flow dynamics occurring on a microscopic scale in the tumour tissue and their potential consequences on macroscopic flow and drug transport. Our model can be further extended to predict drug distribution in tumour tissue, and this will be reported separately.

METHODS

Tumour model

Tumour induced angiogenesis

The tumour model is generated to provide an explicit representation of the tumour vasculature at a microscopic scale. The extravascular space is assumed to consist of uniformly distributed cancer cells and extracellular material. To generate the explicit tumour vasculature, Anderson and Chaplains mathematical angiogenesis model is implemented (Anderson and Chaplain 1998). The model describes the movement of endothelial cells through random motility, chemotaxis in response to tumour angiogenic factor gradients(TAF) and haptotaxis in response to fibronectin. The original 2D angiogenesis model has been extended to 3D to describe the movement of endothelial cells along each orthogonal axis in a 3D space. The nondimensional equations, endothelial cell density (n),fibronectin concentration (f) and tumour angiogenic factor concentration (c), which describing vascular growth are given below:

where subscripts l , m and w denote position on the grid anddenotes dimensionless timestep.

Initial conditions are defined for a spherical tumour centred in the domain with a non-dimensionless radius of 0.25 surrounded by normal tissue. Hence, TAF concentration is expected to be highest in the tumour region as it is secreted by tumour cells with concentration decreasing as distance from the tumour increases. Initial TAF concentrationis defined as follows:

Fibronectin is assumed to be highest in regions in and around parent vessels. Anderson and Chaplain justified this by referring to experimental observations on increased extracellular molecules around parent vessels(Anderson and Chaplain 1998). This relationship is adapted in our model and is described by the following:

A tumour with a necrotic core of radius 0.08 can be generated by simulating a negative TAF gradient at the center with a defined size as follows:

The endothelial cells are unlikely to move down the concentration gradient of TAF. Eq. 5 and 6 are used to initialise the model and a spatiotemporal evolution of the endothelial cell distribution is obtained using Eq. 4.Figure 6A shows the initial TAF distribution for the nonnecrotic and necrotic core tumours. Figure 6B shows the spatio-temporal distribution of endothelial cells where the concentration is zero in the necrotic region defined.

The endothelial cell distribution field is combined with a discrete model to track the movement of tip endothelial cells and resulting vascular network formed.The migration of the tip endothelial cells is determined by seven probability coefficients P0–P6which correspond to the probability of no movement (P0),movement along the x-axis (P1,P2), movement along the y-axis (P3,P4), and along the z-axis (P5,P6). The magnitude of the probability coefficients depends on local TAF and fibronectin concentrations. At each time step the seven probability coefficients are calculated to produce seven ranges between each coefficient. A random number is generated between 0 and 1 and the movement of the endothelial cell depends on the range where the random number falls. Hence, the movement of the tip endothelial cell is based on three components:random, chemotactic and haptotactic movement.Branching and anastomosis are common features in tumour vasculature, and these are represented following the approach of Anderson and Chaplain(Anderson and Chaplain 1998) which is summarized here. Branching from a sprout is dependent on three factors: (1) branching age — a sprout must reach a mature state, measured by a threshold branching age,for it to branch; (2) space — there must be space for the sprout to branch into, checked by ensuring there are no other sprouts around the tip; (3) the number of endothelial cells at the tip — an endothelial cell density threshold must be satisfied for the sprout to branch. If the above three conditions are met, the branching of a sprout is given as a probability that is dependent on the local TAF concentration. The higher the TAF concentration, the greater the chances are for a sprout to branch. Overall, this model describes the movement of a sprout tip from the parent vessels where there is little branching due to the age and as the sprouts reach the tumour, the increased age combined with the high TAF concentration increases the chance of vessel sprouting.Anastomosis or merging of vessels is described simply when a sprout moves into a space occupied by another sprout. The sprout that continues to exist is chosen at random. Combining the discrete method with the fields generated for endothelial cell density, TAF and fibronectin, vascular network geometries can be generated as seen in Fig. 7. Figure 7B highlights the avascular core region in the necrotic tumour that results from the presence of a negative TAF gradient at the tumour core.

Fig. 6 A Initial endothelial cell distribution. B Endothelial cell distribution after 40 days

Fig. 7 A Capillary network formed after 40 days for a non-necrotic and necrotic tumour. B Core region of non-necrotic and necrotic tumour

Vascular network regulation

In order to apply the fluid flow model to the vascular network generated from the angiogenesis model, the geometry of the network must be regulated and further defined. First, nodes are assigned to coordinates where tip merging or branching occurs and where the tips initially sprout and end. Then further nodes are created along the path that the tip moves. A connection matrix logging the connectivity between the nodes is formed based on the data from the raw vascular network generated. A depth-first search algorithm is then used to move down from the initial node, creating segments based on the connection matrix. To define the diameter of the vessels, a maximum diameter is set at the initial parent vessels which are assumed to be maintained for subsequent vessels down the network, except in the case where branching or merging occurs. When this occurs the diameter of the daughter vessel is assumed to decrease in size. The change in radius from a vessel of the nthgeneration to the (n+1)thgeneration can be modelled as monotonously decreasing as follows:

where Dminand Dmaxare defined minimum and maximum diameters, and m is a coefficient less than 1.Applying the vascular regulation algorithm to the networks generated from the angiogenesis model yields the vascular geometries shown in Fig. 8.

Fig. 8 Generated capillary network for non-necrotic and necrotic tumour

Vascular remodelling framework

Blood vessels in tumour tissues are characterized by their heterogeneous nature where they can exhibit varying diameters and permeabilities. Blood vessels naturally adapt in response to mechanical and metabolic stimuli in order to ensure the tissue is well perfused and supplied with nutrients and oxygen. The mechanisms of vessel adaption have been described by Pries et al., which provides a framework to model changes in vascular structure in response to stimuli(Pries et al. 2001). The parameters in their model were fitted specifically to the normal tissue which was studied in their work. Cellular abnormalities in the tumour vessels such as large inter-endothelial gaps limit signal conduction, whilst the lack of pericyte coverage reduces the ability to undergo diameter changes. The impaired structural adaption in addition to the lack of experimental data on tumour vessels makes the Pries model less applicable in this case.Hence, a simple framework is employed in this study,which was developed by Vavourakis et al. to model changes in vessel radius, wall thickness and pore radius(Vavourakis et al. 2017). In their work, vessel adaptivity is described by a single variable, the remodelling timewhich is governed by the shear stress on the vessel and is given by Eq. 9:

wheretm−Tis the time value for a vessel to reach the upper limit of remodelling if τ ≥ τref, andtm−0is the time required for vessel remodelling if shear stress is zero. Here, we setto 100 d andtm−Tto 10 d based on the empirical value adopted by Vavourakiset al.(Vavourakiset al. 2017). The model assumes that a vessel becomes fully remodelled as it moves from a poorly perfused state to a well perfused state which is defined by the reference shear stress, τref. Vascular segment transit time (VSTT) is used as an indicator for hypo perfusion as it takes into account blood flow and vessel length in each segment. Kamoun et al. defined hypo-perfused vessels as those where velocities are below 0.05 mm/s, hence for an average vessel of 0.2 mm in length, the VSTT is approximately 4 s(Kamoun et al. 2010). The reference shear stress can be evaluated in terms of VSTT as follows:

where Rmaxis the maximum radius. Wall thickness,,can be assumed to increase linearly as a function of time and is given as

where wmaxand wminare the maximum and minimum wall thickness. The vessel pore size is also modelled as a function of time and shear stress and is expressed as

where rpminand rpmaxare the minimum and maximum pore size respectively. Initially when a vessel is generated, it is assigned a radius Rmin, pore sizerpmaxand wall thickness wmin. This follows physiological behaviour as when a vessel initially sprouts and is in an immature state, it is highly porous with thin walls. The hydraulic conductivity of each vessel is determined as a function of wall thickness, pore radiusand the fraction of vessel wall occupied by pores

Fluid flow model

Fluid flow in the tumour model can be described in three parts, vascular flow, transvascular flux and interstitial flow. Blood flow is modelled using Poiseuille’s law which describes flow rate,, through a straight cylindrical segment as a function of radiusand pressure drop. Transvascular flux,, is described by starling’s filtration law whilst interstitial flow is described by Darcy’s law through porous medium.The governing equations used to model these flows are summarised as follows:

Applying Pozrikidis’s method (Pozrikidis and Farrow 2003), flow is assumed to be conserved along the vessel length to couple vascular and interstitial flow,

Laplace's Eq. 18 for interstitial pressure is reformulated using the boundary integral formulation as in Pozrikidis’s work to give the interstitial pressure at field pointin terms of a combination of a single- and a double-layer potential defined over the surface of the vasculature. The pointis set at the surface of vasculature and the interstitial pressure is assumed to be independent of angular pressure giving,

Substituting Eq. 15 and 16 into Eq. 19 and assuming a constant radius for each vessel segment gives

Numerical methods

The vasculature is discretized intocylindrical segments where j = 1, 2, 3, ··,and each vessel assigned a radius, lengthand vascular hydraulic conductivity

For a single vessel with uniform radius, Eq. 21 is discretized using the second order finite difference method as in Pozrikidis’s work to obtain the midpoint vascular pressure in terms of the interstitial pressure and end-point vascular pressures,

Eq. 26, 27 and 28 provide a linear system of equations which is solved subject to arterial and venous boundary pressures and interstitial tumour surface pressure.Initially, inlet and outlet capillary nodes are assigned arterial and venous pressure parterialand pvenousrespectively, while interstitial pressures at the surface of the capillary midpoints are set to. Eq. 26, 27 and 28 are used to generate a tridiagonal system of equations which is solved in MATLAB to calculate nodal capillary pressures and interstitial pressure at the surface of the vessel segment.

Once the nodal capillary pressures and the interstitial pressure at the surface of the segments are determined,the interstitial pressure in the tissue space is solved using the discretized boundary integral formulation:

The interstitial space is discretized and the singleand double-layer potentialsandare used to evaluate the influence of vessel segments on each discretized tissue point. Using the influence coefficients, vessel topology, nodal capillary pressure and interstitial pressure at the vessel surface, the contribution of the vessel towards the pressure in the interstitial space at pointis calculated, where= (x,y,z). The interstitial pressure at pointcan then be calculated as a sum of these contributions and the influence of the tumour surface pressure. Rearranging Eq. 22 gives a formulation for the transmural pressure termand is substituted into Eq. 29 to give the interstitial pressure at pointas follows:

Computational method

The angiogenesis and fluid models were implemented in MATLAB and computations were performed using a 64-bit Intel (R) Core (TM) processor (3.40 GHz) with 64 GB RAM. The computational domain was 6 × 6 ×6 mm3, which was discretized into a structured grid with uniform spacing of 150 μm along each orthogonal axis. Using this hardware and discretization parameters, each simulation of vascular network took approximately 90 s. The system of linear equation generated in the fluid flow model was solved in MATLAB to resolve vascular and vessel surface IFP. The tissue space was discretized into a structured grid with a uniform spacing of 85 μm where at each point the IFP was determined as described in Eq. 30.

Model parameters

The parameters used to model the tumour geometry are given in Table 3. In the experimental work that Anderson and Chaplain based their angiogenesis model on, the average distance from the tumour to the parent vessel was 1–2 mm, hence a length scale ofL= 1.5 mm is assumed in this study.

Table 3 Tumour geometry parameters

Vessel hydraulic conductivity is based on the ability of the vessel to allow fluid movement across the vessel wall. Capillary filtration measurements in rat implanted tumours showed hydraulic conductivity values ten times higher than those reported in normal tissue.Vessels in the normal and tumour tissue region were prescribed vessel hydraulic conductivities. A summary of key parameters used in fluid flow simulations are given in Table 4.

Table 4 Fluid flow parameters

Morphological and hemodynamic analysis

In order to validate the tumour network generation model and the fluid flow model, a list of morphological and hemodynamic parameters were calculated and then compared with literature values obtained from empirical studies. For morphological analysis, the properties, tumor volume (Vt), vascular density (Vd),length density (LD), vessel surface to volume ratio(S/V), and maximum extravascular diffusion distance(R), were determined (Kimet al. 2012; Stamateloset al.2014):

Compliance with Ethical Standards

Conflict of interest Moath Alamer and Xiao Yun Xu declare that they have no conflict of interest.

Human and animal rights and informed consent This article does not contain any studies with human or animal subjects performed by any of the authors.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use,sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.