Spatial Metrics of Tumour Vascular Organisation Predict Radiation Efficacy in a Computational Model

Intratumoural heterogeneity is known to contribute to poor therapeutic response. Variations in oxygen tension in particular have been correlated with changes in radiation response in vitro and at the clinical scale with overall survival. Heterogeneity at the microscopic scale in tumour blood vessel architecture has been described, and is one source of the underlying variations in oxygen tension. We seek to determine whether histologic scale measures of the erratic distribution of blood vessels within a tumour can be used to predict differing radiation response. Using a two-dimensional hybrid cellular automaton model of tumour growth, we evaluate the effect of vessel distribution on cell survival outcomes of simulated radiation therapy. Using the standard equations for the oxygen enhancement ratio for cell survival probability under differing oxygen tensions, we calculate average radiation effect over a range of different vessel densities and organisations. We go on to quantify the vessel distribution heterogeneity and measure spatial organization using Ripley’s L function, a measure designed to detect deviations from complete spatial randomness. We find that under differing regimes of vessel density the correlation coefficient between the measure of spatial organization and radiation effect changes sign. This provides not only a useful way to understand the differences seen in radiation effect for tissues based on vessel architecture, but also an alternate explanation for the vessel normalization hypothesis.

Intratumoural heterogeneity is known to contribute to poor therapeutic response. Variations in oxygen tension in particular have been correlated with changes in radiation response in vitro and at the clinical scale with overall survival. Heterogeneity at the microscopic scale in tumour blood vessel architecture has been described, and is one source of the underlying variations in oxygen tension. We seek to determine whether histologic scale measures of the erratic distribution of blood vessels within a tumour can be used to predict differing radiation response. Using a two-dimensional hybrid cellular automaton model of tumour growth, we evaluate the effect of vessel distribution on cell survival outcomes of simulated radiation therapy. Using the standard equations for the oxygen enhancement ratio for cell survival probability under differing oxygen tensions, we calculate average radiation effect over a range of different vessel densities and organisations. We go on to quantify the vessel distribution heterogeneity and measure spatial organization using Ripley's L function, a measure designed to detect deviations from complete spatial randomness. We find that under differing regimes of vessel density the correlation coefficient between the measure of spatial organization and radiation effect changes sign. This provides not only a useful way to understand the differences seen in radiation effect for tissues based on vessel architecture, but also an alternate explanation for the vessel normalization hypothesis.

Author Summary
In this paper we use a mathematical model, called a hybrid cellular automaton, to study the effect of different vessel distributions on radiation therapy outcomes at the cellular level. We show that the correlation between radiation outcome and spatial organization of vessels changes signs between relatively low and high vessel density. Specifically, that for relatively low vessel density, radiation efficacy is decreased when vessels are more Introduction It is increasingly recognised that an important aspect of cancers is their heterogeneity [1]. This heterogeneity exists between patients, between different tumours within a single patient [2], within the differing cellular populations in a single tumour and even at the genetic scale between cancer cells originating from the same ancestor [3]. In particular, microenvironmental heterogeneity is becoming widely accepted as a key factor in tumour progression and response to therapy [1]. Nutrients, growth factors, extracellular matrix and other cell types are all part of the normal tissue that surrounds and pervades a solid tumour and has been shown to vary widely across different tumour stages and types. This is, in part, due to the dynamic and heterogeneous interplay between the tumour and its microenvironment.
Radiation biologists have, for many years, understood the importance of cell biological and microenvironmental factors on radiation response. Current radiation therapy dose planning, however, largely neglects this information and is, instead, based on years of clinical experience using intuition and trial and error. As such, there remains limited tailoring of dose planning to an individual patient's tumour. With the advent of modern quantitative histologic [4] and biological imaging methods [5], however, this paradigm is poised to change.
Research in this area over the last decade [6] has sought to understand the macroscopic spatial distribution of hypoxia within tumours using non-invasive imaging. This information has then been utilised to develop spatially heterogeneous dose plans to improve tumour control. For example, Malinen et al. [7] inferred average oxygen concentrations from radiocontrast concentrations measured by Dynamic Contrast Enhanced (DCE) Magnetic Resonance Imaging (MRI) in a dog sarcoma. Other work to understand the effects of radiation in individual patients has utilized MRI scans in combination with mathematical models of tumour growth. These models have incorporated heterogeneity in cell type by considering a two compartment spatial partial differential equation (PDE) model, separately accounting for proliferation and motility, without consideration of oxygen effects [8] to explain different radiation responses measured by changes in tumour size over time. More recently, cellular automaton models of stem cell driven tumours, comprised of populations with differences in proliferative phenotype [9] were used to compare tumour response to a range of spatially heterogeneous radiation dose plans or to different sequencing of radio-and chemo-therapeutic strategies [10]. What is lacking to date, however, is research into how tissue level microenvironmental heterogeneity can affect radiation response, and how this could be inferred from patient data.
Heterogeneity in tumour oxygenation, in particular the occurrence of hypoxia, is a wellknown cause of radiation therapy failure [11,12]. Work done in vitro to understand differences in radiation effect due to oxygenation differences have been valuable, and have established an empiric relationship called the Oxygen Enhancement Ratio (OER) [13] to understand how radiation efficacy varies with oxygen concentration. These studies do not, however, allow us to understand the effects of radiation in vivo, as they do not consider the heterogeneity of oxygenation at the microscopic, cellular scale.
Beyond the macroscopic changes in oxygenation, vessel and cellular density and metabolism, a number of theoretical studies have suggested that the local microscopic heterogeneity in oxygenation can vary widely in space and time in tumours [14][15][16] and healthy tissues [17] alike. In addition a large body of work has sought to understand an apparent paradox of therapy directed at angiogenesis, the process of new vessel creation [18]. In short, it was thought that by blocking a cancer's ability to create new vessels, it would be possible to starve the cancer of nutrients and oxygen, quickly leading to its demise. While effective anti-angiogenesis drugs have been developed, this promise never came to fruition. The leading hypothesis to explain this failure is termed the 'vascular normalization hypothesis' [19], which suggests that these drugs, which block vascular endothelial growth factor (VEGF), do not simply inhibit new vessel production, but instead work by normalizing the vascular bed in question by pruning out ineffective vessels and enhancing flow in others. This hypothesis, while well supported by experimental work, has not yet been able to fully explain results at the clinical scale [20]. These studies highlight an opportunity to improve our understanding of how spatio-temporal oxygen dynamics at the cellular scale affect tissue-level response to therapy. To this end, we develop a computational, hybrid cellular automaton model of a tumour growing within a surrounding normal tissue in a vascularized domain which we use to investigate whether spatial statistics gleaned from measures of vascular organisation can be used to predict radiation efficacy. By identifying broad relationships between in silico tissue architecture and radiation response, we aim to progress toward a translatable method of radiation plan optimization using information extracted from biopsies.
The remainder of this paper is structured as follows. We first elucidate our methods by describing the underlying rules and parameters governing the cellular automaton as well as the method of calculating oxygen transport and uptake. In the results section we describe simulation results concerning healthy tissue growth in regularly arranged and then heterogeneous vascular architectures. We then describe tumour growth and invasion in regular architectures and follow this with observations concerning the distribution of oxygen tension in different possible vascular organisations. We then develop a metric, based on Ripley's L function [21], which allows us to quantify and correlate these patterns with radiation response. We end with a discussion of how this metric reconciles some of the difficulties with the vascular normalisation hypothesis, and suggest how the metric might be used in the clinic to personalise radiation dose planning.

Methods
We consider the effect of a heterogeneous microenvironment through the inclusion of en face blood vessels modelled as point sources of oxygen, and through competition of tumour cells with an initial field of healthy cells. While the assumption that blood vessels are point sources, and lie en face may be a simplification of the true biological complexity, it serves as a reasonable, and commonly used [15,16,22], starting point for modelling. We begin by creating a hybrid cellular automaton (HCA) model [23] in which we describe cells by individual agents whose states are updated over synchronous discrete time steps of fixed duration on the timescale of the cell cycle (automaton time steps), and which occupy sites on a two-dimensional square lattice representing a slice of tissue. The size of the lattice spacing is chosen such that each automaton element is approximately the size of a single cell. A second, identical lattice is created on which we approximate the continuous concentration of a freely diffusible molecule, representing oxygen which is updated on a finer timescale (oxygen time step-see supplemental information for further description of the relationship between these timescales). While the cells and oxygen distribution are updated on separate lattices, each influences the other. The feedback between the cells and the microenvironment is captured through a partial differential equation (PDE) governing oxygen transport and consumption. Although the parameter values in our model are drawn from data on a particular cancer type, the primary brain tumour glioblastoma, most of the underlying model assumptions are likely to apply to many other vascularised solid tumours.

Cellular automaton model of cell behaviour
Cell fate decisions in our model are determined by a number of microenvironmental and cell state-specific thresholds and values. The order in which cells are chosen to decide their fate is computed in a random fashion so as to avoid any order bias. To determine the position of cells within our domain, we tessellate the continuous domain on which the PDE is defined into squares of size Δx × Δx to arrive at a regular lattice occupied by both cells and vessels. While all cells are assumed to be the same size and shape (a single lattice site), we model two cell types, cancerous and healthy cells, each of which can divide to produce two identical daughters. The probabilities and thresholds for these cell fate decisions are listed in Table 1 and Fig 1, respectively, and will be discussed individually in the coming sections.
Oxygen consumption and transport. With the advent of single-cell techniques for measuring oxygen consumption rates, quantitative data concerning the differential oxygen consumption rates of different types of cancer cells are now available. It is widely accepted that, in general, cancer cells are more metabolically active than healthy cells because of the Warburg shift [24] or other metabolic abnormalities, and we thus use a modification of normal consumption. The continuous component of our model describes the distribution and consumption of nutrients. While it is clear that many nutrients are of biological importance, in this model we focus on oxygen as it is central to the DNA-damaging effects of photon radiation therapy [25]. Blood vessels, which are each modeled as point sources of oxygen, occupying one lattice site, are placed randomly throughout the lattice at the start of a given simulation, at a specified spatial density Θ. While in reality, vessels can be different sizes, they tend to be normally distributed around a mean size, and a variety of metrics including number of vessels, total vessel area and vascular fraction have been used interchangeably [26,27]. As we are interested in the dynamics over the time scale of radiation treatment (order days), in all simulations we neglect vascular remodeling and angiogenesis, processes which are triggered by hypoxia and mediated by hypoxia inducible factor-1α (HIF-1α) and VEGF, among other molecules, and have been well studied [28,29]. Each vessel is assumed to carry an amount of oxygen equal to that carried in the arterial blood [30]. This oxygen is then assumed to diffuse into the surrounding tissue and be consumed by normal and tumour cells. The spatiotemporal evolution of the oxygen field is described by the reaction-diffusion PDE @cðx; tÞ @t ¼ D c r 2 cðx; tÞ À f c ðx; tÞ; ð1Þ where c(x, t) is the concentration of oxygen at a given time t ! 0 and position x 2 (0, NΔx) 2 nV, where Δx is the length of a typical cell, N is the number of lattice sites on a side of the N × N domain and V = {x 1 ,x 2 , . . ., x v } denotes the set of points occupied by blood vessels, which are v in number. We define D c to be the diffusion coefficient of oxygen, which for simplicity we assume to be constant, corresponding to linear, isotropic diffusion. The function f c (x, t) denotes the cellular uptake of oxygen, whose form is given by tÞÞ if a cell of type i occupies a CA lattice site at position x at time t; where i 2 {H, T} signifying healthy (H) and tumour (T) cells, respectively. Here μ i is defined as the cell type-specific oxygen consumption rate constant which modulates r(c), the oxygendependent consumption rate, which we assume to have Michaelis-Menten form where r c and K m denote the maximal uptake rate and effective Michaelis-Menten constant, respectively. We supplement Eq (1) with the following initial and boundary conditions. We begin with a spatially uniform oxygen distribution c(x, 0) = c 0 , a positive constant, with the entire domain, (0, NΔx) 2 nV, occupied by normal cells. To simulate carcinogenesis, we replace the cell in the central lattice site with a tumour cell. Vessels are placed throughout the domain at a prescribed density and pattern, depending on the situation being simulated. To approximate the constant value of the oxygen concentration in the arterial blood we impose c(x, t) = c max at each point, x 2 V, where there is a vessel. We further impose no-flux boundary conditions at the edge of the domain, where n is the unit normal outward vector to the domain boundary.
Proliferation. Cell proliferation is governed by a complicated set of subcellular processes, and has been modelled in great detail [31,32]. The cell cycle has many levels of complexity including checkpoints, specific temporal sequences and many biophysical processes that alter the cell shape and response to microenvironmental cues and stresses at different times in the cycle [33]. As we are interested in a phenomenon occurring at the much larger, tissue, scale, we will encapsulate this complexity into a simple rule set and threshold values on microenvironmental (oxygen) conditions and neighbourhood occupation, a common simplifying assumption in models of this type [23,34].
At each time step of the cellular automaton (see S1 Text, "Time scales and updates") we visit each cell and check whether the oxygen concentration c at that cell's location exceeds a fixed threshold value c p . If the cell is a healthy cell, then we also check if there is at least one empty lattice site neighbouring the cell, and if a cancer cell, at least one lattice site is not filled by other cancer cells. If either of these conditions are not met, then we consider the cell to be quiescent at this time step. If both of these conditions are met, then we generate a uniform random number r 2 [0, 1] and check if this number is less than a fixed, model specific threshold given by p H or p T if the cell is a 'healthy' or tumour cell, respectively. If it is, then we execute cell division for that cell, placing one daughter cell in the parent cell's location and the other daughter cell in a randomly chosen empty neighbouring lattice site (considering normal cells as empty spaces if the dividing cell is a tumour cell). Under optimal conditions then, the duration of each cell's cycle is described by a geometric random variable.
Cell death. The main mechanisms of cancer cell death considered in this study occur due to severe hypoxia and radiation therapy. Healthy cells can also die if they are selected for replacement by dividing cancer cells, a situation which would normally be ascribed to a cancer cell's greater fitness, allowing it to out-compete healthy cells for space or nutrients locally, or possibly secondary to an acid-mediated event [22,24].
Cells under extreme hypoxic conditions are often found to undergo autophagy (directly translated as 'self-eating'), a state in which they become resistant to nutrient starvation [35], and cells are known to die over different time scales and by different mechanisms (apoptosis versus necrosis) depending on the magnitude and duration of the hypoxic insult. While these differences have been shown to affect tumour growth [36], as this is not the main focus of this model, we will simplify this scenario by assigning a rate, p d , for cell death at each cellular automaton update, when under extreme hypoxia (i.e. c < c ap where c ap 2 (0, c max ] is a fixed, model-specific positive constant).
To model the effect of radiation on heterogeneously oxygenated tissues, we begin with the linear-quadratic model of cell survival: where α and β are radiobiologic parameters associated phenomenologically with cell kill resulting from 'single hit' events (α) and 'double hit' events (β) respectively, and n represents the number of doses of d Gy of radiation. This metric, while originally derived to fit in vitro data [37], has since been widely used to also describe the clinical efficacy of radiation [38,39]. The basic interaction driving DNA damage from radiation is mediated either by free electrons or by free radicals formed by the interaction of photons with water. The DNA damage instantiated by these mechanisms can be repaired by reduction by locally available biomolecules containing sulfhydryl groups [25]. This damage can be made more permanent, however, by the presence of molecular oxygen, which binds with the DNA radical, forming a 'nonrestorable' organic peroxide. This damage 'fixation' by oxygen creates a situation in which the damage now requires enzymatic repair, which occurs on a much longer timescale [25]. Regardless of the mechanism of radiation sensitisation by oxygen, the effect is well documented, and has been modelled by replacing the parameters α and β in Eq (5) with functions of the oxygen concentration using the oxygen enhancement ratio (OER): where α max and β max are the values of α and β under fully physiologically oxygenated conditions and α(c), β(c) are, respectively, the values of α, β at oxygen concentration c. We define the OER as a function of the oxygen concentration by using the empirically established relation [13,40] OERðcÞ which is solved individually for α and β (see Table 1 for parameters).
Quiescence. In response to a lack of oxygen or other nutrients or external mechanical cues indicating overcrowding, cells enter a temporary state of quiescence during which no division occurs [41]. We model this process through an oxygen threshold (c < c p where c p 2 [c ap , c max ] is a fixed, model-specific positive constant) below which cellular division is not possible and by the spatial constraint that the cell is quiescent if there are no empty neighbouring lattice sites (using a Moore neighbourhood). In the case of a cancer cell, we require that at least one space be empty or inhabited by a normal cell, a situation in which we assume that the cancer cell will replace the normal cell if chosen for division. When the conditions for quiescence are no longer present, the cell's state of quiescence is reversed.
Normal tissue approximation. The invasion of cancer cells into healthy tissue has been modelled extensively [34,[42][43][44]. While this process is not the focus of our study, we nevertheless must consider the healthy tissue surrounding our growing tumour, as it plays an important role in modulating the local oxygen concentration. Previous investigations have modelled tumour spheroid growth [45], or assumed a constant influx of oxygen from the boundary of the tissue, but as we seek to understand how heterogeneous vascularisation affects the growth of a tumour, we must also incorporate normal tissue effects [14]. To do this, we assume that the entire tissue is initially comprised of normal cells, which can divide if they have sufficient space, which consume oxygen at a rate of μ H r(c) and are subject to hypoxic cell death. In addition, normal cells can die through competition with cancer cells.
Cell competition. The ability to invade into normal tissue is one of the 'Hallmarks of Cancer' [46]. Several mechanisms for this have been proposed and modelled, including degradation of extracellular matrix by secreted proteases [42,47] and mediated by excreted lactic acid [43,44]. As this specific interaction is not our focus, we make the simplifying assumption that normal cells do not spatially inhibit cancer cell proliferation, and a cancer cell can freely proliferate and place a daughter onto a lattice location where a normal cell is located, causing the normal cell to die, either through an acid-mediated or simple fitness based competition mechanism [48]. Normal cells, however, are spatially inhibited by cancer cells and other normal cells.
Vascular dynamics. While in healthy tissue vessels are normally patent, as tumours grow the vasculature can experience significant spatio-temporal changes on several timescales. There continues to be significant theoretical work to describe angiogenesis, the relatively long timescale process of new vessel formation and pruning [23,28,49], as well as the shorter timescale biophysical processes governing vascular occlusion either through physical forces [50,51] or microemboli [52]. In this study, we seek only to understand the effect of heterogeneous oxygen concentrations on tumour growth and progression [14]. Therefore, we consider only the simplest scenario, in which there is no feedback from the environment to the vessels.

Coupling models in the hybrid cellular automaton
Having described each of the constituent parts of the model, we now describe how they are coupled to form the HCA. We consider a two-dimensional lattice of size N × N, with each lattice element identified by coordinates (iΔx, jΔx) where i, j 2 {1, 2, . . ., N}.
As schematised in Fig 1, at each automaton time step every cell (chosen in random order) in the domain is subject to a series of decisions based on the current automaton state and local oxygen concentration, c(x, t). The specific checks that each cell undergoes, regardless of its phenotype, are: 1. compare c(x, t) to c ap and c p ; 2. consume oxygen at rate μ i r c , die, or remain quiescent; 3. check number of free neighbouring lattice sites; 4. determine proliferative behaviour: if proliferative constraints are met (space and oxygen), then consume additional oxygen for proliferation [53] and place daughter cell in randomly chosen empty neighboring lattice site, otherwise become quiescent.
Cell-microenvironment interaction. As we endeavour to understand the effects of a heterogeneous oxygen distribution on healthy tissues and tumours, we must consider a number of parameters that govern a cell's behaviour in response to these environmental cues. Specifically, we consider the hypoxia threshold, c ap , below which cells will die, and be removed from the system. As cell fate decisions are made on a much slower time scale than oxygen diffusion, there is a potential for update bias which would entail cell fate decisions made long after the conditions for these decisions were met. To reduce this, we will ascribe a probability per automaton update, p d , to this death and removal, and update the HCA more frequently than the cell cycle time. We further consider a baseline oxygen consumption required for cellular processes, r c , and a threshold for cells to be proliferatively active, c p [53].
Non-dimensionalisation and parameter estimation. In order to place the model in an appropriate spatial and temporal scale, we non-dimensionalise the system. We rescale time by a typical cell cycle time, τ = 16 h [54], and lengths by a typical diameter Δx = 50 μm for a glioma cell [55]. We choose to rescale the oxygen concentration to reflect the average oxygen concentration in the arterial blood, estimated at about 80mmHg (5.14 × 10 −13 mol −1 cell −1 ) [56], so that the non-dimensional oxygen concentration at a vessel takes the value 1.
Introducing the non-dimensional variablesx ¼ x=Dx,t ¼ t=t andc ¼ c=c max , we define the new non-dimensional parameters For notational convenience, we henceforth refer to the non-dimensional parameters only and drop the tildes for the aforementioned parameters. See Table 1 for a full list of dimensional parameter estimates and corresponding non-dimensional values. Details of our numerical method for solving Eq (1) are provided in S1 Text. Tumour control probability. The tumour control probability (TCP) is a measure that approximates the probability of eradicating all cells within a tumour and thereby controlling (in this case 'curing') it. While a number of methods, both deterministic and stochastic, have been proposed to define such a measure [62], the most commonly used formulation assumes that the number of surviving cells capable of proliferation after radiation follows Poisson statistics and that the TCP is calculated by where N 0 is the number of cells in the tissue before radiation therapy and SF is the surviving fraction defined by Eq (5).
To compare the effect of radiation on different populations of cells in our model, we calculate the total number of cells surviving a given therapeutic intervention. To this end, we consider the individual survival probability of cells experiencing a given oxygen concentration. Let N k i be the number of cells in the tissue that experience a oxygen concentration in the interval [iΔc, (i + 1)Δc) at time t k = kΔt, for i 2 {0, . . ., M}. Here we define Δc = 0.01, and MΔc = 1, reflecting the fact that the oxygen concentration is bounded by its non-dimensional value in the arterial blood, which is 1. We assume that radiation effect is instantaneous, and represent the time after radiation as t k+1 . The number of cells experiencing an oxygen concentration in the interval [iΔc, (i + 1)Δc) immediately following a single radiation dose d is given by Thus the total number of surviving cells at time t k+1 is given by To calculate the number of surviving cells in a given domain after simulated radiation, we first calculate the individual values of α i and β i for each interval of oxygen concentration using Eqs (6)- (7). We can then compute the surviving fraction after a simulated single 2Gy dose of radiation (d = 2Gy, n = 1) by utilising Eqs (10) and (11) and assuming a distribution of cellular-oxygen given by the cellular automaton at the timestep in question.

Oxygenation of healthy tissue
Regular vascular patterning preserves oxygen distribution shape but not mean. We first model a domain seeded with only healthy cells, and regular spatial distributions of vessels of varying density. Estimates in the literature for 'normal' vascular density span several orders of magnitude and are reported using a number of different measures [26,63,64], so we choose to calculate the actual physiological vascular density of this model by exploring regular vessel spacing and healthy tissue only. In each of these cases, the pattern of the vascular distribution will remain constant, but the spacing between the vessels, and therefore the vascular density (number of vessels / domain size) will change.
As expected, the closer the vascular packing in the domain, the more cells can be supported and the higher the mean oxygen tension (Fig 2). There is a critical value for vessel density below which cells begin to die, and therefore the proportion of the domain inhabited by cells (termed carrying capacity henceforth) drops below unity. This value is between 2.7 × 10 −3 vessels per lattice site and 3.1 × 10 −3 vessels per lattice site for the parameters modelled (Table 1), which is well within an order of magnitude of the vascular area to tumour area ratio reported by Zhang et al. [65].
In addition to the carrying capacity, we investigate the distribution of cellular-oxygen concentration. This distribution provides a non-spatial summary statistic of an otherwise spatial entity: we can understand how many (or what proportion) of cells within the domain exist at certain oxygen concentrations, a measure which will become important later when we consider the effect of radiation. In Fig 2 we note that, while the mean cellular oxygen concentration (mean of the cell-oxygen distribution) varies considerably (from a non-dimensional value of 0.16 to 0.37), the second and third moments of the distribution vary little, with the standard deviation staying between 0.07 − 0.08 and the skewness between 2.23 − 2.45. That is, the distribution of cellular oxygen concentration maintains its variance about its mean as well as its level of symmetry while the vascular patterning is regular.
Irregular vascular patterning results in variation in oxygen distribution shape. We next explore how the oxygen concentration that the population of healthy tissue supported experiences based on changes in density and patterning of vessels in irregularly vascularized domains. While this is a situation that would not likely be seen in the healthy state, understanding the behaviour of the model in this simple situation is valuable before progressing to more complex states. To this end we simulate the resulting healthy tissue growth and maintenance in a variety of randomly generated vascular patterns of the same density (Fig 3). We plot the final, Varying the vascular density with regular spacing affects the carrying capacity and cellular-oxygen distribution in normal tissue. We plot healthy tissue growth and maintenance as we vary the vessel density (decreasing density from top to bottom Θ = 0.0018, 0.0024 and 0.0032). We plot cellular distributions (left) with associated spatial oxygen concentration (middle) and non-spatial distribution of cells versus oxygen concentration (right). These plots represent the system at dynamic equilibrium, in which cell death and birth is balanced across the tissue. The mean of the cellular-oxygen distribution decreases with vessel density (0. 26  stable distribution of healthy cells for a random placement of vessels for a given density (Θ = 0.0018, top; Θ = 0.0024, middle; and Θ = 0.0032, bottom) in the left and middle columns. In the right column, we plot the average distribution of cells at specific oxygen concentrations over ten simulations with different vascular patterns.
As expected, we see an increase in carrying capacity with increasing vascular density, similar to the regularly vascularized scenario, with cellularity (cells / domain size) increasing from 0.557 to 0.945. In these examples of the irregularly vascularised cases, however, we find striking differences in the distribution of cell-specific oxygen concentration (Fig 3, right). While in the regular case (Fig 2) we find conserved cellular-oxygen distribution shapes (second and third moments), in the irregularly vascularized case we find two changes. First, even in the highest vessel density case there is a large peak at c ap and second, the distribution becomes more and more skewed (to the left) as the vessel density decreases (skewness ranging from 0.3820 − 0.9315).

Vascularised tumour growth and invasion
Tumour invasion speed changes with, and can be constrained by, vessel spacing in regular vascular architectures. As a tumour initially invades into healthy tissue, the vasculature that it would encounter would be that of the healthy tissue. While we are ignoring the effects of cell crowding, vascular deformation and angiogenesis, it is of value to understand how the growth rate and the pattern of cancer spread would change given differing vascular spacing in an otherwise regular vasculature. To this end, we seed a series of regularly patterned vascularised domains with increasing regular vessel spacings with a single initiating cancer cell and observe the growth. Note that, to ensure initial tumour growth, we have placed an extra vessel at the centre of the computational domain with the initial cancer cell as in certain regular spacing patterns, the centre would be an area of necrosis and the simulation would end with no cancer cells.
We plot the results of these simulations in Fig 4 along with the individual simulation growth rates. We find four qualitatively different patterns/rates of growth for different vascular densities. First, when there is high vessel density (box 1, green trace, Fig 4), the tumour undergoes rapid growth. This approximates contact inhibited tumour growth with proliferation at the edge only. Even in this growth regime, we begin to see necrotic areas in the bulk of the tumour (box 2, blue trace, Fig 4). The second qualitatively different regime is when the vessel density decreases to Θ = 0.0033, and the spacing of the vasculature is such that it has lost its ability to maintain continuous populations of cancer cells (box 3, red trace, Fig 4). In this regime, we find that the cancer's growth rate begins to slow, and further, that the boundary begins to become irregularly shaped as the cancer cells have to persist in regions of hypoxia long enough to divide. Further, only when the daughters of these cells are placed into the adjacent areas of increased oxygen, will the tumour successfully cross the area of hypoxia, slowing the macroscopic growth. Further reduction in vessel density to Θ = 0.0025 results in areas of necrosis forming in the healthy tissue and a sharp reduction in tumour growth rate (box 4, cyan trace, Fig 4). The final scenario is defined by tissue that is so poorly vascularized that the healthy tissue cannot maintain contiguous populations, and the cancer cannot continue growing beyond the area of initial vessel oxygenation (box 5, black trace, Fig 4), and would need to produce its own vessels to grow further.
Tissue carrying capacity and vessel density. We next investigate in further detail the effect of vessel density on the carrying capacity and cellular-oxygen distribution. We consider a domain of size N = 73 (chosen for ease of comparison for normal vessel spacing) and vary vessel number from very low until we reach saturation of cells (in this case from 3 to 236 vessels). In each simulation, we initialize the domain full of cancer cells and then allow the tissue to experience oxygen dependent birth and death until a dynamic equilibrium is reached, defined as experiencing less than a 1% change in cell number for 50 time steps. We choose this initial condition, instead of seeding the domain with a single cell, because, in this section, we are not interested in growth kinetics, but instead focus on tumour bulk equilibrium characteristics. From this point, we then record and average 100 time steps of cell and oxygen concentration data. For each vessel number, we run 500 simulations, as detailed above, with random vessel configuration and plot the result of each set of simulations and the associated standard deviation (inset) in Fig 5. As expected, we see a monotonic approach to cellular saturation as vessel number increases. When we plot the standard deviation, we find that it is significantly higher in the region that is most physiologically relevant, and that these differences can reach as high as 10% (Fig 5), making cellularity estimates based on vascular density difficult and unreliable. Further, how the Effect on cellular-oxygen distribution in heterogeneous domains. To begin to understand how the patterning of the vasculature affects the cellular-oxygen distribution within our simulations, we consider two cases from the sample of our simulations for two separate vessel densities (24 and 54 vessels). We plot and compare the vascular patterns which support the maximum and minimum cellular populations from two given vascular densities (Fig 6). We plot the minimum cellular population represented in the left column at low (top) and high (bottom) relative vascular densities and the maximum populations in the right column.
In both minimum population cases, there is a large peak of cells near to the hypoxic minimum (c ap ) as well as a smaller peak around an oxygen concentration of 0.5. In the maximum cellularity examples, however, these similarities disappear. In the low vessel density, we see that the hypoxic population is maintained, but we have lost the central, well oxygenated fraction, while in the high vessel density maximum population this is reversed, and we lose the hypoxic population and the entire population is centered around normoxia. This opposite effect (losing the hypoxic fraction in one case and increasing it in the other) would drastically change radiation efficacy, suggesting that vascular organization, not just density, could play an important role in clinical radiation therapy.
Further, we see that in each case, the maximum population is supported by more homogeneous vascular patterns, but that in the two cases the distributions are quite different, and indeed the skewness reflects this. The effect of heterogeneous vascular patterns on surviving cell number As discussed in the previous section, the organisation of vessels in our model can have a significant effect on both the carrying capacity of a domain and also the resulting cellular-oxygen distribution of the cells inhabiting the domain. It is exactly these two values (cell number and surviving fraction) which influence our computation of total surviving cells, a measure critical for comparison across heterogeneous samples, and to understanding radiation efficacy at larger scales, like TCP. How these effects are governed by vessel organisation however, we have not yet elucidated. In this section we will investigate the effects of differing vascular patterns on radiation efficacy. Surviving cell number after radiation varies widely with vessel pattern. To begin to understand the variability inherent in tumour control secondary to radiation effect, we calculate the number of cells surviving after a single simulated dose of 2 Gy of radiation in each simulation and plot the results as a function of vessel density (Fig 7). We previously found that as the vessel density increases, the mean number of cells supported by the domain increases in regular domains, and this holds in irregular domains as well (Fig B in S1 Text). All else being equal, this should translate into an increase in the number of surviving cells after radiation. Indeed, it does translate into a decrease in the surviving fraction after 2Gy, a standard metric used in radiation biology (Fig C in S1 Text). However, we see that the expected relation holds for surviving cell number as well for the low vessel density simulation families, but as the vessel density increases above 0.01, the mean number of surviving cells counter-intuitively begins to decrease. Further, we observe, near this transition point, that the standard deviation of surviving cells within each family of simulations peaks (Fig 7, inset). This highlights the fact that vessel density, and mean oxygen, are insufficient to predict surviving cells, therefore, in order to better predict the number of surviving cells in a domain, we require a measure of vascular organisation.
Vessel density is insufficient to predict radiation effect. We have previously shown that vascular density is a strong predictor for the carrying capacity of the tissue in our model (Fig 5). This is intuitive and, for regularly vascularized domains, as we expect in healthy tissue, is sufficient to explain not only the carrying capacity, but also other measures of the cellular distribution (e.g. the shape of the cellular-oxygen distribution). We have now shown that this measure is not sufficient when we begin to consider heterogeneously vascularized domains, as we know exist in cancer [66], where many possible spatial patterns can exhibit the same density. The mapping from cellular-oxygen distribution to surviving cells is straightforward in the simulations we have presented, as the cellular-oxygen distribution can be measured directly, and the subsequent surviving cells can be calculated using Eq (11). In routinely obtained tissue biopsy specimens however, it is not possible to measuring oxygen concentrations at the cellular resolution. It would therefore be of translational value to define a metric by which cellular-oxygen distribution, or the subsequent radiation efficacy, could be inferred from a readily measurable surrogate, such as the distribution (and density) of microvessels within these samples.

Spatial metrics of vessel organisation correlate with radiation effect
The vascular normalisation hypothesis [67] suggests that radiotherapy should be more efficacious when applied to tissues with normalised vascular beds. To test this in our model, we introduce a spatial metric by which to understand the overall level of spatial heterogeneity in our simulations.
We will utilise the variance stabilized Ripley's L function [21] to measure the deviation from homogeneity in our vessel distributions. This measure, which is a function of distance, describes the average number of points within a given distance of any other point. For a complete description of this measure, see the supplemental information.
Spatial metrics correlate directly with carrying capacity and with radiation response if vessel density is known In order to allow for easy correlation, we first distill Ripley's L function (which is a function of distance) into a single number by taking the mean value from 0 − 19 cell diameters (from adjacent to a distance beyond the ability to affect one another). We calculate the mean Ripley's L function for all of the 500 simulations in each vessel number case and plot it against the associated carrying capacity (Fig C in S1 Text). We find a significant negative correlation for the lower vessel densities which loses statistical significance (p-values not shown) as the domains begin to reach confluence for all but the minority of vessel arrangements. Fig 8 shows scatter plots of Ripley's L function and surviving cells after radiation in 6 representative families of the 500 simulations with increasing vessel numbers. We find that at low vessel densities there is a strong negative correlation between Ripley's L function and surviving cells, but this correlation changes sign at high vessel densities, explaining the counter-intuitive change in surviving fraction after Θ = 0.01 in Fig 7. This suggests that in situations where the vessels are more rarefied, normalising existing vasculature could actually make radiation less effective, in contrast to the vessel normalisation hypothesis.

Discussion
We have used an HCA model of vascular tumour growth in a planar domain to investigate the dependence of cellular oxygenation and radiotherapy efficacy on vascular density and patterning. Our results indicate that simple spatial summary statistics such as Ripley's L function, which could be easily obtained from biopsy images, may predict, together with more standard measures like vessel density, radiation therapy efficacy and the effect of vascular normalisation on this. While our model is parameterised from glioblastoma data, we anticipate these results to be more widely applicable to other cancer types.
Our results corroborate those of Alarcón et al. [14], who also used a HCA approach to studying vascular tumour growth. However, our work differs in several key ways. Specifically, the vasculature in their model lies in plane rather than en face, a choice which was made consciously in our model to ease future translational validation with patient tissue samples. Further, we have ignored several mechanisms of competition in order to focus on defining summary measures of cellular-oxygenation status. These differences aside, Alarcón et al. found, as we did, that heterogeneity on the scale of oxygen concentration affects cell growth, both in overall speed of tumour growth and also in shape of resulting tissue. Further, they reported significant heterogeneity in steady state oxygen concentration across their domains which was dependent on vascular organisation, but this was never explored in terms of radiation effect. A similar finding was reported by Al-Shammari et al. in a biophysical model of healthy muscle tissue [17], this time in a system utilizing en face oxygen sources, as in our work. We have observed similar changes in cellular-oxygen distributions at equilibrium, and indeed, it is these changes that drive our findings concerning radiation effect.
We have seen that the assumption that mean vessel density, and subsequently mean oxygenation, can be used as a surrogate for the number of surviving cells after radiation is insufficient. Further, we found that the relationship between each of our vessel pattern measures and the surviving cells after radiation exhibits a sign change in the mid vessel density range. This change of sign in correlation means that the patterns, within a given vessel density, that take the extreme values of each measure of vascular organization can represent opposite extremes of radiation response. This suggests that if one were to perturb the vasculature, for example toward a more homogeneous distribution, one could induce opposite effects on radiation response, depending on the vessel density.
With Folkman's discovery, in 1971, of a master tumour angiogenesis regulating factor [18], the world thought that the suggested method of blocking this factor, by which it was promised that we could 'starve' tumours of their oxygen supply, would dominate cancer research. It was thought that a cure for cancer would occur in a short time period. However, early trials of single agent anti-angiogenic drugs failed to produce results [20]. Later trials, with combinations of chemotherapy and anti-angiogenic drugs, however, showed promise, but even this was discordant with the leading hypothesis describing the mechanism of anti-angiogenic therapy, which was thought to entirely starve tumours of blood supply.
It was not until 2001, when Jain suggested the 'vascular normalisation hypothesis' [19], that these results could be understood under a single rubric. Jain suggested that anti-angiogenic drugs (sometimes called vascular normalisation therapy, or VNT), instead of entirely blocking new vessel formation, worked to normalise vasculature, pruning inefficient vessels and creating a more regular lattice, thereby improving drug and oxygen delivery. More recent iterations of this hypothesis also include improvement of the efficiency of existing vessels (reviewed by Jain [67,68]).
While this advance in our thinking has provided a way to explain the counter-intuitive results of many trials, it still does not explain why these combinations of anti-angiogenic drugs with chemotherapeutics (or radiation therapy) do not help all patients. Using our simple model system, we have observed a changing correlation with spatial measures and radiation response. This suggests that in certain cases, 'improving' the homogeneity of vascularisation would hurt the radiation effect, whilst in other cases it would help: a heterogeneous response to 'normalisation' of vascular patterning. We show that for certain cases, vessel density held constant, normal vascular patterning can respond either better or worse to radiation therapy (Fig 8).
To translate these conclusions to the clinic would first require biological validation of these hypotheses, to include a careful study of variation in spatial vessel patterning on different spatial scales within the same tumor, or possibly between spatially distinct biopsies. Further, there are a number of model limitations to overcome. Specifically, the assumption that all vessels are en face and all vessels are the same in terms of size and efficiency. Validation could be achieved either through in vivo window chamber experiments, or indirectly through examination of post-radiation surgical specimens. Model development, to include the addition of angiogenesis, biologically derived vessel geometries, vessel size heterogeneity and collapse, could be used to extend these predictions to more realistic situations.
If these limitations were addressed, and validation was achieved, the technology exists currently to take advantage of this new idea. Specifically, macroscopic oxygen concentrations could be inferred from DCE MRI (or other advanced imaging) to create optimised dose plans, as suggested by Malinen et al. [7]. This information about putative vessel density could then be coupled with histologic measurements using automated localisation of vessels [27] and an algorithm to calculate Ripley's L. These two pieces of information could then be incorporated by creating a temporally and a spatially optimised radiation plan whereby the appropriate radiation dose would be delivered before VNT to the area that would suffer after normalisation, and then the final radiation could be delivered after VNT to the area that would benefit from it.
Supporting Information S1 Software. Contains a zipped folder with our java implementation of the mathematical model and bespoke MATLAB analytics as well as a detailed discussion in README.txt describing how to use this software. (ZIP) S1 Text. Contains detailed discussion on our numerical methods of the oxygen transport equation and our statistical methods and results as well as a sensitivity analysis of our mathematical model. (PDF)