Enhanced perfusion following exposure to radiotherapy: A theoretical investigation

Tumour angiogenesis leads to the formation of blood vessels that are structurally and spatially heterogeneous. Poor blood perfusion, in conjunction with increased hypoxia and oxygen heterogeneity, impairs a tumour’s response to radiotherapy. The optimal strategy for enhancing tumour perfusion remains unclear, preventing its regular deployment in combination therapies. In this work, we first identify vascular architectural features that correlate with enhanced perfusion following radiotherapy, using in vivo imaging data from vascular tumours. Then, we present a novel computational model to determine the relationship between these architectural features and blood perfusion in silico. If perfusion is defined to be the proportion of vessels that support blood flow, we find that vascular networks with small mean diameters and large numbers of angiogenic sprouts show the largest increases in perfusion post-irradiation for both biological and synthetic tumours. We also identify cases where perfusion increases due to the pruning of hypoperfused vessels, rather than blood being rerouted. These results indicate the importance of considering network composition when determining the optimal irradiation strategy. In the future, we aim to use our findings to identify tumours that are good candidates for perfusion enhancement and to improve the efficacy of combination therapies.


Introduction
Tumour cells exploit various biological opportunities to disregard the extracellular signals that normally direct an organism's somatic cells to sacrifice themselves for the sake of its germ cells [1].One of these opportunities is a tumour's ability to incite the growth of its own blood supply network in a phenomenon called 'tumour angiogenesis'.Tumour angiogenesis produces vasculature that grows in a disordered manner.Therefore, blood circulation in a tumour is sluggish, which leads to poor perfusion and a lack of oxygen.Poor oxygenation, in turn, attenuates the effect of major oncological therapies.The efficacy of radiotherapy, which is used to treat more than half of all cancer patients, depends on the generation of reactive oxygen species that induce irreparable DNA damage [2].The mechanism governing the increased efficacy of radiotherapy in the presence of oxygen is explained by the 'oxygen fixation hypothesis' [3,4].Upon irradiation, free electrons are generated that mediate irradiation-induced damage, either by directly damaging the macromolecules, resulting in sites of DNA base loss and DNA singlestrand breaks (SSBs), or by interacting with oxygen to form peroxyl radicals, which induce DNA damage, including double-strand breaks (DSBs) that are difficult for cells to repair [5,6].In hypoxia, the lower amount of oxygen results in reduced formation of reactive oxygen species (ROS), including peroxyl radicals.Hence, there is less DNA damage, leading to a greater ability for the damage to be repaired, resulting in increased radioresistance.Importantly, pockets of hypoxia within a tumour can exhibit a lower response to radiotherapy by up to a factor of three [7].Moreover, adaptations to hypoxia can lead to tumour phenotypes with increased chemoresistance and metastatic capabilities [8].For a list of drugs that perform more poorly in hypoxic conditions than in normoxic conditions, see [9].
Tumour growth is primarily monotonic, i.e., the tumours continue to grow and the tumour vasculature also grows continuously after the angiogenic switch has been triggered.Consequently, there is no spontaneous regression of the tumour blood vessels, resulting in a chaotic and disorganised vasculature with limited remodelling and maturation.In the development of normal vessels, one of the important mechanisms regulating their complexity and organisation is the trimming of a vascular network, termed 'vessel pruning', which marks the physiological regression of a subset of microvessels within a growing vasculature [10][11][12].This phenomenon was studied in detail in the postnatal mouse retinal vascular network where it was shown that the vessel density at the sprouting front of a postnatal retina is significantly higher than in the mature vascular network of the adult retina.The main driver of vessel pruning is the vascular endothelial growth factor (VEGF) and VEGF receptor (VEGFR) signalling.However, other signalling pathways activated by fibroblast growth receptor 2 (FGF2), angiopoietin 2 (ANG2), platelet-derived growth factor (PDGF), and delta-like 4 (DLL4) have been shown to work independently of the VEGF/VEGFR pathway [11].Anti-angiogenic therapies take advantage of these signalling axes and can cause a transient induction of vessel pruning in tumours, leading to a more normal-looking and stabilised vasculature.This is called 'vascular normalisation', which involves not only morphological but also functional improvements: reduced interstitial fluid pressure (IFP), reduced tumour hypoxia, and improved penetration of macromolecules from these vessels, which significantly improves the efficacy of some therapies, such as radiotherapy and chemotherapy [10,13].Recently it was shown that radiotherapy can also act as a vascular normalization therapy by preferentially pruning the small, non-perfused vessels [14,15].Importantly, non-perfused vessels and vessels with sluggish blood flow are also predisposed to vessel pruning and changes in blood flow act as regulators of vessel regression or maturation.In the mouse retina, a pharmacological inhibition of vasoconstriction or induction of vessel dilation resulted in reduced blood vessel stability confirming that vessel pruning is dependent on haemodynamic forces [16,17].However, it is not yet clear whether molecular components, such as angiogenic factors, or mechanical factors, such as shear stress, are the primary regulators of vessel regression associated with flow.
In the clinical setting, irradiation of the normal tissue surrounding the tumour cannot usually be avoided.In order to reduce the normal tissue toxicity of radiotherapy, the total radiation dose is divided into several fractions.'Conventional fractionation' has been used in clinics for decades and divides the total radiation dose into 1.8-2 Gy fractions administered 5 days a week up to a total dose of 40-70 Gy.This type of fractionation takes advantage of the biological differences between the tumour and normal tissue, which usually results in less damage to the normal tissue with the same degree of tumour control [5].In contrast to conventional lowdose fractionated therapies, hypofractionation describes a radiation treatment in which the total dose is administered in a smaller number of larger (>2 Gy) fractions, resulting in a shorter overall treatment time.If the treatment is administered as a single dose or in a small number of fractions, usually with a dose of 8-30 Gy per fraction, it is referred to as stereotactic body radiation therapy (SBRT) [18].The use of SBRT has increased over the past decade.This is primarily due to technological advances in image guidance and highly conformal dose delivery, which allow a high dose of radiation to be delivered to the tumour while maintaining an acceptable dose to surrounding normal tissue [5,18].
Experiments conducted to measure the change in tumour perfusion after irradiation of tumours have recorded varying outcomes.A review by Park et al. found that reported results were inconsistent, and that functionality improved and then worsened in human tumour vasculature during the early stages of radiotherapy [19].Conversely, in human tumour xenografts or murine tumours, irradiation resulted in mild to severe damage depending on the dosage, reducing blood perfusion [19].Kim et al. found that stereotactic body radiation therapy decreased perfusion [20].On the other hand, Shibuya et al. found that in the case of cervical cancer, blood flow improved significantly after irradiation [21].Moreover, a study by Bussink et al. found that irradiation led to rapid changes in perfusion, observing increases shortly after irradiation followed by significant decreases [22].Kaeppler et al. used a transgenic CreERt2-td-Tomato mouse model with <95% of endothelial cells fluorescently labelled coupled with in vivo intravital multiphoton microscopy which allowed them to simultaneously image the perfused and non-perfused tumor vasculature to study its response to single-dose and fractionated radiotherapy [15].They used two tumour cell lines-a highly vascularised colon adenocarcinoma (MC38) and a less vascularised melanoma (B16F10)-and found that the former contained a larger proportion of hypoperfused vessels compared to the latter.Moreover, it was observed that the smaller-diameter vessels were more likely to be hypoperfused and also more likely to be pruned following irradiation (i.e., their endothelial cells were more susceptible to apoptosis).Thus, it was concluded that tumour perfusion post-irradiation depends on the density of small (typically hypoperfused) vessels.
An improved understanding of the relationship between irradiation and perfusion would allow for accurately-planned fractionated radiotherapy.However, it remains unclear why certain tumours exhibit increases versus decreases in perfusion after irradiation.The blood vessels that comprise tumour vasculature have heterogeneous architectural characteristics that present a complex system resistant to experimental analysis.Thus, there is a need for computational tools to guide and complement experimental research.
Broadly speaking, tumour vascular architecture has been modelled spatially with synthetic networks that reflect biological characteristics or, more recently, by digitising the geometry of real tumour vasculature.Regular forking networks, constructed by Bernabeu et al. using pathological dimensions gleaned from experimental data, are highly ordered because they are generated using a simple set of deterministic rules and the resulting geometry (vessel diameters, lengths, and arrangement of vessels into the network) is thus very regular [23].Moreover, these networks exhibit strict hierarchy in that every parent vessel subdivides into two daughter vessels that are shorter and thinner than the parent.Non-hierarchical and ordered vasculature has been previously represented as a network with a repeated hexagonal unit, similar to that observed in avian yolk sacs [24].Alarco ´n et al. have used a similar network to couple processes at the intracellular, cellular, and tissue scales [25].Owen et al. employed both hexagonal and disordered networks to show that vascular remodelling could be achieved with a balance of pruning and angiogenesis [26].
Synthetic networks are also commonly generated by simulating the random migration of capillary tips based on a chemotactic gradient.Anderson and Chaplain developed a model using a random walk biased towards higher TAF (tumour angiogenesis factor) levels, which resembled in vivo angiogenic networks [27].Similar models of biased random motility and sprout formation have been employed by Macklin et al. and Shirinifard et al. in 2D and 3D, respectively [28,29].In addition, Stepanova et al. have employed a multiscale random walk model to study endothelial cell dynamics [30].The modelling of angiogenesis still remains an active area of research [31][32][33].
While synthetic networks can replicate many features of real tumour vasculature, some studies have gone a step further and digitised experimentally-acquired tumour vasculature.Grimes et al. used a 3D digitised network to estimate oxygen distribution [34].Grogan et al. used a similar method to compare 2D and 3D representations, while Sweeney et al. used a 3D model to suggest that using realistic vasculature was key to modelling tumour fluid dynamics [35,36].However, such simulations presuppose the knowledge of boundary conditions and the methods for direct observation of flow in individual vessels are limited (especially smaller microvessels) [37].
Topological data analysis (TDA) is an emerging mathematical field that uses topological and geometric approaches to quantify the 'shape' of data [38,39].TDA characterises shape via topological invariants such as connected components and loops at multiple scales.The most prominent method from TDA, persistent homology, has been successfully applied to quantify dynamic characteristics of vascular networks in experimental data [40,41] and to distinguish between synthetic vascular networks produced by different parameter regimes in a mathematical model of tumour-induced angiogenesis [42].Persistent homology tracks topological features across a filtration, e.g., a nested sequence of vascular networks, and outputs information on the persistence of these features in the filtration in the form of a barcode [39].The information contained in the barcode can either be vectorised and integrated with statistical or machine learning tools [43,44] or it can be compressed into a single topological descriptor of the data, e.g., the number of loops across the filtration, as in [41].In experimental data, the change in the total number of loops in vascular networks over time following radiotherapy has shown great variation across different tumours relative to the day of irradiation (single dose) [41].However, this study did not take perfusion into account and, in particular, did not differentiate between different initial structural characteristics of the networks pre-therapy which can affect their response to radiotherapy.More recently, a novel topological descriptor for weighted directed graphs was developed [45].For a weighted directed graph, e.g., a vascular network weighted by vessel diameter or flow time with the directions given by the flow, the method outputs a barcode describing its structure.While this topological descriptor encodes rich information about such a network, it has not yet been applied to real data, is computationally intensive, and would need to be substantially compressed to enable interpretable comparisons of multiple networks.Here, we focus on combining simple topological and geometrical information to describe the structure of vascular networks interpreted as undirected and unweighted graphs, i.e., disregarding information on flow direction and vessel diameters.We then link these structural features to changes in perfusion via flow simulations where flow direction and vessel diameters are taken into account.We are motivated by previous TDA studies of tumour vasculature [40][41][42], in which the number of loops appeared to reflect significant structural changes of vasculature.In contrast to these studies, we do not apply persistent homology, but focus on the total number of loops in the networks as a topological descriptor at consecutive time points.This topological descriptor is readily computable via the Euler characteristic (see the section Topological determinants of perfusion), does not require nestedness of the networks (in contrast to persistent homology), and is thus more suited to the analysis of networks undergoing pruning.None of the previous studies directly explored the link between network topology and perfusion as altered following irradiation and to the best of our knowledge, there exists no experimentally-motivated computational study of the effect of irradiation-induced pruning on network perfusion that considers both geometric and topological architectural features.
In this paper, we aim to characterise the extent to which radiotherapy-induced changes in vascular geometry (e.g., pruning of small vessels) may contribute to increased vascular perfusion.We identify key geometric and topological descriptors characterizing vasculatures with enhanced perfusion following vessel pruning as well as mechanisms by which such enhancement might occur.While [15] provided a broad analysis comparing two tumour cell lines under various radiotherapeutic treatment regimes, we noticed that even within one cell line (MC38) under single-dose radiotherapy, some tumours increased and some decreased their perfusion.We will therefore focus solely on MC38 under single-dose radiotherapy and aim to assess whether certain structural features pre-irradiation can be used as a proxy for predicting perfusion response to radiotherapy.
The structure of this work is summarised in Fig 1 .Having introduced the biological and computational background for this study, we next describe the in vivo experimental study motivating and informing our model in the section Experimental motivation.In this section, we also highlight key correlations between certain architectural (geometric and topological) metrics and changes in perfusion after irradiation.In the section Model overview, we introduce key aspects of our computational model including the proposed architecture of the initial networks, the pruning rules, and the way in which perfusion is measured.We then use our model to investigate the underlying causal links and present two mechanisms of perfusion improvement in the Results section.We conclude with a contextualisation of our findings and a discussion of their implications in the Discussion section.The complete experimental procedure and model description can be found in the Materials and Methods section.Bottom-left panels present representative vascular regions containing perfused (yellow) and hypoperfused (purple) vessels from Day 0 (just before irradiation) for one of the tumours for which irradiationinduced vessel pruning led to a decrease in perfusion (tumour 6; top) and one for which it led to an increase (tumour 1; bottom).Note that the latter vasculature contained many hypoperfused blunt-ended vessels, while the former contained more loops.Bottom-right panels show pruned synthetic (forking) networks exhibiting similar properties with perfused (red) and hypoperfused (blue) vessels.https://doi.org/10.1371/journal.pcbi.1011252.g001

Uncertainty in direct flow simulations
Ideally, one would assess the impact of radiotherapy-induced vessel pruning on network perfusion by performing direct blood flow simulations on networks extracted from microscopy images.However, information about which network nodes serve as inlets and outlets is often limited or non-existent.Even given such information, it would still be challenging to determine whether constant pressure or constant flow rates should be imposed at inlets, and what values these pressures or flow rates should take to obtain quantitative agreement with measurements of blood flow in vivo.Despite receiving more attention from the scientific community in recent years, the problem of reliably assigning appropriate conditions at boundary nodes remains unresolved.In [37], an algorithm was developed whereby the boundary conditions are fitted so that the deviation of pressures and wall shear stresses in microvessels from appropriately selected target values are minimized.This method has been used in recent years to find appropriate boundary conditions for cerebral haemodynamics as well as in tumour perfusion simulations [46,47].However, due to the complexities of the dynamic tumour microenvironment, a reliable estimation of flow boundary conditions remains an active area of research.In this work, we choose a different strategy.First, in the section Understanding perfusion response to radiotherapy and its determinants, we identify key global geometric and topological metrics that can be used to determine whether irradiation-induced vessel pruning improves or impairs the perfusion of vascular networks extracted from tumours.Then, in the following sections, we present a theoretical study of network perfusion and vessel pruning in synthetic networks.In doing so, we aim to confirm the utility of the newly-proposed metrics in distinguishing tumour vascular networks whose perfusion increases or decreases following radiotherapy.

Understanding perfusion response to radiotherapy and its determinants
Exposure to radiotherapy causes DNA damage in endothelial cells, which was found to induce cell death due to apoptosis preferentially in smaller hypoperfused vessels and result in cell cycle arrest in larger vessels which thus remain functional channels for the flow of blood [15,48].Regardless of the exact mechanism driving the cell death, we assume that radiationinduced DNA damage is the dominant cause of vessel pruning with a strong preference for small vessels and we neglect a preferential regression of hypoperfused vessels found in developmental vascular networks [17,49].
We define perfused and hypoperfused vessels as in [15] (see also the section Experimental procedures and data preprocessing for the experimental procedure and the section Perfusion threshold for the minimum differentiating flow rate in our simulations) and similarly use the perfusion fraction (PF) to quantify network perfusion.The perfusion fraction PðtÞ is the ratio of the number of perfused vessels to the total number of blood vessels, i.e.:

PðtÞ ¼ # of perfused vessels at time t # of remaining ðunprunedÞ vessels at time t : ð1Þ
The perfusion fraction is time-dependent because both the total number of vessels and the number of perfused vessels change over time, due to pruning (see the representative example in Fig 2).Experimental observations of irradiated tumours suggest that vessel pruning occurs on timescales that range from hours to days [15,48,51].We thus focus on short-term responses, during the first four days following radiotherapy (Fig 3).We number the studied tumours from 1 to 7. Firstly, we note that the time between irradiation (Day 0 in Fig 3 ) and the first day on which a reduction in vessel number (vessel pruning) was observed varied across the tumours that were studied.While most tumours (1, 2, 3, and 4) exhibited a decrease in vessel number on Day 1, others (5 and 6) experienced a slight increase in vessel number over this time period.We note further that this delayed decrease in vessel number was observed for tumours with a relatively low vessel number and (at the same time) a large average vessel diameter at the time of radiotherapy (see Table 1).We speculate that these tumours were still undergoing extensive angiogenesis and had not developed a fully functioning tumour microvasculature at the time of radiotherapy.We note also that the vessel count for tumour 7 increased up to, and including, Day 4. This tumour's vasculature had the largest number of connected components per size, and the size of its largest connected component was the smallest across all tumours (see Table A in S1 Appendix).These data are indicative of problems in image processing for this tumour and, therefore, we exclude it from our analysis.Secondly, the number of days over which the number of vessels decreased varied between tumours: for some tumours (1 and 3) the decrease only lasted one day, while for others (2, 4, 5, and 6) it lasted two days.This variability is likely due to the complex nature of the tumour microenvironment and uncertainties in the timescales for vessel pruning.The pruning phase was followed by a period characterised by a significant increase in vessel number, likely due to angiogenesis.Based on the above observations, we divide the tumours into two groups A and B so that a tumour belongs to group A or B if its perfusion fraction increases or decreases respectively during the pruning phase, which we define to be the period from Day 0 to the last day on which the number of vessels decreased, prior to the onset of angiogenesis.For a more quantitative comparison, we also define the pruning-induced perfusion difference DP and its relative counterpart D % P as follows: Table 1.Key geometric and topological characteristics of tumour vasculatures on Day 0. Note that all lengths are in μm, all geometric resistances (total, mean, and per loop) in μm −3 , and all resistances with viscosity in cP�μm −3 , where cP (centipoise) is a unit of dynamic viscosity equal to mPa�s).The perfusion increases if D % P > 0 and decreases if D % P < 0. In Table 2, we summarize these quantities and rank the tumours with respect to the relative change in PF (as measured by D % P).

Change in PF DP
It might seem that the response to radiotherapy can be predicted by whether or not the tumour is initially highly-perfused.Tumours that are initially highly-perfused (e.g., 4 and 6) tend to have impaired perfusion post radiotherapy while those that are initially hypoperfused tend to have improved perfusion post radiotherapy (e.g., 1 and 2).However, even in our small dataset, we can find a counterexample: tumour 5 experienced a low perfusion fraction on Day 0, and an even lower value post-pruning.The particular value of the perfusion fraction (on any day of the experiment) is expected to be strongly dependent on the number, location, and strength of the inlet vessels which are often unknown.However, in the absence of such information and provided the inlet vessels are not pruned, it may still be possible to predict the perfusion response of tumour vasculature to radiotherapy (as measured by the relative change in the perfusion fraction D % P) based on geometric and topological characteristics of the vasculature on Day 0 (i.e., just before the irradiation).Next, we carefully investigate the relationship between key characteristics of the tumour vasculature on Day 0 and the perfusion response.
Network size (vessel count).In general, the resistance of a vascular network to flow increases with the number of vessels.Therefore, one might naively expect that pruning vessels in a large network should always increase its perfusion and, consequently, that irradiation should improve perfusion at a rate proportional to the number of vessels in the network.Most of our networks are consistent with this principle.However, the PF decreased following irradiation for tumour 4, even though it has a large network (see Table 1).We conclude that the vessel count alone is insufficient to predict perfusion response to radiotherapy.
Geometric determinants of perfusion.We now consider geometric determinants of flow resistance in individual vessels.One typically assumes that the flow rate Q through a cylindrical vessel, subject to a pressure drop Δp across its length L, satisfies Poiseuille's law: where d denotes vessel diameter and μ eff is the effective viscosity of blood [23,52].This parameter depends in a nonlinear way on the vessel diameter and haematocrit (see [53]).Determining the haematocrit distribution within a network is challenging due to uncertainty in the locations and strengths of the inlets, the lack of consensus about the functional form describing haematocrit splitting (this remains an active research area [23,[54][55][56]), and the highly coupled and nonlinear nature of the haematocrit and blood flow.Therefore, for the sake of computational efficiency, we impose a uniform haematocrit H = 0.45 in all networks (the effect of varying H is documented in Section S1.5 of S1 Appendix).Substituting this value into standard formulas obtained in [53] (see Eqs ( 17)-( 19)), we arrive at a relationship describing how the effective blood viscosity depends on vessel diameter (also known as the Fåhraeus-Lindqvist effect) [57,58].Note that for consistency with existing mathematical models [53], all length scales in Eqs ( 17)-( 19) are nondimensionalised with respect to 1 μm.Substituting this relationship into Eq (3), we obtain an explicit expression for the resistance R in a vessel of length L and diameter d.As a further simplification, we propose the quotient R geom ¼ L=d 4 as the simplest proxy for vessel resistance involving only the key geometric parameters.At the network scale, we define the following measures for the total resistance of a network: where the sum is over all network vessels i.Finally, we define the following proxies for mean vessel resistance: Due to the fourth-power dependence of the resistance on diameter, we expect that the average vessel diameter in a network will strongly impact its resistance.Table 1 confirms this: networks with small mean diameters typically have high mean resistance (both R and R geom ).Furthermore, the mean diameter and mean resistances appear to be good predictors of radiotherapy response: low values of mean diameter (high values of mean resistance) correspond to group A, and vice versa.We note that neither the mean length nor the total resistance (which neglects network size, and decreases when vessels are pruned) distinguish between the two groups of tumours.In summary, for the 6 studied tumours, the mean diameter d and mean resistance R distinguish between the two groups, with group A having, on average, thinner (higher resistance) vessels than group B. Topological determinants of perfusion.We use Betti curves [59,60] to track the number of loops in vascular networks during radiotherapy-induced vessel pruning.The Betti numbers β 0 and β 1 refer to the number of connected components and the number of loops in a network, respectively.Via the Euler-Poincare ´formula [60,61], the Euler characteristic X of a network, which is given by its number of edges (vessels) N V and nodes N N , i.e., X ¼ N N À N V , is directly connected to its Betti numbers: X ¼ b 0 À b 1 .Given β 0 , N V , and N N , β 1 can therefore be computed as [60]: To obtain a Betti curve for β 1 , we track the number of loops throughout the pruning process on a network.As in the previous section, we also introduce metrics normalised by the number of vessels, i.e.: We note from Tables 1 and A (the latter in S1 Appendix.)that β 0 , b 0 , and β 1 do not distinguish between the two groups of tumours (A and B), whereas b 1 -which we will call loops per vessel -provides a (weak) distinction between the two groups: tumours with fewer loops per vessel belong to group A. To elucidate why the loops per vessel on Day 0 might impact the perfusion response, we use Eqs ( 6) and ( 7) to write: In what follows, we consider a network containing N I/O/ST nodes of degree 1 (inlets, outlets, and tips of angiogenic sprouts that have not yet anastomosed) and N B nodes of degree 3 (bifurcation points).For simplicity, we ignore nodes of degree 2 (vessels subdivided into segments), 4 (trifurcations), and higher as these seldom appear in our networks.We then have N N = N I/O/ ST + N B and, counting vessels twice by looping over all network nodes, we have N V = (N I/O/ST + 3N B )/2.Using these relations, Eq (8) can be simplified to read: where N I=O=ST denotes the normalised number of degree-1 nodes.We see that small values of b 1 result from large numbers of nodes of degree 1, and vice versa.As the degree of a node relates to its connectivity, vessels adjacent to degree-1 nodes are unlikely to play a significant role in network perfusion and pruning such vessels is desirable.Moreover, while some of these are inlet or outlet nodes, the majority are likely to represent tips of angiogenic sprouts (blunt ends) that have not yet anastomosed.Consequently, the relevant vessels are hypoperfused.Pruning these sprouts increases the perfusion fraction simply by reducing the denominator in Eq (1), whereas pruning vessels that connect bifurcation points (degree = 3) may disconnect parts of the network that were previously connected.This is confirmed in the bottom-left panels in Fig 1 : the vasculature for which the PF increased (tumour 1; bottom) contains more angiogenic sprouts (i.e., is less inter-connected, with lower b 1 ) on Day 0, than the one for which the PF decreased (tumour 6; top).Taken together, these results provide a possible explanation for why tumours with low b 1 increase their PF, and vice versa.

Metrics combining geometry and topology.
A full understanding of the determinants of network perfusion is still lacking.It requires more detailed knowledge of the diameter and length distributions within a network and their connectivity.The intricate (nonlinear) nature of blood rheology further complicates the situation.A single measure (geometric or topological) is unlikely to contain complete information about the perfusion fraction.To assess the impact of radiotherapy on the perfusion fraction, one needs to consider the distribution of high-and low-resistance loops that are being pruned from the network.To this end, and inspired by the observations above, we now propose new metrics that combine network topology with its geometric properties in the form: Table 1 shows that these metrics can distinguish between the two groups of tumours (particularly R b , which accounts for the dependence of viscosity on the vessel diameter).To our knowledge, these are the first global metrics to combine the geometry and topology of vascular networks.
Taken together, our analysis suggests that networks containing large numbers of high-resistance vessels and angiogenic sprouts are likely to increase their PF following radiotherapy.In other words, we postulate that irradiation-induced vessel pruning will increase the PF of vascular networks containing large proportions of thin and blunt-ended vessels-we will next test this hypothesis using synthetic networks.

Model overview
Our analysis of the experiments revealed that vascular architecture influenced the outcome of pruning: tumours with lower mean vessel diameters (d) and fewer vascular loops per vessel (b 1 ) exhibited an increase in perfusion (D % P > 0) when irradiated.We then used our computational model to better understand the causal relationship between the architectural properties of a vascular network and how its perfusion changes following radiotherapy (Fig 4).In brief, we constructed multiple networks to replicate the characteristics of interest (d and b 1 ) exhibited by the biological dataset.We then simulated blood flow through the networks and measured the change in perfusion as vessels were successively pruned.The design of the model is detailed below, with the parameters in Section S1.8 of S1 Appendix.

Network design
We used forking and hexagonal networks to represent hierarchical and non-hierarchical vasculature respectively (Table 3).We then simulated blood flow through the networks.Using simple symmetric geometries allows us to isolate the influence of the different geometric and topological factors (e.g., vessel diameters, lengths, and network loops) under consideration.While the geometries of the networks are simple, they reflect several properties of biological vasculature.Here, we describe the forking networks since they form the focus of our results.The design of the hexagonal network can be found in Section S1.3 of S1 Appendix.In the forking network, each vessel divides into two daughter vessels until 7 generations are created (counting the inlet vessel as generation 0).After the seventh generation (generation 6), the network converges symmetrically into a single outlet vessel.The effect of varying the number of generations is documented in Section S1.6 of S1 Appendix.Note that throughout this work, for simplicity, we consider any straight-line segment to be an individual vessel (as opposed to defining vessels by their two endpoints being bifurcation nodes, inlets, outlets, or blunt ends).
Vessel dimensions.The diameters of the two daughter vessels (d A and d B ) are related to the diameter of the parent vessel (d parent ) via Murray's law [64]: Therefore, the inlet vessel diameter d inlet modulates the thickness of all subsequent vessels.Following [23], we assume that the vessel lengths L are proportional to vessel diameters (d) so that: Based on the mean vessel length to mean diameter ratios in our experimental data, we fix λ = 4.The branching angle at a bifurcation is dictated by the y-extent of its constituent vessels.The y-extent of a vessel in generation i � 1 is the length of its projection onto the y-axis (V i ) in a manner that places the following limit on the spatial extent of the network [23]: We set V 1 = 0.9L 1 for our simulations to allow the network to extend to a sufficient degree along the y-axis, in line with [23].If we further impose that two daughter vessels of any parent vessel must have the same diameter (d A = d B ), then the network geometry is fully specified.This forms our reference network (Fig 5).Our reference network must be modified to reflect the properties of interest (d and b 1 ) in the biological networks.Therefore, we offset the mean vessel diameter of the synthetic networks to match the minimum, maximum, and average values of d observed across the biological vasculatures (see Section S1.2 of S1 Appendix for details on how we did so).Throughout our diameter variations, all vessel lengths and branching angles remain identical to the reference network.Therefore, Eq (12) does not apply in these cases.Variations in d do not result in variations in b 1 within a network.Therefore, we tracked b 1 as we pruned the forking network and measured the corresponding changes in perfusion.In effect, we considered every vessel removal as a way of generating a new network with a different b 1 .Network heterogeneity.To ensure heterogeneous diameters in the network, we fixed the diameter of one of the two daughter vessels (d A ) at each bifurcation to be α-times (α > 1) thicker than the other (d B ): In the homogeneous case (α = 1.0), both daughter vessels have the same diameter.The diameter of a vessel in any generation can be calculated as: where d i is the diameter of a vessel in generation i on a path that features n thick thick vessels originating from an inlet vessel with diameter d inlet .Since every permutation of thick and thin vessels on a path exists in the network, it makes no difference which daughter vessel is α-times thicker than the other at each bifurcation.We choose the vessel that extends upwards from a bifurcation to be thicker than the lower vessel.As a result, vessels at the top of the network are the thickest in their generation and those at the bottom are the thinnest (Fig 5).

Pruning
Radiotherapy preferentially prunes smaller vessels and the extent of pruning increases with the dose [14,15].Therefore, we simulated increasing magnitudes of dosage by removing vessels individually in order of increasing diameter (Fig 6).Recall that the constant λ dictates the ratio between vessel lengths and diameters in the forking network.Therefore, pruning by diameter effectively prunes by length as well, but provides more data points than the 7 distinct lengths present in the forking network.Removing one vessel at a time allows us to examine pruning at a greater temporal resolution than obtainable experimentally.If two vessels have the same diameter then the vessels are pruned in order of their Vessel ID as detailed in Section S1.2 of S1 Appendix.The effect of pruning multiple vessels of equal diameters simultaneously is documented in Section S1.7 of S1 Appendix.Due to the lack of hierarchy, isolated vessels may remain in the hexagonal network during the course of pruning.Note that our model can only be used to mimic the effect of single-dose radiotherapy over relatively short time scales, as a reliable model of a fractionated radiotherapy over the course of multiple weeks would require one to take into account significant angiogenic growth and vascular remodelling.

Perfusion threshold
Instead of identifying vessels as perfused based on fluorescence, we track the flow directly and introduce the perfusion threshold (Q min ) as the minimum flow rate at, or above, which vessels are considered perfused and below which vessels are considered 'hypoperfused'.This threshold serves as a proxy for the sensitivity with which our experimental apparatus can detect the signal coming from perfused vessels.In the forking network, we set the value of this threshold Q min to be 3 × 10 −12 m 3 s −1 which, for the default model parameters and the three offsetting scenarios of the least heterogeneous network (α = 1.0), yields initial PFs of 0.24 (minimum), 0.50 (average), and 1.00 (maximum).In this way, we cover a large range of initial PFs spanning almost the entire interval (0, 1) which was also observed in real vasculatures (initial PFs ranging from 0.14 to 0.97).Similarly, we set Q min = 3 × 10 −13 m 3 s −1 for the hexagonal network.This value allows us to cover a sufficiently large range of initial PFs.

Perfusion fraction
We compute the PF using the experimental formula (Eq (1)).Guided by Eq (2), we quantify improvements in perfusion for synthetic networks by introducing D max P, the maximum value attained by DP during pruning, and D max % P, the maximum percentage change in the PF, where: D max P ¼ P max À P 0 and D max % P ¼ In Eq (16), P 0 and P max denote, respectively, the initial and maximum value of the PF during pruning.

Results
In line with our analysis of actual tumours, we found that the largest increases in PF in our synthetic vascular networks were associated with lower values of d and b 1 .Additionally, we identified two mechanisms that can effect a positive change in perfusion post-irradiation.The vascular remodelling that occurs in vivo makes it difficult to isolate the effects of each mechanism.By contrast, synthetic networks offer no such barrier to analysis.We detail our inferences from these simulated networks in the following subsections.
While our forking networks are inherently hierarchical, the same (high) level of order might not be present in dysregulated tumour vasculatures.In order to gain insight into the role that hierarchy plays, we also examined the effect of non-hierarchical structures modelled as hexagonal networks and the impact of diameter heterogeneity.We found that hierarchy was conducive to perfusion enhancement and that heterogeneity resulted in an increased and sustained D % P response.For a full discussion, see Section S1.4 of S1 Appendix.

Two mechanisms can increase the perfusion fraction post-irradiation
One can deduce from Eq (1) that the perfusion fraction (P) can increase not only if the number of perfused vessels increases after pruning, but also if the number of hypoperfused vessels decreases ( Fig 7).We conclude that two mechanisms can increase the value of D % P when a hypoperfused vessel is pruned: • Mechanism 1: When a vessel is pruned, blood flow may be rerouted causing one, or more, hypoperfused vessels to become perfused and increasing the number of perfused vessels (Fig 8).• Mechanism 2: Removing hypoperfused vessels increases the proportion, but not the number, of perfused vessels.
Fig 9 shows how both mechanisms act on a network during pruning.In the original, unpruned network, several hypoperfused vessels have blood flow rates that are close to the perfusion threshold while others have flow rates that are well below the perfusion threshold.When the least-perfused vessels are pruned, their flow is rerouted and the flow rate in the remaining vessels (several of which were previously hypoperfused) increases above the perfusion threshold.A decrease in the number of hypoperfused vessels alone would correspond to a network that becomes more efficient following radiotherapy.However, it would not necessarily correspond to an improvement in oxygen distribution.

Mean vascular diameter determines the relative contribution of the two perfusion improvement mechanisms
Having identified two mechanisms that can increase the perfusion fraction, we next examined how the mean vascular diameter in the unpruned networks (d) affects the percentage change in the perfusion fraction (D % P).As in the biological experiments, we found that lower values of d generated larger increases in D % P and higher values of D max % P (Fig 10).However, we also examined the two mechanisms of improvement in isolation and found that higher values of d were conducive to flow rerouting, while lower values tended to improve D % P through Mechanism 2.
A low d results in an increase in D % P for two reasons.Firstly, thinner networks have lower starting perfusion fractions (P 0 ) than thicker networks ( Fig 11).Therefore, any change in P is greater relative to the network's initial P 0 .P 0 is lower in thinner networks because they offer greater resistance to blood flow (Fig 12).Secondly, thinner networks have a greater proportion of hypoperfused vessels that can be pruned before any perfused vessels are pruned.In these networks, P can increase a greater deal as a result of pruning only hypoperfused vessels.
On the other hand, vascular architectures with a high d are more susceptible to increases in D % P through blood flow rerouting ( Fig 10).For these architectures, flow rates in the pruned (hypoperfused) vessels are closer to the perfusion threshold.Therefore, these vessels require a smaller increase in blood supply (via rerouting) to turn them into perfused vessels.These vessels also provide less resistance to the blood that is rerouted into them once pruning begins.

A drop in the proportion of loops precedes a rise in perfusion
Thus far, we have compared networks that differ in terms of their initial, unpruned architecture.Since pruning produces blunt ends (thereby changing the number of loops as discussed in the section Topological determinants of perfusion), the pruned networks can serve as Considering the central portion of this network, panels (B)-(G) present the progression of flow distributions in the process of pruning, noting that the maximum of the colourmap is set equal to the perfusion threshold (i.e., vessels that are not dark red are hypoperfused).Having pruned (B) the 20 thinnest vessels, the network contains pairs of daughter vessels which experience flow rates slightly below the perfusion threshold (Q min = 3 × 10 −12 m 3 s −1 ).Pruning one vessel of such a pair is likely to result not only in (C) fewer hypoperfused vessels but also in more perfused vessels due to flow rerouting.Between (D) and (E) and between (F) and (G), we primarily prune proxies for networks that initially contain different proportions of angiogenic sprouts.Contrasting the number of loops per vessel (b 1 ) with the change in perfusion fraction (D % P), we observed two distinct results (Fig 13).
In the early stage of pruning (less radiation-induced cell death), D % P typically increases regardless of the behaviour of b 1 .In this stage, pruning removes the thinnest vessels, which tend to constitute hypoperfused loops or blunt ends.Thus, D % P increases through both blunt ends which increases the perfusion fraction (D % P) via Mechanism 2. Between (E) and (F), D % P dramatically decreases because perfused paths are disrupted.https://doi.org/10.1371/journal.pcbi.1011252.g009

Fig 10. The impact of the mean diameter and the diameter heterogeneity on the perfusion response. Series of plots (A-I)
showing, for the forking network, the influence of the mean diameter (d) and the diameter heterogeneity (α) of the unpruned network on how perfusion (D % P) changes during pruning.An increase in the number of perfused vessels is evidence of flow being rerouted into vessels that were previously hypoperfused.An increase in D % P without an increase in perfused vessels is evidence of Mechanism 2. Regions in which Mechanism 1 is primarily active are highlighted in grey.https://doi.org/10.1371/journal.pcbi.1011252.g010mechanisms of improvement.In the later stages of pruning (more radiation-induced cell death), however, peaks and troughs in b 1 correspond to those in D % P. In particular, a reduction in b 1 often precedes an increase in D % P. When a vessel that separates two loops is pruned, the loops merge into one and the total number of loops decreases.At the site of the pruned vessel, a blunt end may be created.Subsequent pruning favours the removal of this blunt end due to its small diameter by virtue of previously being connected to the smallest vessel (recall Murray's law, Eq (11)).Pruning the blunt-ended vessel does not change the total number of loops.However, it increases the proportion of loops (b 1 ) and the resulting reduction in the hypoperfused vessel count leads to an increase in D % P. Therefore, a decrease in b 1 indicates that the network now contains a greater proportion of blunt ends (and vice versa).Thus, we find that D % P increases when networks with a low value of b 1 are pruned.We also infer that b 1 indicates when an improvement seen in D % P is caused by the pruning of blunt-ended vessels with no flow, rather than by rerouting.

Discussion
In this paper, we used a computational model to investigate the influence of vascular architecture on the outcome of radiation-induced pruning.As per [15], we measured these outcomes in terms of overall network perfusion, as defined by the PF (Eq (1)).In our experimental data, we found that perfusion typically improved in vascular networks with lower mean diameters (d) and fewer loops per vessel (b 1 ).It may seem counter-intuitive for network perfusion to increase after vessels are removed in a model that excludes angiogenesis.However, our synthetic networks mirrored our observations of the experimental data.We modelled the influence of radiotherapy by assuming that pruning proceeds from the smallest to largest diameter vessels.Using our in silico model, we elucidated the mechanisms underlying the perfusion response, including two mechanisms of improvement in the PF.We have also shown that different architectural features are susceptible to different mechanisms of improvement in the PF.In particular, we found that networks with low average diameters tend to exhibit improvements in PF via the reduction of hypoperfused vessels, while networks with large average diameters are more prone to rerouting (see the section Mean vascular diameter determines the relative contribution of the two perfusion improvement mechanisms).In the former, large doses of radiotherapy (i.e., more vessels pruned) seem to improve the PF, while rerouting seems to occur in the latter with smaller doses of radiotherapy ( Fig 10).
Therefore, our model represents a first step towards understanding the different perfusion responses found in the literature, where some studies reported increased perfusion post-irradiation and some did not [19].Moreover, our exposition of the different mechanisms of PF improvement may explain why Bussink et al [22] observed an increase followed by a decrease in tumour perfusion.Bussink et al. attributed the increase in perfusion to the redistribution of blood and the refilling of previously non-functional vessels (in line with our model's predictions).
While we are not aware of any experimentally informed computational models of radiotherapy-induced pruning, several studies have modelled the effect of vessel pruning due to low wall shear stress.For example, when modelling angiogenesis, Owen et al. found that poorlyperfused vessels tended to be pruned and that new vessels were formed by 'stealing' flow from other vessels [26].Although we pruned our synthetic networks in order of increasing vessel diameter, we note that the thinnest vessels are likely to have the lowest blood flow.Further, our models also exhibited flow rerouting, with previously hypoperfused vessels becoming perfused following vessel pruning.Ste ´phanou et al. simulated vessel pruning based on flow rates to investigate the rate of drug uptake [65].Pruning vessels randomly sometimes resulted in improved rerouting of blood flow, and increased drug uptake [65].Pruning low-flow vessels, on the other hand, did not improve drug uptake but improved delivery speed as blunt-ended vessels were removed [65].In our study, the removal of blunt ends is also unlikely to result in improved oxygenation, although a full investigation is reserved for future work.

Experimental implications
Our findings also have several implications for experimental methodologies used to study radiation-induced pruning.For one, our results in the section Two mechanisms can increase the perfusion fraction post-irradiation have shown that the PF is not an ideal measure of a tumour's perfusion, since it fails to account for hypoperfused vessels that are pruned.Furthermore, the PF offers no indication of a vascular network's ability to sufficiently deliver oxygen in tissue.Theoretically, a region could feature a single perfused vessel and have a PF of 1 despite the fact that the vessel may not cover a sufficient area or carry enough red blood cells to maintain a normoxic environment.Experimentally, Bussink et al. found discrepancies between perfusion and hypoxia [22].Monitoring the PF alone cannot tell us which mechanism (as discussed in the section Two mechanisms can increase the perfusion fraction post-irradiation) is at work during pruning.Discerning the mechanism is important because it tells us whether or not blood flow has improved in the remaining vessels.Future experiments could be used to validate the model's findings.For example, if a vascular network has a low average diameter and a low number of loops per vessel, our study suggests that the network is likely to show an improvement in perfusion after irradiation.To avoid the pitfalls of using the perfusion fraction metric, it would be useful to measure flow as total blood flow rate or in terms of the tissue area under normoxic oxygen tension.
In the section A drop in the proportion of loops precedes a rise in perfusion, we showed that the number of loops per vessel may act as an indicator of improvements in PF via the reduction of hypoperfused vessels.In Table 1, we classified biological tumours into groups A and B based on the DP on a fixed day.However, the non-monotonic response of P in our simulations implies that this classification is subject to the sampling point (Table 4).In other words, a network belongs to group A or B depending on the number of vessels pruned from it.Interestingly, for any fixed α, the lower the mean diameter, the more often one observes an increase in PF.Similar correlations were found in real vasculatures.
We also computed for forking networks a compound metric discussed in the section Metrics combining geometry and topology: the resistance per loop.As expected from Eq (10), over the course of pruning, the trend of R geom b follows that of b 1 in an inverse fashion (Figs 13   and 14).However, R geom b also encodes information on the mean diameter (via resistance) and thus allows us to distinguish between networks of varying mean diameters.Consistent with real vasculatures, we observe that networks with initially higher R geom b exhibit greater potential for perfusion improvement, as measured by D max % P (compare with Fig 10).This leads us to propose the metric R geom b for future investigations as a predictor of perfusion improvement, useful especially when d and b 1 individually give contradictory predictions.This might be relevant, for example, for vasculatures with low mean diameters and a large proportion of loops, even though such a combination of characteristics did not occur within the small tumour sample studied here.These speculations need to be confirmed by future work and the mechanism of improvement must also be discerned.
While the focus of this study is the generation of mechanistic insights into enhanced perfusion, we note that the clinical implications of our model are limited by current vascular imaging technologies.High-resolution methods, such as the multiphoton microscopy used in this study, are limited to a penetration depth of approximately 1 mm, requiring the use of a window chamber [66].In a clinical setting, vascularised tumours may be located below this depth.On the other hand, techniques such as MRI and CT have previously been used to study vascular function in humans but are limited to a resolution of 100 microns, while our model shows that the redirection of blood flow is prevalent in smaller vessels [67].Thus, our predictions from synthetic networks would require clinical imaging technology to advance considerably before their prognostic value in humans can be evaluated.

Future work
In this work, we used forking and hexagonal networks to illustrate the role of hierarchy in vascular architecture.However, both network types are regular (i.e., symmetric) and might not adequately represent all features of dysregulated tumour vasculature.As a next step, one should consider irregular networks, such as those generated by the Voronoi tessellation of a planar region [68,69].Tumour vasculature is further known to contain many leaky and tortuous vessels which will also need to be accounted for in the model.
We performed our blood flow simulations under the assumption that the haematocrit (fraction of blood made up of red blood cells) distributes itself evenly in all vessels.In reality, however, the spread of haematocrit at vessel junctions is more complex.The haematocrit in each vessel affects its flow rate (Eq (17)), as well as its capacity to carry oxygen.Additionally, our colleagues have hypothesised that irradiating a tumour with a larger proportion of small vessels might be more likely to lead to improved perfusion by altering the proportion of haematocrit splitting [15].Therefore, it will be important to contrast different haematocrit splitting rules sourced from the literature [23,70,71].Different splitting rules will give rise to different oxygen distributions in the tissue.These distributions are a key determinant of the efficacy of standard oncological treatments, such as chemotherapy, which can also be modelled within Microvessel Chaste [62,63] as a chemical transport process.The death of tumour cells caused by radiotherapy or chemotherapy lowers net oxygen consumption, which further increases tissue oxygen levels.The inclusion of an oxygen consumption term proportional to the local density of live tissue cells will be an important next step in the precise modelling of tissue oxygen perfusion.We have already highlighted that the PF is an imprecise metric to describe the oxygenation status of the tissue.Once we explicitly simulate oxygen diffusion, we will need to adopt and develop new spatial metrics that better describe the spread of oxygen.With a model of oxygenation in place, we will then consider the dynamic nature of the tissue perfused by the vasculature.Cells that comprise the tissue will have varying uptake rates, and will divide and die based on the amount of oxygen available to them.
Pruning can also be conducted stochastically, provided a reliable relationship between the vessel diameter and the probability of its pruning is found, and the pruning rules can be modified to reflect radioresistance when the region is anoxic.Let us note that the existing model can also be adapted to better represent the regression of hypoperfused vessels as observed in developmental vascular networks by pruning vessels in order of increasing flow rate [17,49].Since the vasculature itself is also dynamic, future models must factor in structural adaptation and the secretion of angiogenic factors.Random sprouting models of tumour angiogenesis may be useful here [72].Percolation models would also serve as a good representation of the statistical distribution of avascular spaces found in tumours [73].Moreover, percolation models have been used to replicate observed oxygenation conditions in mice tumour xenografts [73].Such a study may also involve modifying haematocrit splitting rules to include higherorder splitting, since tumour vasculature is not limited to bifurcations [23].
Finally, we may leverage our insights from synthetic networks and apply them to a real biological network obtained via the methods outlined in [35].Provided that inlet and outlet vessels with appropriate boundary conditions are identifiable, we may then simulate blood flow and the pruning of vessels in image-derived networks.However, the implications of this study remain limited due to a relatively small experimental dataset associated with a single tumour cell line.Trends identified in this work will thus need to be tested using different cancers and biological models.

Experimental procedures and data preprocessing
Ethics statement.All animal experiments were conducted in accordance with the United Kingdom Animals (Scientific Procedures) Act 1986 as amended (Amendment Regulations 2012 [SI 2012/3039]), under the authority of a UK Home Office Project License (PPL 30/2922 and PCDCAFDE0), with local ethical approval from the University of Oxford Animal Welfare and Ethical Review Panel, as previously described [15].
Experimental procedures.Abdominal imaging window implantation.This procedure was based on a previously described method [74].Transgenic mice with fluorescent protein tdTomato expressed in endothelial cells on C57Bl/6 background were prepared in a surgical unit, and administered with inhalational anaesthesia and pre-operative analgesics.Body temperature and respiration rates were monitored throughout the procedure.A 1 cm cut was made along the abdominal midline approximately 5 mm underneath the sternum followed by blunt dissection around the cut to separate the connective tissue from the skin.A custommade imaging window frame (Workshop at the Department of Oncology, Oxford University) was fitted underneath the skin.Continuous sutures were used to secure the skin to the window frame.Approximately 2.5 × 10 5 MC38 cells (murine colon adenocarcinoma cells) stably expressing eGFP in 5 μL containing 30% of Matrigel and 10% of Evan's blue dye were injected under the connective tissue and above the abdominal muscle layer.The chamber was then flushed with water to lyse non-injected cells by osmotic shock, tapped dry with sterile cotton swabs and flooded with saline.A cover glass glued on the chamber's lid was secured onto the window frame.The animals were then placed onto a heat mat for post-operative recovery, and their health and tumour growth was monitored by visual examination.
Treatment regimes.Animals with tumours approximately 100 mm 3 growing in the window chamber were administered radiation treatment.Prior to the radiation treatment, mice were anaesthetised under inhalation with isoflurane and placed in an imaging-guided small animal radiation research platform (SARRP) irradiator (Xstrahl Ltd).A Cone Beam CT scan (computerised tomography) of each mouse was obtained and the treatment was planned using Muriplan (Xstrahl Ltd).The SARRP was used to deliver 15 Gy of X-rays (220 kVp copper filtered beam with HVL of 0.93 mmCu) to the tumour at 2 Gy per minute.Dosimetry of the irradiator was performed as previously described [75].We refer to the start of treatment as Day 0.
Intravital two-photon imaging.Mice were imaged for 7 days after radiation treatment with a Zeiss LSM 880 microscope equipped with a respiratory monitoring system.The stage and atmosphere were heated to 37˚C.To label perfused vessels, Qtracker705 Vascular Labels (0.2 μM, ThermoFisher Scientific) were infused intravenously using a motorised pump at a rate of 0.84 μL�min −1 .A mode-locked MaiTai laser tuned to 920 nm was used to simultaneously excite eGFP, tdTomato, and Qtracker705.The Qtracker705 signal was acquired through a BP700/100 filter with a non-descanned detector.GaAsP detectors were used to acquire the signal of tdTomato selected by a BP 650/45 filter and the eGFP selected by a BP525/50 filter.Images were acquired in Z-stack tile scans with a pixel size of 0.823 μm and an image size per tile of 512 × 512 × 5 in x, y, and z, respectively.A water immersion 20× objective made for UV-VIS-IR transmission with a numerical aperture of 1.0 was used.Representative examples of vascular networks from our experiments are displayed in Fig 2 .Data preprocessing.The biological networks were obtained by multiphoton intravital 3D imaging [76] and consisted of 3D stacks of images of tumour blood vessels.Skeleton files were extracted from the imaging data by combining two segmentation models and taking their geometric average.The skeletons were then pruned (see reference [77], p. 165, for a full description).We extracted blood vessel networks from skeleton files using the method VesselTree from UNET_CORE.VESSEL_ANALYSIS in the Python code package UNET-CORE [78].The extracted networks consist of points on vessel branches (multiple points per vessel branch including branching points) which represent the network nodes, and the vessels that connect them which constitute the edges of the network.VesselTree also enables us to extract network features such as vessel diameters and lengths.

Simulation methods
Blood flow.We assume that the blood flow rate Q in a vessel of length L and diameter d is determined by Poiseuille's law (Eq (3)).Following [23,53], the effective blood viscosity in Eq (3) is expressed as: " " where H is the discharge haematocrit, d is the vessel diameter, μ p is the plasma viscosity, m 45 ¼ 6e À 0:085d þ 3:2 À 2:44e À 0:06d 0:645 ; ð18Þ and: As discussed in the section Geometric determinants of perfusion, we assume that haematocrit is distributed evenly within the network (H = H inlet = 0.45) to yield a linear flow problem.With signed flow rates Q i for each vessel i (i.e., each edge connected to a node), we also impose conservation of blood at each network bifurcation node: In each (unpruned) network, the leftmost nodes serve as inlets and all other vessels with one detached node serve as outlets.Inlet nodes are assigned a pressure (p inlet ) of 3333 Pa (� 25 mmHg), while outlet nodes are assigned a pressure (p outlet ) of 2000 Pa (� 15 mmHg) [79].Software development and post-processing.We performed our simulations using Microvessel Chaste (version 3.4.3),an add-on to the open-source simulation package Chaste (version 2020.1)[62,63].We added custom functionality to the base version of Microvessel Chaste, including new network generators and pruning functions.We generated synthetic networks detailed above and pruned them vessel-by-vessel in order of increasing vessel diameter.At each pruning step, we updated the flow distribution, and exported the perfusion fraction and the network itself in the .vtkformat, accounting for all vessels that had not yet been pruned.We used Python scripts to generate from the .vtkfiles adjacency matrices weighted by vessel diameter and vessel length, from which we calculated key geometric metrics, such as mean diameter and mean geometric resistance.Based on a MATLAB package [80], we wrote Python scripts to import the adjacency matrices from above, calculate both the number of loops and loops per vessel at each pruning stage using Eq (6), and finally divide the mean resistance by the latter.

Fig 1 .
Fig 1. Chart summarizing the key components of the present study.Bottom-left panels present representative vascular regions containing perfused (yellow) and hypoperfused (purple) vessels from Day 0 (just before irradiation) for one of the tumours for which irradiationinduced vessel pruning led to a decrease in perfusion (tumour 6; top) and one for which it led to an increase (tumour 1; bottom).Note that the latter vasculature contained many hypoperfused blunt-ended vessels, while the former contained more loops.Bottom-right panels show pruned synthetic (forking) networks exhibiting similar properties with perfused (red) and hypoperfused (blue) vessels.

Fig 2 .Fig 3 .
Fig 2. Representative vascular architectures from our experiments[15].A single dose of 15 Gy X-rays was delivered to a tumour and changes to its vascular structure were monitored over the course of 4 days (reproduced from[50]).In the top row, endothelial cells are in cyan, while perfused vessels (qDot705) are in red.Every panel is a 2D representation of a 3D image in the form of a Z-stack approximately 300 μm tall.The depth of field for a single Z-slice was 2 μm.The middle and the bottom row show images of endothelial cells and perfused vessels, respectively.The scale bar corresponds to 100 μm.https://doi.org/10.1371/journal.pcbi.1011252.g002

Fig 5 .
Fig 5. Blood flows from left to right through a single inlet and a single outlet in the forking network.L i and V i denote the length and y-extent of a vessel in generation i, respectively.Note that our implementation of network heterogeneity results in vessels at the top of the network being the thickest in their generation and vessels at the bottom being the thinnest.https://doi.org/10.1371/journal.pcbi.1011252.g005

Fig 7 .
Fig 7. The two mechanisms of P increase.When a forking network (α = 1.1, d ¼ 28:50 μm) is pruned, perfusion (D % P) can increase, shown in (A), via Mechanism I when blood flow is rerouted from pruned vessels to other vessels, as seen in (B) where the number of perfused vessels increases and the number of hypoperfused vessels decreases.In contrast, perfusion (D % P) can also increase, shown in (A), via Mechanism II when hypoperfused vessels are pruned without any net increase in the number of perfused vessels, as evidenced in (B) by the constant number of perfused vessels.https://doi.org/10.1371/journal.pcbi.1011252.g007

Fig 11 .Fig 12 .
Fig 11.Perfusion response as measured by P. The initial perfusion fraction (P 0 ) depends on the vascular architecture (A-C).Note that, for a given diameter heterogeneity (α), the perfusion fraction (P) of all mean diameters (d) converges during the later stages of pruning, when the networks only contain thicker (low-resistance) vessels.Panel (D) shows the pruned architecture at a stage common to all values of d for α = 1.2 (147 pruned vessels).Here, we observe that all remaining vessels are either perfused or have blunt ends (no flow).No flow rerouting can occur in the hypoperfused vessels.Thus, subsequent pruning produces the same changes in P, regardless of the initial d.Note that the vessel colour indicates the local flow rate, with the range of the colourbar adjusted so that perfused vessels are dark red and vessels with little or no flow are dark blue.https://doi.org/10.1371/journal.pcbi.1011252.g011

Fig 13 .
Fig 13.The proportion of loops.A drop in the number of loops per vessel (b 1 ) precedes a positive change in the perfusion fraction (D % P) in the forking networks (d ¼ 33:65 μm) across all values of α (A-C).Minor differences in b 1 arise between heterogeneities because α changes the order of pruning by virtue of changing the diameters of vessels.https://doi.org/10.1371/journal.pcbi.1011252.g013

Fig 14 .
Fig 14.Resistance per loop.The resistance per loop (R geom b ) follows trends inverse to those observed for the loops per vessel (b 1 ) for all values of α (A-C), but can also distinguish between networks of varying mean diameter.https://doi.org/10.1371/journal.pcbi.1011252.g014

Table 2 . Perfusion fractions on Day 0 and the final day of the pruning phase, the pruning-induced perfusion dif- ference, and its relative counterpart for six studied tumours.
Tumours are ordered based on the relative change in PF as measured by D % P.

Table 4 . The classification of forking networks into groups A and B varies based on the sampling point, i.e., how many vessels have been pruned when
D % P is measured.Cells have been shaded green or red to represent positive or negative values of D % P.