{'fraction_non_alphanumeric': 0.04047442226202438, 'fraction_numerical': 0.017517801756148706, 'mean_word_length': 4.804132860040568, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modeling the impact of spatial oxygen heterogeneity on radiolytic oxygen depletion during FLASH radiotherapy

It has been postulated that the delivery of radiotherapy at ultra-high dose rates ("FLASH") reduces normal tissue toxicities by depleting them of oxygen. The fraction of normal tissue and cancer cells surviving radiotherapy depends on dose and oxygen levels in an exponential manner and even a very small fraction of tissue at low oxygen levels can determine radiotherapy response. The effect of FLASH on radiation-induced normal and tumour tissue cell killing was studied by simulating oxygen diffusion, metabolism, and radiolytic oxygen depletion over domains with simulated capillary architectures. Two architectural models were used: 1.) randomly distributed capillaries and 2.) capillaries forming a regular square lattice array. The resulting oxygen partial pressure distribution histograms were used to simulate normal and tumour tissue cell survival using the linear quadratic model of cell survival, modified to incorporate oxygen-enhancement ratio (OER) effects. Tumour cell survival was found to be increased by FLASH as compared to conventional radiotherapy, with a 0-1 order of magnitude increase for expected levels of tumour hypoxia, depending on the relative magnitudes of radiolytic oxygen depletion and tissue oxygen metabolism. Interestingly, for the random capillary model, the impact of FLASH on well-oxygenated (normal) tissues was found to be much greater, with an estimated increase in cell survival by up to 10 orders of magnitude, even though reductions in mean tissue partial pressure were modest, less than 7 mmHg for the parameter values studied. The presence of very small nearly hypoxic regions in otherwise well-perfused normal tissues with high mean oxygen levels resulted in a greater proportional sparing of normal tissue than tumour cells during FLASH irradiation, possibly explaining empirical normal tissue sparing and iso-tumour control results.


Introduction
Within the oxygen depletion hypothesis, it is commonly argued that FLASH radiotherapy reduces normal tissue toxicity by decreasing the level of oxygen in these tissues, thus rendering them more radioresistant due to the OER effect [1][2][3][4][5]. Low linear-energy-transfer (LET) radiation damages DNA primarily through the generation of free radicals. Oxygen is needed to "fix" this damage, rendering it irreparable; the lower the oxygen levels in tissue, the lower the impact of radiation [6,7]. This effect is usually quantified by the OER, defined as the ratio Dhypoxic/D(p), where Dhypoxic is the dose needed to achieve a given level of cell killing under maximally hypoxic conditions and D(p) is the dose needed to achieve the same level of cell killing at an oxygen partial pressure p. Viewed as a function of partial pressure, the OER is equal to one at zero tissue partial pressure, p = 0. As the oxygen partial pressure is increased, OER(p) initially increases very rapidly, nearly reaching its maximum value (between 2.5 and 3) by 10-15 mmHg [8,9], but changes slowly above this range. The mean oxygen partial pressures ̅ of most normal tissues range from ~20 to 50 mmHg [10]. Based on the slow evolution of OER(p) with respect to p at partial pressures above ~ 10 mmHg, one might thus expect that a substantial (≳ 10-40 mmHg) reduction in partial pressure by radiation would be needed to significantly alter cell killing in normal tissues and hence, reduce normal tissue toxicities. Recently, this has been used to argue that the oxygen depleting effects of FLASH cannot be responsible for its apparent reductions in toxicity, due to comparatively small reductions in oxygen levels [11,12].
From both in vitro [13][14][15] and in vivo [16] experiments, radiation has been shown to deplete oxygen in a dose-dependent manner by 0.1-0.42 mmHg per Gy of delivered radiation, defining the rate of radiolytic oxygen depletion (ROD). Here, ∆ ̅ is the change in the mean oxygen partial pressure in the system being studied and D is the administered dose. ROD occurs whenever the dose delivery time is smaller than the characteristic time scale needed to restore oxygen levels depleted by radiochemical reactions [13]; in vivo, this is the time (~1 s) needed for oxygen to diffuse out of capillaries into surrounding tissue and dose must be delivered in ≲ 1 s to observe these effects [17]. In preclinical in vivo experiments, the normal tissue toxicity-reducing effects of FLASH have been observed after single-fraction doses between 10 and 50 Gy [18][19][20][21]. Based on measured values quoted above, the lower end of this range is expected to give rise to a reduction in mean tissue oxygen partial pressure by 1-4 mmHg, consistent with measured values by Cao and colleagues, who reported FLASH-induced reductions of ~1 and 2 mmHg in tumours and normal tissues, respectively, after a single 20 Gy dose [16]. These values are much smaller than the magnitude estimated previously (≳ 10-40 mmHg) that would naïvely be needed to significantly alter the OER in normal tissue. This led Pratx and colleagues to propose that normal tissue hypoxic stem cell "niches" (≲ 10 mmHg) are responsible for the radiation response of normal tissues [2,17], since such cells would be expected to be in the range of oxygen levels where changes induced during FLASH would significantly influence cell survival.
This argument hints at an important point: even in nominally well-perfused normal tissues, the vascular network is irregular [22][23][24][25][26], and while more efficient at delivering oxygen than the typical tumour vasculature [27], oxygen distributions still exhibit a substantial degree of heterogeneity [25,[28][29][30]. Consequently, in well-oxygenated normal tissues (mean partial pressures ̅ ≳ 20 mmHg), there may be hypoxic regions where the local partial pressure p ≲ 10 mmHg and hence, FLASH susceptible. Even if these are proportionately small compared to the normoxic regions, radiation-induced cell survival varies exponentially with the effective dose •OER, and the response of tissue to radiotherapy will be determined by these most hypoxic cells.
The heterogeneity of oxygen partial pressures within tissues is thus a potentially important factor in determining the impact of FLASH. In this work, we simulated its effect by solving the oxygen reaction-diffusion equation [31,32] including the effect of ROD [2,4] over domains containing capillaries (oxygen sources) distributed randomly and on the vertices of a square lattice, thereby simulating two extremes of geometric irregularity (and hence, oxygen heterogeneity). The resulting oxygenation maps were used to calculate target cell survival in normal and tumour tissues using the linear quadratic model, modified to account for oxygen sensitizing effects. These simulations were run with and without ROD, allowing us to calculate the differential cell survival between FLASH and conventional dose rate radiotherapy (herein referred to simply as "conventional radiotherapy"). Our work thus builds on the calculations of Pratx et al. [2] and Rothwell et al. [4], who carried out related simulations for a single cylindrical capillary ("Krogh geometry").
To connect these cell survival results with empirical preclinical studies showing reduced normal tissue toxicities and apparent iso-tumour control [18,33,34] using FLASH, we calculated the ratio of the FLASH and conventional radiotherapy doses needed to achieve isocell survival for normal and tumour tissues. These were compared with empirical dosemodifying factors for iso-toxicity, which have been measured in several tissues, including lung, brain [19,33], and skin [35]. The limitations of equating cell survival with toxicity are considered in the Discussion.

Cell survival analysis
The differential impact of FLASH and conventional radiotherapy was quantified by estimated normal and tumour tissue cell survival fractions (SF's). SF's were calculated using the linear-quadratic (LQ) model for a single fraction of radiotherapy (following most preclinical experimental FLASH regimens), modified to account for OER effects by convolving the LQ expression for survival fraction with the time-dependent distribution ( , ) of oxygen partial pressure values within the tissue: Here, is the biologically effective dose corresponding to a dose , modified by the oxygen partial pressure-dependent OER. Equation (2) generalizes a well-known expression for the oxygendependent linear quadratic formula [8,32] to the case where the oxygen distribution ( , ) is time-dependent, as happens during FLASH radiotherapy, here taken to be delivered between a time = 0 and = .
The partial pressure dependence of the oxygen enhancement ratio was approximated as [8,36,37] OER(p) = • OER max + + where OER max = 3 is the maximum OER value, achieved when ≫ . The partial pressure at which OER achieves half its maximum value is estimated to lie between 1.9 and 3.28 mmHg [8,37]. α and β in Equations (2) and (3) quantify the linear and quadratic dependencies of radiosensitivity with respect to dose; values for these parameters and all others used in the calculation of SF are shown in Table 1.
It was assumed that Equation (2) applied equally well to normal and tumour tissues, with different values of α and β. For tumours, the ratio α β ⁄ was chosen to be 10 Gy, while for normal tissues, it was taken to be 3 Gy, standard values for early-and late-responding tissues [38]. Because the values of α under oxic conditions are highly variable between tumour [39][40][41] and normal [41] tissues, we simulated normal and tumour cell survival using a range of values: 0.2, 0.3, and 0.4 Gy -1 , consistent with empirical in vitro values. Equation (2) also assumed that the density of "target" cells (e.g., tumour clonogens) is uniform with respect to oxygenation.
The heterogeneity of oxygen partial pressures within a tissue was quantified by the distribution ( , ), which represents the probability that a given point in space will have an oxygen partial pressure , a time after the start of FLASH. As noted in the Introduction, even though the mean oxygen partial pressure ̅ = ∫ ∞ 0 ( , 0) before FLASH is delivered is sufficiently large in normal tissues that OER( ̅ ) ~ OER max , any irregularities in the spatial distributions of capillaries means that there may be regions at low partial pressure, (~0, 0) ≠ 0. These radioresistant regions, for which OER ≈ 1, will ultimately determine SF owing to its exponential dependence on • OER( ) in Equation (2). In what follows, we will suppress the time variable when discussing the pre-FLASH equilibrium oxygen partial pressure histogram: ( ) ≡ ( , ≤ 0).

Oxygen partial pressure distributions
Oxygen partial pressure distributions ( , )'s were calculated by solving the oxygen reactiondiffusion equation for the oxygen partial pressure ( ⃗, ) at position ⃗ on a two-dimensional domain, subject to ( ⃗, ) = for ⃗ on the surface of capillaries (radius ), with the capillary oxygen partial pressure. Details of the two capillary architecture models are described below.
Equation (5) combines the effects of oxygen diffusion (diffusivity O 2 ) with metabolism (the second term on the right-hand side) and radiolytic oxygen depletion (third term on the righthand side). Oxygen metabolism is the same as the Michaelis-Menten form used by Daşu and collaborators [42] and Petit and colleagues [32], with max the maximum metabolic rate and = 2.5 mmHg. The radiolytic oxygen depletion (ROD) term was assumed to similarly obey Michaelis-Menten kinetics; it reduces to the form used by Pratx and collaborators [2] in the limit that ROD → 0. In vitro studies have discerned no deviation from a constant ROD as a function of oxygen partial pressure down to several mmHg [13][14][15], and yet ROD must vanish as the partial pressure does; correspondingly, we chose a small value for the partial pressure ROD = 1 mmHg at which the ROD is half its maximum value. We solved Equation (5) with and without the ROD term to simulate the effect of FLASH on oxygen distributions.
Following the approach of Weiss and others [2,13,43], Equation (5) assumes that ROD is simply proportional to the dose rate ⁄ , where is the dose delivered in a time , and 0 is the constant of proportionality. More generally, the kinetics of the radiochemical reactions responsible for ROD and oxygen fixation of radiation-induced DNA damage are properly described by second-order reaction equations [4,44]. As shown by Ling, however, the results of calculations using these second-order reaction equations approach those using the simpler kinetics in Equation (5) for radiation pulses lasting longer than the time scales describing oxygen fixation ( fix ) and ROD kinetics ( ROD ), both on the order of several microseconds (μs) [44,45]. (At the same time, as noted in the Introduction, the pulse duration must be smaller than the time needed for oxygen to diffuse away from capillaries for there to be substantial ROD.) For pulsed irradiation, we further require that the time between pulses is much longer than fix , ROD . Most recent experiments, involving either pulsed or continuous irradiation well-satisfy all these conditions (see Table 1 in the review paper by Wilson and colleagues [34]).
It was further assumed that the capillary oxygen partial pressure was unaffected by FLASH, based on the fact that the hemoglobin oxygen dissociation rate ~20 s -1 [46] is much faster than typical values of 1⁄ . We thus expect that ROD in blood is immediately compensated for by hemoglobin giving up oxygen to maintain a near-constant partial pressure.
Random two-dimensional capillary domains have been used previously to simulate the oftendysfunctional tumour vasculature [32,42]. In contrast, normal tissues are often assumed to have a well-organized vasculature, optimized to meet metabolic demands, and modelled as a regular array of non-interacting cylinders (Krogh geometry) [22], including in FLASH studies [2,4]. There is a growing appreciation, however, that the normal tissue vasculature exhibits a substantial degree of irregularity, and that this irregularity may be important for understanding tissue physiology [23,25]. The hallmarks of such irregularity are regions of hypoxia and, related to this, a long tail in the distribution n(DNC) of the minimum distances between sampled points in the tissue extravascular space and the nearest capillary, i.e., the "distance-to-nearest-capillary" (DNC). This distribution as well as the closely-related minimum intercapillary distance distribution have been studied extensively in both normal [24-26, 28, 47-50] and tumour [24,48] tissues and both distributions indeed exhibit such a tail at long distances (> 50 m), indicating the likely presence of small regions of poorly oxygenated regions, even in healthy normal tissues [25,28]. This has been confirmed in measurements of normal tissue oxygen partial pressure distributions, which show the presence of anoxic regions [ (~0) ≠ 0] even in nominally well-oxygenated tissues such as liver, brain, and spinal muscle [30].
We calculated oxygen partial pressure distributions ( , ) as well as DNC histograms for both randomly distributed capillaries and capillaries distributed on the vertices of a square lattice to represent the two extremes of vascular irregularity used in the literature. Based on the similarity between ( ) and DNC histograms calculated using the random capillary model and empirical tumour and normal tissue results, we argue in the Discussion that the random capillary model is a better representation of the vasculature in both these tissue types, although results are presented for both. We emphasize that we are not asserting that the normal tissue vasculature is random, only that it exhibits sufficient irregularity that the predicted oxygen partial pressure distributions and architectural metrics arising from a random capillary model better match with empirical results than e.g., a model representing a regular array of capillaries.
Parameter values used in our simulations are shown in Table 1. Of note, simulations were carried out using single-fraction doses of 20 Gy delivered in T = 0.4 s, two representative values of ROD, G0 = 0.2 and 0.4 mmHg/Gy, a capillary oxygen partial pressure = 40 mmHg [4], and maximum oxygen metabolic rates max of 5, 15, and 40 mmHg/s for both normal and tumour tissues. Justification for these choices is given below.
Recent in vivo preclinical experiments have investigated the effects of FLASH irradiation using single-fraction doses spanning 10 to 30 Gy [16,20,21,51]; we used 20 Gy as a representative value, also the value used in the experiments of Cao et al. [16], the results of which are most comparable to our own. Assuming = 0.4 s gives a dose rate of 50 Gy/s, which is a typical dose rate at which normal tissue sparing effects have been observed [34].
We use the notation G0 to distinguish the rate of ROD that enters the oxygen reactiondiffusion equation, Equation (5), from the empirically observed values G, defined in Equation (1). As noted by Pratx and Kapp, the two values only coincide when oxygen diffusion and metabolism are ignored [2]. We use G0 = 0.2 and 0.4 mmHg/Gy, approximately encompassing the range of measured values in vitro [13][14][15], where oxygen distributions are fairly uniform spatially and hence, diffusion is less relevant.
The choice of = 40 mmHg corresponds to the expected venous oxygen tension [10,52]. Although partial pressures drop across capillaries from ~ 90-100 mmHg in arteries to 35-45 mmHg in veins [10,52], radiation response will depend on the least well-oxygenated components of the tissue, i.e., near the venule ends of capillaries. In the Discussion, we discuss the impact of this assumption on our results in some detail, arguing that alternative choices of (or any other of the values used in our simulations) do not change the conceptual results of our study. This also includes the maximum rates of oxygen metabolism max , assumed here to be 5, 15, and 40 mmHg/s. The value 15 mmHg/s is a classic value for tumours [32,42,53,54], based on measurements of human tumours by Warburg [53]. Individual tumours exhibit a range of metabolisms, however [55]. Similarly, normal tissues exhibit a range of metabolic activities, with liver and kidney generally exhibiting a higher metabolic rate than the average tumour, brain and intestines exhibiting comparable levels, and resting skeletal muscle exhibiting lower levels (see Ref. [30] as well as Table SM1 in the Supplementary Materials). Further discussion of tumour and normal tissue metabolic levels is given in the Supplementary Materials Section 1, where we argue that 5-40 mmHg/s encompasses an appropriate range for most normal tissues.
For a specified value of the areal capillary density (= 2 ⁄ , where is the number of capillaries and = 1 mm is the size of the simulation domain), Equation (5) was solved using a finite-element technique, having randomly distributed the capillaries over the domain or having placed them on the vertices of a regular square lattice with spacing ≈ √ ⁄ . Zeroflux Neumann boundary conditions were imposed on the partial pressure ( , ⃗⃗⃗ ) along the edges of the domain boundaries. For each capillary configuration, Equation (5) was solved for conventional RT (approximated as ( ) ⁄ • 0 = 0; from here on, we simply indicate conventional RT as 0 = 0) and for FLASH RT ( 0 = 0.2 and 0.4 mmHg/Gy). The equilibrium (static) conventional RT solution was used as the initial condition ( , ⃗⃗⃗ 0) for the corresponding FLASH calculation, carried out for 0 ≤ ≤ . Figure 1 shows the results of one such calculation for a random capillary array. Following the approach of Ref. [32], for the random capillary arrays, the calculation was repeated at least 100 times for different random capillary placements but fixed , to ensure adequate statistics in the calculation of SF [Equation (2)]. For calculations involving the lattice capillary array, only a single calculation was needed for each . At discrete timepoints t = 0, ∆ , 2∆ , . . . , separated by the increment ∆ , the resulting spatial oxygen partial pressure maps ( ⃗, ) were sampled over the 1mm × 1mm spatial domains, excluding the outermost regions < 100 m from the boundaries, to generate oxygen partial pressure distributions ( , ) for each value of for conventional and FLASH radiotherapy. This exclusion was done to avoid finite-size effects arising from the zero-flux boundary conditions. For each simulation and each timepoint, 900 × 900 points were thus sampled over the domains, corresponding to one partial pressure data point per m 2 . These datasets were combined to generate a total dataset of partial pressure values for each and timepoint of size ≥ 8.1 × 10 7 (= 900 × 900 × 100). Values of spanned 8 mm -2 to 637 mm -2 . Examples distributions are shown in Figure 2 for three representative values of for t = 0 and t = . The sampled partial pressure datasets were used to calculate SF, approximating the integral as a sum [Equation (6)] over the number of sampled points and a trapezoidal time integration, where ( ) is the ith element of the sampled histogram dataset at timepoint . ≡ ∆ ⁄ is the number of timepoints used in the integration. Because the spatial datasets were large, we used only = 5 timepoints. At each point , ⃗⃗⃗ ( ⃗, ) decays smoothly with time and the numerical error associated with this relatively coarse integration step size is small. Survival fraction curves were computed as functions of the mean tissue partial pressure without ROD (equivalently, pre-FLASH), similarly discretized as shown in Equation (7): Here, (0) is the ith element of the sampled histogram dataset without ROD (i.e., 0 = 0). We also calculated the change in mean partial pressure over the course of FLASH; i.e., the difference in partial pressure between the end (t = T) and beginning (t = 0) of FLASH.

Dose modifying factor
To connect the SF results with empirical normal tissue toxicity and tumour control results, we defined the dose modifying factor (DMF) as where conv, iso-SF is the single-fraction dose of conventional radiotherapy needed to achieve the same survival fraction ("iso-SF") as 20 Gy of FLASH. Having calculated the SF for 20 Gy and the FLASH oxygen distribution for a given capillary density using Equation (6), Equation was then solved iteratively using the corresponding conventional dose-rate oxygen distribution to find the dose conv, iso-SF that resulted in the same SF. The DMF was calculated for both tumour (α β ⁄ = 10 Gy) and normal tissues (α β ⁄ = 3 Gy). Results are shown for the representative metabolic rate max = 15 mmHg/s and both simulated rates of ROD.

FLASH-induced reduction in oxygen partial pressures
FLASH lowered mean tissue oxygen partial pressures by an amount dependent on the initial mean oxygen partial pressure, the value of G0, and the metabolic rate cmax; see Figure 3. Tissues with low mean partial pressure ̅ values were largely refractory to FLASH, with ∆ ̅ approaching zero as ̅ did. With increasing partial pressure, ∆ ̅ grew, reaching a maximum when ̅ ≳ 20 mmHg. For the random capillary array model and intermediate metabolic rate max = 15 mmHg/s, maximum ∆ ̅ values were ~ 3 mmHg for G0 = 0.2 mmHg/Gy and ~ 7 mmHg for G0 = 0.4 mmHg/Gy, respectively. For the simulated dose of 20 Gy used in our calculations, these correspond to the maximum ROD rates G slightly below G0. As can be seen in Figure 3, ∆ ̅ decreased with increasing metabolic rate, since the magnitude of the ROD effect depends on the relative sizes of ( ⁄ ) • 0 and max ; see Equation (5). Results for square lattice capillary configurations were quantitatively similar, and only results for the random capillary model are shown.

Impact of FLASH on tumour and normal tissue cell survival
For both models of capillary architecture studied, FLASH led to a pronounced increase in cell survival fraction (SF), depending on the mean tissue oxygen partial pressure ̅ values. Example plots of SF's versus ̅ for α = 0.3 Gy -1 and α β ⁄ = 3 Gy are shown in Figure 4 for the three metabolic rates studied (cmax = 5, 15, 40 mmHg/s), a dose of 20 Gy, and the random (top) and lattice (bottom) capillary architectures. SF curves for the other simulated values of α and β were qualitatively similar and are shown in the Supplementary Materials Section 2.
At low ̅ values (hypoxic tissue), tissues were again found to be refractory to the impact of FLASH, and the SF curves for conventional and FLASH radiation coincided in the limit ̅ → 0. For tumours (α = 0.3 Gy -1 and α β ⁄ = 10 Gy), the difference in SF's between FLASH and conventional radiotherapy remained ≲ 1 order-of-magnitude for ̅ ≲ 10 mmHg (the expected range for most tumours), using cmax = 15 mmHg/s; see Figure SM1 in the Supplementary Materials.
With increasing ̅ , the SF curves separated; for the lattice capillary geometry, a maximum separation was achieved for ̅ ~5-20 mmHg (depending on parameter values), while for the random capillary model, the curves separated monotonically with increasing ̅ . For the random capillary model, the separation was substantial, with normal tissue cell survival increasing by as much as 10 orders of magnitude, depending on the relative magnitudes of G0 and cmax ( Figure  4).

Dose modifying factors
This differential cell-killing manifested itself in tumour and normal tissue dose-modifying factors ≥ 1, rising from 1 at ̅ = 0 to maximum values approaching 1.7 (for cmax = 15 mmHg/s and G0 = 0.4 mmHg/Gy); see Figure 5. Similar to the survival fraction results, although the magnitude of the maximum DMF was comparable for the random and lattice capillary models (not shown), the maximum DMF for the lattice model was achieved at relatively low mean partial pressure values, ̅ ~ 5 −10 mmHg, before decreasing rapidly. In contrast, for the random capillary model, DMF continued to grow with increasing ̅ , achieving a maximum for ̅ ≳ 20 mmHg, the expected range of oxygen tensions for normal tissues. To the extent that tumours have lower ̅ than normal tissues, the impact of FLASH was thus greater for normal tissues than for tumours. Figure 5 shows representative DMF curves for tumour (α β ⁄ = 10 Gy) and normal (α β ⁄ = 3 Gy) tissues using α = 0.3 Gy -1 , cmax = 15 mmHg/s, and different values of the rate G0 of ROD. cmax = 15 mmHg/s was assumed here to represent both tumour oxygen metabolism as well as normal tissue at an intermediate metabolic rate; see Supplementary Materials Table SM1.

Discussion
By depleting tissues of oxygen, FLASH radiotherapy was found to increase cell survival because of the oxygen enhancement ratio (OER) effect. The difference in cell survival between FLASH and conventional radiotherapy grew by orders of magnitudes with increasing tissue oxygenation, providing a possible explanation of the observed reductions in normal tissue toxicity and apparent iso-tumour control [18,33,34], since tumours are often found at lower oxygen levels than normal tissues, with mean oxygen tensions typically below 10 mmHg [56][57][58][59][60][61]. For completely anoxic tissue (i.e., oxygen partial pressure exactly zero), FLASH cannot further deplete oxygen from tissue and the survival curves for FLASH and conventional radiotherapy coincide.
The substantial impact of FLASH on cell survival in well-oxygenated tissue is the central result of this manuscript, contradicting the naïve expectations that oxygen sensitization effects based on mean tissue partial pressures should be small [2], and that oxygen must be depleted completely for the ROD mechanism to contribute significantly to toxicity reductions [11,12]. The relative increase in radioresistance with increasing oxygenation is a result of the heterogeneity of the oxygen distribution in tissues. A single dose of 20 Gy kills all welloxygenated cells (SF~10 -20 for OER = OER max using e.g., α = 0.3 Gy -1 and α β ⁄ = 3 Gy) and hence, the survival fraction is essentially the fraction (~0) of maximally radioresistant cells at ultra-low oxygen levels. FLASH increases cell survival by increasing the number of cells in this population. Because (~0) is small for well-oxygenated tissues (see e.g., the middle and right panels of Figure 2), the increase in (~0, > 0) during FLASH is proportionately greater for such tissues and the relative cell survival-FLASH versus conventional radiotherapy-likewise grows with increasing mean partial pressure. The exponential nature of cell survival after radiation means that even a very small fraction of hypoxic cells (e.g., less than 0.5% of the non-FLASH irradiated tissue represented in the right panel of Figure 2 is estimated to have partial pressures less than 5 mmHg) will have a large impact.
To different extents, this behaviour arose in both the capillary geometry models that we used, the random capillary model as well as the lattice model. However, the effect was far more robust using the random capillary model. For the lattice model, the window of mean tissue partial pressures where FLASH was predicted to result in a substantial sparing effect was substantially reduced compared to the random capillary geometry since the lattice geometry eliminates the possibility of any tissue being overly distant from capillaries. On geometric grounds, for the lattice geometry, (~0, ) will be nonzero only when the maximum distanceto-nearest capillary (i.e., the point in the middle of the lattice unit cell) is greater than the critical distance (̃m ax ) ≡ √4 O 2̃m ax ⁄ past a single capillary beyond which the oxygen tension vanishes [53]. Equivalently, when the capillary density satisfies Here, ̃m ax = max + 0 ( ⁄ ) is the approximate effective metabolic rate, accounting for the effects of FLASH [In Equation (5), FLASH is formally equivalent to a change in oxygen metabolism in the limit where the Michaelis-Menten parameters are identical, ROD = ]. For larger capillary densities (equivalently, larger mean tissue partial pressures, ̅ ), (~0, ) will be zero for all ≤ , and the effect of FLASH on cell survival vanishes: lattice geometry: In contrast, for the random capillary model, (~0, ) decays exponentially at large mean tissue partial pressure, but remains nonzero, yielding an exponential separation of the survival curves with increasing tissue oxygenation ( ): random capillary geometry: Here, is a constant of order unity; see the Supplementary Materials Section 3 for a derivation of this result. These respective behaviours-the absence of anoxic regions in the lattice model above a critical capillary density and the convergence of the survival curves for the lattice model and exponential separation for the random capillary model-are apparent in Figures 2  and 4.
The question of which model-random capillary or lattice array-better describes the actual normal and tumour vascular architectures is thus an important one, insofar as only the random capillary model predicts a robust increase in the impact of FLASH with increasing tissue oxygenation: as evident from Equation (12), the choice of parameter values , O 2 , c , max , and 0 , used in our simulations will affect the rate at which the survival curves separate with increasing oxygenation, but not the existence of such a therapeutic advantage. In contrast, for a regular lattice capillary array, unless Equation (10) is satisfied, the advantage will not exist. Note that to satisfy the metabolic demands of the majority of a given normal tissue, one would generally expect that ≳ 1 2 2 ⁄ and hence, Equation (10) may not be satisfied for normal tissues, with (~0) vanishing for physiologically relevant parameter values. As noted above, it is the fact that (~0) is nonzero-albeit it, very small-in the random capillary model even at high mean tissue oxygen tensions, ̅ ≳ 20, that produces the large separation of survival curves, with FLASH producing a significant tissue sparing effect. Empirical studies have indeed found small but nonzero (~0) even in well-perfused normal tissues such as liver and brain where ̅ ≳ 20 mmHg [30], in sharp contrast to the ( )'s calculated using the lattice capillary model (see e.g., the lower panel in Figure 2), for which (~0, ) is identically zero for large capillary densities, consistent with Equation (10).
This behaviour supports the conclusion that, in terms of the predicted oxygen distributions at least, the random capillary model may be a better model of the normal and tumour tissue vasculature than a regular array model. Further evidence for this can be found by comparing calculated distance-to-nearest capillary (DNC) histograms. As noted in the Introduction, vessel morphometry studies have quantified this distribution and the closely-related minimum intercapillary distance distribution in normal [24-26, 28, 47-50] and tumour [24,48] tissues. A key hallmark of vascular irregularity is a "tail" in these distributions at long distances (≳ 50 m), indicating the presence of small regions of poor vascularization. In Supplementary Materials Section 4, we show the results of a calculation of n(DNC) for a well-perfused simulated normal tissue with a random vasculature and compare it with results for tissue with a lattice of capillaries. The former is found to result in an n(DNC) that qualitatively matches empirical normal tissue results, while the latter fails to reproduce the characteristic tail. Hence, we conclude that the random approximation is a likely an adequate model of normal and tumour tissue vasculature to predict the impact of FLASH on both these tissues. We emphasize, however, that differences in vascular architectures amongst normal tissues (e.g., brain versus skin) will result in different normal tissue sensitivities to radiation, and a better understanding of vascular irregularity as well as oxygen heterogeneity in different tissues will help improve our assessment of the magnitude of FLASH's impact.
The fact that the dose-modifying effects of FLASH are only felt by cells at low local oxygen partial pressures led Pratx and colleagues to hypothesize that normal tissue stem cells in hypoxic "niches" are primarily responsible for the radiation response of normal tissues to radiation [17]. Our tissue oxygenation distribution simulations suggest that such hypoxic niches may be present in all tissues, even nominally well-vascularized ones. Whether or not the target cells (normal tissue stem cells in the work of Pratx et al.) preferentially populate lowoxygen regions or are uniformly distributed with respect to oxygenation (as assumed in this work), the high oxygen enhancement ratio of cells in these niches will protect them from radiation, characterized by the large enhancement in cell survival fraction.
The magnitude of the impact of FLASH-both on the ROD-induced change ∆ ̅ in mean oxygen partial pressure and the change in cell survival fractions-depended on our chosen parameter values. Although a comprehensive assessment of the range of parameter values for which FLASH can have a substantial effect has been undertaken by Rothwell and colleagues, who also used a more detailed model of FLASH kinetics [4], we make two observations here. Using dimensional analysis, the impact of FLASH depends on the magnitudes of the two dimensionless parameters As noted previously, is expected generally to be on the order unity for normal tissues, ensuring that most of the tissue is at nonzero oxygen partial pressures (note that this does not preclude the possibility that (~0) ≠ 0, only that it is a small number, much less than one). For all our choices of parameter values corresponding to expected normal tissue mean partial pressure values, ̅ ≳ 20 mmHg, varied from ~ 0.5 to ~ 3, indicating that our parameter choices were physiologically reasonable. Hence, irrespective of our precise choices of parameter values-, O 2 , c , and max -the large separation between survival fraction curves shown in Figure 4 is expected to hold for normal tissues, defined here as tissue for which ̅ ≳ 20 mmHg, but more generally as tissue for which ≳ 1. This leaves as the only "free" parameter, unconstrained by physiology. It represents the ratio of the ROD-induced temporal rate of oxygen depletion and the rate of oxygen metabolism. If this quantity is small, then ROD represents a small perturbation on top of the baseline oxygen metabolism; conversely, if it is large, ROD substantially renormalizes the effective oxygen metabolism. Our smallest simulated value-corresponding to max = 40 mmHg/s and 0 = 0.2 mmHg/Gywas 0.025. Use of this value gave rise to a maximum ∆ ̅ value of ~ 2 mmHg (see Figure 3), yielding maximum values ≡ ∆ ̅ ⁄~ 0.1 mmHg/Gy of the empirical rate of ROD, smaller than 0 . Our largest value of -corresponding to max = 5 mmHg/s and 0 = 0.4 mmHg/Gywas 4. Use of this value gave rise to a maximum ∆ ̅ value of ~ 6.5 mmHg, corresponding to = 0.33 mmHg/Gy, again slightly below 0 . In comparison, empirical rates of ROD span 0.1-0.42 mmHg/Gy [13][14][15][16]. Hence, our choice of values for max and 0 largely encompassed observed rates of ROD.
The fact that ∆ ̅ was predicted to vary with respect to ̅ (Figure 3) is a consequence of oxygen heterogeneity and is an important differentiating feature between in vitro and in vivo experiments, since oxygen should be much more homogeneously distributed in the former. In several in vitro experiments, ∆ ̅ was found to be nearly independent of the pre-FLASH mean oxygen partial pressure ̅ down to several mmHg [13][14][15]. In contrast, in a recent in vivo experiment, Cao et al. found ∆ ̅ ~ 1 mmHg for tumours (pre-FLASH ̅ ~ 10 − 14 mmHg) but ∆ ̅~ 2.3 mmHg for normal tissues (pre-FLASH ̅~ 20 − 35 mmHg) [16], for 20 Gy. This is close to the behavior shown in Figure 3 for e.g., G0 = 0.2 mmHg/Gy and max = 15 mmHg/s. As with our survival fraction results, this behaviour can be understood as arising from the persistence of a fraction (~0, ≤ ) ≠ 0 of the tissue that is refractory to FLASH, resulting in the slow increase in ∆ ̅ with respect to ̅ in our simulations, likely mimicking the in vivo situation. In contrast, the relative homogeneity of oxygen in vitro means that there is no such fraction for mean partial pressures above several mmHg's and the decrease in mean partial pressure is independent of ̅ .
Because of the up to 10 order-of-magnitude increase in cell survival fractions for welloxygenated tissues ( ̅ ≳ 20 mmHg; Figure 4) using the random capillary model, for the range of parameters studied by us, FLASH doses would need to be increased by 10%-70% to achieve the same survival fraction as conventional radiotherapy ( Figure 5). This is comparable to the range of empirical iso-toxicity dose modifying factors, ~1.2-1.8, observed pre-clinically [18,19,33,35,62,63] (see Ref. [34] for a review of DMF's). Note that some of the higher (>1.4) empirical dose-modifying factors established iso-toxicity using larger doses of FLASH radiation, 25 -75 Gy [18,34,35,63], which are expected to result in greater radiolytic oxygen depletion than the 20 Gy used in this work.
The effective dose-modifying impact of FLASH radiotherapy on normal tissues should be contrasted with the apparent iso-efficacy of FLASH and conventional radiotherapy in terms of preclinical tumour volume kinetics [18,33,34]. Although our results predict that FLASH should have a substantially smaller impact on poorly oxygenated tissues such as tumours, the apparent absence of any effect in preclinical experiments is striking. It is likely that the effect of FLASH on tumours is multi-factorial, producing an enhancement in radioresistance via ROD, but perhaps benefitting from elevated immune response [34,64,65] or differences in radiochemistry as compared to normal tissues. An interesting recent proposal by Spitz and colleagues has noted that differences between levels of labile iron and metabolism between normal and tumour tissues may give rise to different levels of radiation-induced organic hydroperoxides, and hence, DNA damage [3]. More work needs to be done, however, to verify the apparent iso-efficacy for more clinically relevant doses and fractionation schedules, and for tumours exhibiting different levels of measured hypoxia.
Beyond the modelling assumptions detailed in the Methods section and the parameter sensitivity analysis presented above, a major assumption of our work concerned the relevance of linear-quadratic model cell survival calculations for normal tissue toxicity. Other computational studies have likewise considered the impact of FLASH on the OER and survival [2] as well as just the OER [2,4]. A major potential limitation of these studies, including our own, is the assumed correlation between cell survival and normal tissue toxicity. Normal tissue toxicities are believed to be multi-factorial, possibly depending on immune response [66] and damage by radiotherapy-induced free radicals [67], in addition to radiation-induced DNA damage and subsequent cell killing. There is no clear relationship between cell survival and the probability of normal tissue toxicity. Despite this, iso-toxicity fractionation regimes have been shown to be well-described by the linear quadratic model [38], providing evidence that target cell survival is likely to play a major, if not the dominant, role in toxicity.
A recent alternative explanation for the potential large difference in cell survival curves (FLASH versus conventional) at high oxygen tensions has been proposed by Petersson and colleagues [68]. Neglecting the spatial heterogeneity of oxygen considered in the present work, they developed a novel model of ROD kinetics and found a maximum separation between mean oxygen enhancement ratios at relatively high oxygen partial pressures, ≳ 20 mmHg [68]. This is important since, as with our work, it provides a mechanism by which oxygen enhancement ratio effects can be significant even in nominally well-oxygenated tissues.

Conclusions
Modelling the heterogeneity in oxygen partial pressures that are expected to arise in tissues, we demonstrated that FLASH radiotherapy can substantially increase target cell survival as compared to conventional radiotherapy, even though the reduction in mean tissue oxygen levels was not nearly large enough to completely deplete tissue of oxygen. This effect was more pronounced for better-oxygenated normal tissues with proportionately smaller hypoxic fractions than hypoxic tumours, possibly explaining the observed differential impact of FLASH on normal tissue toxicities and tumour response. Our results highlight the importance of quantifying intra-tissue oxygen heterogeneity to fully assess the impact of FLASH on tissue radiation response.    The rate of oxygen metabolism for various normal tissues in units of μlO2/gtissue/min are shown in Table SM1. The corresponding rates in units of mmHg/s are also shown, using the same conversion factor applied by Thomlinson   As noted by Vaupel and colleagues [8], tumours generally exhibit a level of oxygen metabolism that is intermediate between the fastest oxygen-metabolizing tissues such as liver and the slowest metabolizing tissues, such as skeletal muscle at rest. Although considerable variation exists in the metabolic rate between different tumours [4,9], the "average" tumour also has an oxygen metabolic rate comparable to that of intestine and brain.
Based on these literature values, we use the metabolic rates cmax = 5, 15, and 40 mmHg/s in the simulations described in the main text, encompassing the expected rates for normal and tumour tissues (Table SM1).

Survival fraction curves for different values of and ⁄ .
Survival fraction curves for the random and lattice capillary model and different rates of ROD and oxygen metabolism are shown in Figure 4 in the main text using the radiobiological parameters ⍺ = 0.3 Gy -1 and ⍺/β = 3 Gy. Results for different parameter values are qualitatively similar. In Figure SM1, we show results for the random capillary model, metabolic rate cmax = 15 mmHg/s, and both rates of ROD studied. Note that the curves are qualitatively similar to each other, a result of the SF for D = 20 Gy being closely approximated by the fraction of tissue near zero oxygen partial pressure, irrespective of the values of and ⁄ [see Equation (SM3)].

Analytic formula for differential cell survival curves in the random capillary model
For a random distribution of capillaries in two dimensions, the fraction HF Λ of tissue that has an oxygen tension less than a value Λ can be approximated from Poisson statistics: Here, is the distance from a capillary at which oxygen levels are reduced to Λ by diffusion and metabolism, is the areal capillary density, O 2 is the diffusivity of oxygen, c is the capillary oxygen tension, max is the rate of oxygen consumption, and Λ is a constant of order unity that depends on the chosen threshold. For e.g., Λ = 0, Λ = 4 in two spatial dimensions [1]. Increasing Λ, Λ will decrease from this value. The exponent in Equation (SM1) represents the expected fractional value of tissue that is within a distance Λ of a capillary. HF Λ is thus the fraction of tissue that is more than Λ from the nearest capillary and hence, hypoxic.
Choosing Λ ~ 0-5 mmHg to coincide with the tension below which the OER changes most rapidly and further assuming that only these most hypoxic cells survive radiotherapy, the survival fraction is approximated by HF Λ : To the extent that the Michaelis-Menten parameters and ROD are not very different in magnitude, FLASH radiotherapy can further be approximated (see Equation 5 in the main text) as a reduction in metabolism, max → max + 0 ( ⁄ ) ≡ ̃m ax