There is increasing evidence that high doses of radiotherapy, like those delivered in stereotactic body radiotherapy (SBRT), trigger indirect mechanisms of cell death. Such effect seems to be two-fold. High doses may trigger an immune response and may cause vascular damage, leading to cell starvation and death. Development of mathematical response models, including indirect death, may help clinicians to design SBRT optimal schedules. Despite increasing experimental literature on indirect tumor cell death caused by vascular damage, efforts on modeling this effect have been limited. In this work, we present a biomathematical model of this effect. In our model, tumor oxygenation is obtained by solving the reaction–diffusion equation; radiotherapy kills tumor cells according to the linear–quadratic model, and also endothelial cells (EC), which can trigger loss of functionality of capillaries. Capillary death will affect tumor oxygenation, driving nearby tumor cells into severe hypoxia. Capillaries can recover functionality due to EC proliferation. Tumor cells entering a predetermined severe hypoxia status die according to a hypoxia-death model. This model fits recently published experimental data showing the effect of vascular damage on surviving fractions. It fits surviving fraction curves and qualitatively reproduces experimental values of percentages of functional capillaries 48 hours postirradiation, and hypoxic cells pre- and 48 hours postirradiation. This model is useful for exploring aspects of tumor and EC response to radiotherapy and constitutes a stepping stone toward modeling indirect tumor cell death caused by vascular damage and accounting for this effect during SBRT planning.

Significance:

A novel biomathematical model of indirect tumor cell death caused by vascular radiation damage could potentially help clinicians interpret experimental data and design better radiotherapy schedules.

Developments in radiotherapy treatment delivery have led to a situation where high-dose fractions can be safely delivered. This has boosted the use of stereotactic body radiotherapy (SBRT) and stereotactic radiation surgery (SRS; refs. 1–3). There is an intense debate on the radiobiology of large radiotherapy doses (4–9). It is becoming clear that, in addition to direct cell death due to damage to the cell DNA, high doses may cause indirect cell death. Such effect seems to be 2-fold: (i) high doses may trigger an immune response against the tumor (8–15), and (ii) may produce vascular damage, causing cell starvation and death (6, 8, 9, 16). A seminal experimental study showed that damage of capillaries triggered tumor cell death and increased control (17). Although there is evidence of vascular damage at high doses (18) and induced indirect cell death (16), the role of these effects on tumor control is still debatable: A critical analysis by Brown and colleagues (6) showed that control obtained with SBRT can be explained by dose escalation, with no need for new biology, and recent preclinical studies claimed that endothelial cell (EC) death does not play a significant role in tumor control after investigating genetically modified tumor mouse models with hyper-radiation–sensitive ECs (19, 20). There is also evidence of radiation-induced bystander effects contributing to the overall radiotherapy response (21).

Despite the complexity and current uncertainty about the involved mechanisms, it would be highly desirable to have reliable radiobiological models to help experimentalists, and to assist clinicians to design SBRT/SRS optimal schedules, similarly to the use of the LQ model in conventional radiotherapy. Serre and colleagues (22) have developed biomathematical models for the immune response triggered by radiotherapy, which reproduce experimental data and may prove useful for a better understanding of this effect (23). Despite the increasing experimental literature on indirect tumor cell death caused by vascular damage, efforts on modeling this effect have been limited. Recently, Merrem and colleagues at DKFZ have developed an analytical model of capillary death based on the inactivation of functional subunits (24). Paul-Gilloteaux and colleagues (25) have developed an in silico tumor response model that includes EC damage and its effect in radiotherapy response. In this work, we present a dynamical model of indirect tumor cell death triggered by damage to the tumor vascular system. It includes the interplay between oxygenation and direct tumor cell death, and EC death/proliferation, which can induce capillary death/recovery. Capillary death will cause tumor hypoxification, which competes with radiation-induced reoxygenation, and eventually may lead to indirect tumor cell death due to severe hypoxia. The model has been qualitatively validated by fitting experimental data reported in ref. 16. The model seems useful to explore aspects of tumor and EC response to radiotherapy and its effect on tumor control, and aims to be a stepping stone toward fully modeling indirect tumor cell death caused by high-dose radiotherapy.

Overview of the model: simplifications and limitations

We have developed a dynamical model of indirect tumor cell death caused by damage to the tumor vascular system. The flowchart of the model is presented in Fig. 1. The preirradiation oxygenation of the tumor is calculated by solving the reaction–diffusion equation with a finite-element-method (FEM), with a capillary network created according to a certain vascular fraction and capillary radii distribution (refer section “Oxygen distribution in the tumor and FEM solver”). This oxygen distribution is evaluated: Regions with severe hypoxia are treated like a preexisting necrotic core, where cells are dead, local consumption is adjusted accordingly, and the oxygen distribution is recalculated. Then, oxygenation-dependent α and β values are calculated by including oxygen enhancement ratios (OER). The tumor is irradiated, and tumor cells die according to the linear–quadratic (LQ) model, which is referred to as direct death (refer section “Tumor cell death, kinetics, and reoxygenation”). Doomed cells do not instantly, though, but follow a given death kinetic function after a mitotic delay. Non-damaged cells can proliferate and repopulate the tumor. Radiotherapy also kills ECs forming the vasculature of the tumor. If the death of ECs in any given capillary fulfills some predetermined conditions according to a capillary-death model, the capillary becomes non-functional and stops oxygenating the tumor, which we will refer to as capillary death (refer section “EC and capillaries' death”). This will affect tumor oxygenation, driving tumor cells in the vicinity of non-functional capillaries into severe hypoxia. This effect competes with reoxygenation caused by tumor cell death. Capillaries may recover function in time due to ECs proliferation (refer section “EC and capillaries' death”). Tumor cells entering a predetermined severe hypoxia status are tagged and killed according to a hypoxia-death model (refer section “Indirect tumor cell death due to severe hypoxia”). Finally, direct and indirect death, and overall surviving fractions are evaluated.

Figure 1.

Sketch of the flowchart of the model. Step 1, computation of baseline oxigenation map (t = 0) with the finite element method (FEM) solver; Step 2, dose delivery, computation of surviving fraction of tumor cells (including OERs) and ECs, SFD and SFe, respectively; Step 3, evaluation of direct and indirect cell death of tumor cells at time t; Step 4, evaluation of EC death, proliferation, and capillary death/recovery at time t; Step 5, recomputation of oxigenation map; Step 6, update of oxygenation status of cells; Step 7, if the time t has reached the maximum simulation time (TMAX), the simulation stops and data are saved, otherwise t = t+dt and return to Step 3.

Figure 1.

Sketch of the flowchart of the model. Step 1, computation of baseline oxigenation map (t = 0) with the finite element method (FEM) solver; Step 2, dose delivery, computation of surviving fraction of tumor cells (including OERs) and ECs, SFD and SFe, respectively; Step 3, evaluation of direct and indirect cell death of tumor cells at time t; Step 4, evaluation of EC death, proliferation, and capillary death/recovery at time t; Step 5, recomputation of oxigenation map; Step 6, update of oxygenation status of cells; Step 7, if the time t has reached the maximum simulation time (TMAX), the simulation stops and data are saved, otherwise t = t+dt and return to Step 3.

Close modal

For the sake of computational simplicity, we assume that the vascular fraction is homogeneous within the tumor (thus we can simulate only one voxel of the tumor). Also, we restrict our study to single-fraction treatments. Other limitations/simplifications of the model include: (i) The reaction/diffusion equation is solved in 2D, because this problem can be solved with a precise FEM (whereas a 3D-FEM would be unmanageable), and the 2D problem has been shown to be a good surrogate of the 3D problem for particular vascular networks. Nonetheless, the 2D approximation may miss some effects caused by the chaotic vascular networks of tumors; (ii) we use simple tumor vascular networks: We use realistic distributions of capillary radii, but the 2D solver implicitly assumes the vascular network to consist of parallel capillaries, and we assume the same length for all capillaries, oxygen partial pressure was considered equal in each capillary, 40 mm Hg, but real tumors may present a heterogeneity of values; (iii) There are very limited data on the dose–volume response of capillaries. We have considered a parallel/serial model, combining a model of capillary collapse in which accumulation of damage (EC death) leads to capillary pressure decreasing below the critical closing pressure, with a more geometrical, wire-severance model, accounting for clustered damage inspired on (24). Even though we obtain good fits with this model, the dose–response of capillaries might be different. (iv) Cell diffusion and tumor shrinkage have not been included at this stage, which, although a limitation, should be less important for single-dose schedules; (v) We model direct death with the LQ model, but at high doses repair saturation may limit the application of the LQ model (even though this effect is expected to be small compared with indirect death effects; refs. 4, 5).

Oxygen distribution in the tumor and FEM solver

The stationary state of the reaction–diffusion equation for oxygen distribution is (26, 27):

formula

In the above equation u(x, t) is the distribution of the oxygen, x and t denote spatial coordinates and time, Δ is the Laplacian operator, D the diffusion coefficient, and g is the oxygen consumption rate, which follows the Michaelis–Menten kinetics (28):

formula

where gmax is the maximum oxygen consumption rate, and k is the oxygen pressure for half-maximum consumption rate. The variable grel(x, t) accounts for changes in consumption due to death of cells, and will be defined later. We solve this equation in two spatial dimensions, which was shown to be a good surrogate of the 3D problem for vascular networks consisting of parallel capillaries (28). We have implemented a Matlab-based FEM solver. Blood capillaries are assumed to be cylindrical, with a heterogeneous distribution of radii. There are experimental studies showing lognormal and exponential distributions of radii of tumor capillaries (29–31), and we have implemented both options in our model. Minimum and maximum radii limits are also imposed according to experimental studies (Supplementary Table S1). Blood capillaries are randomly positioned in a voxel according to a vascular fraction, vf, defined as the ratio of capillary area to total area (28). The overlap of capillaries is explicitly avoided when generating the geometry. A steady 40 mm Hg oxygen partial pressure is associated to capillaries by imposing Dirichlet conditions on their boundaries. Homogeneous Neumann conditions are applied to the flux of u at the outer boundary of the voxel. If the function of a capillary is destroyed by radiotherapy at any time, that capillary is removed from the set of Dirichlet boundary conditions and it does not behave as an oxygen source anymore. The geometry is finely meshed, typically >10,000 nodes for a 1 × 1 mm pixel (Supplementary Fig. S1).

More detailed information about the simulation of oxygen distribution is provided in the Supplementary Materials.

Tumor cell death, kinetics, and reoxygenation

Direct tumor cell death is modeled with the LQ model, with no repair saturation:

formula

where SFd is the surviving fraction, d the radiotherapy dose, and α and β depend on the oxygenation of the cells as in ref. 32:

formula
formula

αox, βox, OERαm, and OERβm are the α-β values and the maximum enhancement ratios under fully aerobic conditions. Km controls the slope of the curve.

The delivery of a given dose will create a population of doomed cells, Nd = N0(1-SFd), whereas the population of undamaged cells will be N = N0 SFd, where N0 is the total number of tumor cells before irradiation. These population will evolve with time as follows:

formula
formula

Doomed cells that have been lethally damaged will progressively die according to a given kinetic model, whereas undamaged cells will be capable of proliferation. In this work, we have used exponential cell death and proliferation terms to describe this effect as in (33). Both cell death and proliferation start after a time τ, the mitotic delay caused by radiotherapy. The proliferation rate can depend on the level of oxygenation of cells, u, but there are some conflicting experimental reports on this issue, with some works (34) pointing to an acceleration of proliferation under hypoxic conditions, and some pointing the opposite (35). In this work, for the sake of simplicity we ignore possible dependences of the proliferation rate on the oxygenation level, but we assume proliferation arrest for cells under severe hypoxia, as defined in section “Indirect tumor cell death due to severe hypoxia.”

Doomed cells that disappear from the system, as well and newly generated cells are taken into account at each time step to rescale oxygen consumption in Eq. A.

EC and capillaries' death

ECs conforming blood capillaries are assumed to die following the LQ model (OERs are not considered as these cells are well oxygenated):

formula

We consider that damaged ECs do not die instantly after irradiation, but follow a given kinetic curve. According to some experimental results (36–38), we have modeled such kinetics as a linear decrease with time. In addition, non-doomed ECs can proliferate. We assume this proliferation to start after the end of the death kinetics, and to follow a Gompertz-like function. Therefore, we can write the evolution of the number of ECs, expressed as the surviving function of initial ECs, as:

formula

In this equation, SFe(t) shows the evolution of the relative number of surviving ECs with time. The above equation is used to stochastically evaluate death/recovery of ECs at each time step during the simulation.

Capillaries are assumed to have cylindrical geometry with length L and radius R. For simplicity, we assume a network of parallel capillaries, with all capillaries having the same length, but different radii according to a probability distribution. ECs of size dx*dy form the capillary wall. Thus, each capillary can be represented as a matrix with size dependent on L, R, dx, and dy (Supplementary Fig. S2). A matrix describing the state of each EC (1 = alive, 0 = death) is created at the beginning of the simulation. After irradiation, each EC is randomly tagged as lethally damaged or undamaged, according to Eq. H. Cells tagged as lethally damaged are then monitored during the simulation, and at each time step they are killed according to a probability distribution constructed from Eq. I. It seems important to notice that tumor cells can infiltrate capillaries and form part of the capillary wall. Although this effect could be incorporated in our models by tagging certain elements of the capillary matrix as tumor cells rather than ECs, and using the relevant surviving fractions to evaluate tumor cell death, we have ignored it for the sake of simplicity.

We present a model of capillary collapse that is based on the classical parallel/serial radiobiological response of organs–tissues to radiotherapy. The model has two contributions:

Parallel contribution.

Death of ECs increases capillary permeability, which reduces internal capillary pressure and can lead to capillary occlusion if the pressure goes below the capillary critical closing pressure (39). We have modeled this by using classical fluid mechanics. The flow in a capillary i is given by:

formula

In capillary bifurcations the flow is split according to the radii of the capillaries (40). Similarly, we consider that a dead EC causes a leak that takes away part of the capillary flow, proportional to a parameter r0. The change in flow is,

formula

and the flow when considering all dead ECs in a capillary is:

formula

This causes a decrease of the pressure in the capillary, which can be written using Poiseuille law as:

formula

If the pressure in the capillary, given by 40-ΔPi (40 mm Hg is the pressure that we assume in every functional capillary), gets below the critical closing pressure, the capillary will collapse. The critical closing pressure is given by the external interstitial pressure, PI, minus the wall pressure, PW, which is the excess pressure that the capillary wall can hold (41).

formula

Serial contribution.

If damage is clustered in a small section of the capillary (which we call functional sub-unit, FSU), we will consider that this leads to loss of capillary function. FSUs are defined as cross sections with n cells thickness. The capillary is killed if (and when) any FSU conforming it has been depleted (notice the overlap of different FSUs; capillary of length m cells will have m-n+1 FSUs of length n). This model is similar to that presented in (24), which instead used rectangular-shaped FSUs to model capillary death.

Lost capillary function can be restored due to EC proliferation, which was described above. Newly generated cells fill the EC state matrix, and the function of capillaries is restored if these new cells cancel out the capillary death conditions. No new capillaries are created due to EC proliferation, only recovery of already existing capillaries is allowed. Capillary function can also be restored by changes in the interstitial pressure, which we consider to be dependent on the number of cells in the voxel:

formula

where N(t) and Nd(t) are the number of tumor cells at time t, and N0 is the initial number of tumor cells.

Indirect tumor cell death due to severe hypoxia

A model of indirect tumor cell death due to severe hypoxia has been developed on the basis of experimental results (42, 43). We consider that cells below a certain oxygen threshold (uth) can die due to the effect of severe hypoxia. The probability of death depends on the time that the cells suffer severe hypoxia. This relationship is well described by an exponential law (43):

formula

In the above equation, SFH represents the surviving fraction of a population of cells exposed to severe hypoxia during a time th (t1/2H = log(2)/μth). Again, we use the symbol to denote that this surviving fraction refers to when all lethally damaged cells have died. Cells exposed to severe hypoxia do not die instantly. In fact, experimental results show that cells continue to die even after normal oxygenation has been restored (42). We have used a simple linear kinetic model to account for the death of cells that are (or have been) exposed to severe hypoxia:

formula

Time t' in the above equation is measured from the instant when cells enter in severe hypoxia. tH is the time at which we consider that every lethally damaged cell has died. Nodes entering/leaving severe hypoxia are tagged at each step of the calculation, and cells are killed according to the model above. There is one particularity with the implementation that is worth noticing here: As there is no way to know a priori how long cells will stay in severe hypoxia, the value of th (and SFH) in Eq. Q is reevaluated at each step of the simulation. We have seen numerically that this limitation does not affect much the shape of the curve.

Direct and indirect cell death of tumor cells will affect oxygen consumption and reoxygenation, as discussed in sections “Oxygen distribution in the tumor and FEM solver” and “Tumor cell death, kinetics, and reoxygenation.” Oxygen consumption is thus scaled as:

formula

Most relevant experimental parameters

There are conflicting experimental reports on the radiosensitivity of ECs, which is a critical parameter of our model: In vitro clonogenic assays point out a high radiosensitivity (αe ∼ 0.18 Gy; ref. 44), but in vivo experimental analyses show highly radio-resistant cells (αe ∼ 0.007 Gy; refs. 36, 37, 45). We have investigated different radiosensitivities of ECs. Li and colleagues (37) observed progressive death of ECs after a large single dose of radiotherapy until 16 hours postirradiation, whereas Peña and colleagues (36) observed a peak of apoptosis of ECs 12 hours postirradiation, and much lower values of apoptosis 24 hours postirradiation. Therefore, a value of te around 15 to 20 hours seems appropriate. Different lognormal and exponential distributions of capillary radii have been obtained from several publications (29–31). Regarding EC proliferation, doubling times in the range of 2–10 days have been reported (18).

Parameters a and b in the capillary death model can be calculated for ideal capillaries as a function of physical characteristics of blood and capillaries. However, we will consider them free parameters, with values ranging around a∼104 s−1 and b∼10−8 mm Hg m s. Values in the μm range were considered for r0, and different values of n were investigated according to the radiosensitivity of ECs. Interstitial pressure and wall pressure were taken as 10 and 3 mmHg, like in ref. 41.

There is not much experimental evidence on the value of uth. Some studies claim that below 1 mmHg cells can enter cell cycle arrest and eventually die (46, 47). On the other hand, Ljungkvist and colleagues (43) observed hypoxia turnover for cells with oxygen levels below 10 mmHg (cells marked with pimonidazole and CCI103F). We have not used such high thresholds, as they would result in a large percentage of cells affected by hypoxia death at baseline for standard vascular fractions. Values of t1/2H have been mainly obtained from (42, 43). The latter study reported values ≈15 to 50 hours, whereas the former observed progressive cell death 24 hours after cells had returned to normoxic status, therefore the term tH in our model to account for such effect. We have used tH = 3t1/2H. We have considered oxygen-dependent radiosensitivities for tumor cells (32), and a half-life of doomed cells (t1/2d = log(2)/λ) of 120 hours (33, 48). Reported OERαm and OERβm values are usually in the 2 to 3 range. We have used a single OERm for α and β sensitization, which is supported by studies like (49), and set its value to 2.5. Mitotic delays around 1 day are commonly reported, so we have taken that value.

Song's experimental data

In ref. 16, tumors grown in mice were irradiated with single doses (10–30 Gy), and surviving fractions evaluated at different times postirradiation. They observed a marked decrease of the surviving fraction with time, which they attributed to indirect tumor cell death caused by vascular damage. Further analysis showed reduced blood perfusion and increased hypoxia in the tumor microenvironment. We have extracted SFs at different times postirradiation and used our model to fit the data. We also extracted data of Hoechst 33342 and pimonidazole uptake preirradiation and 48 hours after 20 Gy irradiation.

To fit experimental data, surviving fractions at each time t in our model were obtained by considering both dead and lethally damaged cells at that time (we assume that lethally damaged cells will not proliferate when plated). Therefore, this will include all cells killed or lethally damaged by radiotherapy, which is delivered at time t = 0, and cells killed or lethally damaged by severe hypoxia after irradiation. On the other hand, we have considered percentage of areas positive for Hoechst reported in ref. 16 as a surrogate of capillary function. To obtain values that we can compare with experimental pimonidazole positive areas, we assume that pimonidazole can diffuse a certain diffusion distance, dpimo, from functional vessels, and that areas within that range with living cells and u<10 mmHg will uptake pimonidazole (Supplementary Fig. S3; ref. 43).

We simulated a tumor with a homogeneous 4% vascular fraction and a lognormal distribution of capillary radii with minimum and maximum values of 2.8 and 26.7 μm, taken from ref. 30. The rest of the parameters took the values reported in the previous section, or where qualitatively optimized within those ranges to obtain good fits to experimental data. Parameter values are reported in Table 1.

Table 1.

List of experimental and modeling parameters of used to fit experimental data in Fig. 2 

Tumor cellsECs and capillaries
DefinitionSymbolValueDefinitionSymbolValue
Alpha (LQ) for tumor αox 0.22 Gy−1 Alpha (LQ) for ECs αe 0.009 Gy−1 
Mitotic delay τ 24 h Maximum survival time for lethally damaged EC te 20 h 
Proliferation rate of tumor cells μ 0.0058 h−1 Proliferation rate of ECs λ 0.0046 h−1 
Maximum oxygen enhancement ratio OERαm OERβm 2.5 Length of capillaries L 300 μm 
Sensitization curve slope parameter Km 3.28 mmHg Size of capillary cells dx, dy 10 μm 
Half-life of radiation lethally damaged cells t1/2d 120 h Capillary size distribution  Lognormal 
Reaction/diffusion equation Radii of capillaries (min, max)  ≈3-27 μm 
Definition Symbol Value Capillary leakage due to EC death r0 5 μm 
Maximum oxygen consumption rate gmax 15 mmHg s−1 Parameter linking capillary flow and radius a 4 × 10−14 s−1 
Oxygen pressure for half-maximum consumption rate k 2.5 mmHg Parameter linking flow and pressure b 1.7 × 1016 mmHg μm s 
Oxygen pressure in capillaries ucapillary 40 mmHg Interstitial pressure PI 10 mmHg 
Oxygen diffusion coefficient D 2 × 10−9 m2 s−1 Capillary wall excess pressure PW 3 mmHg 
Severe hypoxia death length of the function subunit n 
Definition Symbol Value Maximum diffusion distance of pimonidazone dpimo 125 μm 
Severe hypoxia oxygen threshold uth 1.5 mmHg    
Half-life of cells in severe hypoxia t1/2H 20 h    
Maximum surviving time for cells lethally damaged by hypoxia tH 3 t1/2H    
Tumor cellsECs and capillaries
DefinitionSymbolValueDefinitionSymbolValue
Alpha (LQ) for tumor αox 0.22 Gy−1 Alpha (LQ) for ECs αe 0.009 Gy−1 
Mitotic delay τ 24 h Maximum survival time for lethally damaged EC te 20 h 
Proliferation rate of tumor cells μ 0.0058 h−1 Proliferation rate of ECs λ 0.0046 h−1 
Maximum oxygen enhancement ratio OERαm OERβm 2.5 Length of capillaries L 300 μm 
Sensitization curve slope parameter Km 3.28 mmHg Size of capillary cells dx, dy 10 μm 
Half-life of radiation lethally damaged cells t1/2d 120 h Capillary size distribution  Lognormal 
Reaction/diffusion equation Radii of capillaries (min, max)  ≈3-27 μm 
Definition Symbol Value Capillary leakage due to EC death r0 5 μm 
Maximum oxygen consumption rate gmax 15 mmHg s−1 Parameter linking capillary flow and radius a 4 × 10−14 s−1 
Oxygen pressure for half-maximum consumption rate k 2.5 mmHg Parameter linking flow and pressure b 1.7 × 1016 mmHg μm s 
Oxygen pressure in capillaries ucapillary 40 mmHg Interstitial pressure PI 10 mmHg 
Oxygen diffusion coefficient D 2 × 10−9 m2 s−1 Capillary wall excess pressure PW 3 mmHg 
Severe hypoxia death length of the function subunit n 
Definition Symbol Value Maximum diffusion distance of pimonidazone dpimo 125 μm 
Severe hypoxia oxygen threshold uth 1.5 mmHg    
Half-life of cells in severe hypoxia t1/2H 20 h    
Maximum surviving time for cells lethally damaged by hypoxia tH 3 t1/2H    

Sensitivity analysis

We followed the Sobol methodology to calculate the parametric sensitivity of the model: First-order and total effect Sobol indices were calculated following ref. 50. More information on this methodology is provided in the Supplementary Materials.

Implementation

The model was implemented in different routines in Matlab (The Matworks). The model is computationally demanding, taking around 15 minutes to perform a simulation for one-dose value with Tmax = 120 hours and dt = 2 hours.

Fit to Song's experimental data

Fitting SFs at t = 0 to the LQ model showed effective α and β values of 0.13 Gy−1 and 0.002 Gy−2, respectively, but with very large uncertainties (∼100%) due to noisy experimental data. As those values are for tumors with heterogeneous oxygenation, and the input parameters of our model are αox and βox, fitted values were scaled up to account for hypoxia radio-resistance. In addition, as the ratio α/β seems to be very large, we have taken βox = 0 to reduce the number of free parameters. Regarding the radiosensitivity of ECs, we have explored radiosensitivities around in vivo experimental results (αe∼0.007 Gy; refs. 34, 35, 40), and for the same reason as above we have taken βe = 0.

Figure 2 shows the fit of our model to the surviving fraction versus time data, with parameters αox = 0.22 Gy−1, αe = 0.009 Gy−1, t1/2H = 20 hours, and proliferation rates of μ = 0.0058 h−1 and λ = 0.0046 h−1 for tumor cells and ECs, respectively. Best fitting parameters are reported in Table 1. We report mean and standard deviations of 25 simulations, varying parameters according to a Gaussian distribution with standard deviation equal to 10% of the mean, to study the effect of uncertainties in the parameters.

Figure 2.

Surviving fractions at different times postirradiation with different doses. Experimental data taken from ref. 16. Experimental data points (o) and model results: mean (solid lines) and mean ± 1 SD of 25 simulations considering Gaussian uncertainties of the parameters, with 10% SD around best fitting parameters (dashed lines).

Figure 2.

Surviving fractions at different times postirradiation with different doses. Experimental data taken from ref. 16. Experimental data points (o) and model results: mean (solid lines) and mean ± 1 SD of 25 simulations considering Gaussian uncertainties of the parameters, with 10% SD around best fitting parameters (dashed lines).

Close modal

The decrease in surviving fraction with time is due to EC death and capillary collapse, which causes severe hypoxia in the tumor. Surviving fractions of ECs with these parameters are 0.91, 0.87, 0.83, and 0.76 for 10, 15, 20, and 30 Gy respectively, which results in severe damage to the capillary network (as shown in Fig. 3), leading to important indirect cell death due to vascular damage. With the employed parameters, capillary death is dominated by the parallel contribution. There is vascular damage even for 10 Gy, but it does not result in important indirect cell death because the remaining capillaries can maintain adequate oxygenation. Capillary functionality is progressively recovered in time, due to EC proliferation.

Figure 3.

Functionality of capillaries versus time postirradiation obtained with our model. Only results for one simulation with best fitting parameters are reported for the sake of graphical clarity. The experimental point shown was taken from ref. 16 and corresponds to 48 hours postirradiation with 20 Gy.

Figure 3.

Functionality of capillaries versus time postirradiation obtained with our model. Only results for one simulation with best fitting parameters are reported for the sake of graphical clarity. The experimental point shown was taken from ref. 16 and corresponds to 48 hours postirradiation with 20 Gy.

Close modal

The percentage of functional capillaries 48 hours after irradiation with 20 Gy is 14% ± 9%, not far from the ratio of Hoechst positive areas of 8% ± 3% extracted from ref. 16. On the other hand, in our simulation pimonidazole positive areas increase from 20% ± 5% before irradiation to 25% ± 11% 48 hours after irradiation, although the difference is not significant. Experimental pimonidazole positive area increases from 17% ± 7% to 21% ± 3.5% in ref. 16.

In Fig. 4, we show a qualitative representation of the spatio-temporal evolution of oxygenation, and surviving fractions associated to direct (dose) and indirect death (severe hypoxia) in our model, for an irradiation with 15 Gy. At t = 0 the voxel has a heterogeneous distribution of oxygenation, including a necrotic core where u < uth. Twenty-four hours after irradiation most capillaries are nonfunctional, and many regions are in severe hypoxic conditions. Cells in those regions are already dying, creating a heterogeneous map of indirect cell death (lower surviving fractions around nonfunctional capillaries). On the other hand, the map of direct cell death also presents heterogeneities, due to heterogeneities in oxygenation at the time of irradiation. Ninety-six hours after irradiation many damaged capillaries have recovered function due to EC proliferation. That, and the death of many cells, leads to a fairly well oxygenated voxel. On the other hand, heterogeneities in direct, and especially indirect cell death are more noticeable.

Figure 4.

Qualitative representation of the spatiotemporal evolution of oxygenation (left), and surviving fractions associated with direct death (middle) and indirect death (right). White circles are capillaries. Notice the different colormap scale in the middle panels to highlight the effect of OERs on the surviving fraction.

Figure 4.

Qualitative representation of the spatiotemporal evolution of oxygenation (left), and surviving fractions associated with direct death (middle) and indirect death (right). White circles are capillaries. Notice the different colormap scale in the middle panels to highlight the effect of OERs on the surviving fraction.

Close modal

Sensitivity analysis

The most sensitive parameters of the model are: r0, αox, αe, vf, L, a, b. The model is mostly sensitive to parameters characterizing cell death (αox, αe), and those characterizing capillary death (αe, a, b, L), and therefore indirect hypoxia driven cell death. The strongest sensitivity corresponds to the parameter r0, which is expected due to the term r03 in Eq. K. Sobol indices are reported in Supplementary Table S2.

Indirect tumor cell death caused by vascular damage seems to play a role in high-dose hypofractionated radiotherapy. Modeling of this effect has not been extensively addressed in the literature despite the importance it might have for hypofractionated schedules. In the present study, we have developed a dynamic model of indirect cell death caused by tumor vascular damage after radiotherapy. Our model can qualitatively fit experimental data presented in ref. 16, in particular the decrease of surviving fraction with time observed in that study, which was attributed to indirect cell death. In our model, this trend is caused by EC and capillary death, which causes severe hypoxia in the tumor and induces indirect cell death. Interestingly, the fit of the SF versus time data also leads to a good match with other experimental parameters reported in ref. 16. The percentage of functional capillaries 48 hours after irradiation with 20 Gy obtained with our model is not far from the experimental ratio of Hoechst-positive areas; the percentages of area with u< 10 mm Hg and within diffusion distance from a functional vessel before irradiation and 48 hours after irradiation obtained with our model are similar to the percentages of pimonidazole positive areas at the same times.

Indirect cell death in our model is highly sensitive to the radiosensitivity of ECs and capillary death modeling, as shown by sensitivity analysis. There are conflicting experimental reports on EC radiosensitivity, with reported values ranging from highly radiosensitive to highly radioresistant (36, 37, 44, 45). In this work, we have presented fits obtained with highly radioresistant ECs (αe ∼ 0.01), values that are similar to those reported in in vivo experiments. However, we have explored different regimes, and similarly good results can be obtained with moderate radiosensitivities (αe ∼ 0.1), if the capillary-death parameters are adjusted accordingly. Also, there are very limited data on the dose-volume response of capillaries. We have considered a parallel/serial model, combining a model of capillary collapse due to capillary pressure decreasing below the critical closing pressure due to EC death, with a more geometrical, wire severance model, inspired by ref. 24. Model parameters can be adjusted to give more weight to one effect or another when more data on capillary response becomes available.

The set of experimental data used for validation was particularly valuable for that purpose, as the authors reported the kinetics of cell death associated to vascular damage, as well as indirect measurements of capillary functionally and evolution of oxygenation pre- and postirradiation. More experimental data like that presented in ref. 16 would be very important for further validation of the model.

Dosage and fractionation of SBRT/SRS schedules are empirically set, as experimental data and models suggesting optimal doses and number of fractions are scarce. It would be highly desirable to have reliable radiobiological models to help clinicians to design SBRT/SRS optimal schedules, similarly to the role of the LQ model in conventional radiotherapy. After proper experimental validation, one could envision models like the one presented here and experimental studies guiding the selection of optimal doses and fractionations at the population level, or even guiding the individualization of therapy (for example according to functional imaging maps of baseline hypoxia and tumor vasculature, or according to resistant vascular phenotypes of some patients, which could be used to identify patients more likely to benefit from hypofractionation). This is still far from achievable with the current understanding of indirect cell death. As models like this are probably too complex and computationally demanding to be directly employed in the clinic, an interesting approach would be to use such complex models and results obtained from experiments to develop simple phenomenological extensions of the LQ model. An interesting hypothesis that may be worth further investigation is that those effective models might lead to better fits of the dose–response curves of hypofractionated treatments than the classical LQ model.

It seems important to point out that some studies claim that indirect cell death due to vascular damage does not play a significant role in the outcome of SBRT/SRS treatments (6, 19, 20). It seems important to highlight that further experimental and modeling investigations seem necessary to shed light upon this controversy. Biomathematical models may be useful to interpret such data.

Indirect cell death due to radiation-induced immune response also seems to play a role, an effect that has been investigated by modelers (22, 23). Models like the one presented here could be coupled with those models to develop complete models of indirect cell death in radiotherapy and radio-immunotherapy. Also, as EC damage may be involved in radiation-induced immune response (15), models like this could prove useful in studying radiation-induced immune response.

No potential conflicts of interest were disclosed.

Conception and design: J. Pardo-Montero

Development of methodology: P. Rodríguez-Barbeito, P. Díaz-Botana, A. Gago-Arias, M. Feijoo, S. Neira, J. Guiu-Souto, Ó. López-Pouso, J. Pardo-Montero

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): P. Rodríguez-Barbeito, P. Díaz-Botana, A. Gago-Arias, M. Feijoo, J. Pardo-Montero

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): P. Rodríguez-Barbeito, P. Díaz-Botana, A. Gago-Arias, J. Pardo-Montero

Writing, review, and/or revision of the manuscript: P. Díaz-Botana, A. Gago-Arias, M. Feijoo, J. Guiu-Souto, Ó. López-Pouso, A. Gómez-Caamaño, J. Pardo-Montero

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): P. Rodríguez-Barbeito

Study supervision: J. Pardo-Montero

J. Pardo-Montero acknowledges the hospitality of Instituto de Física, Pontificia Universidad Católica de Chile, where part of this work was carried out (funded by Concurso VRI 2017 Visitas de Investigación de Profesores Extranjeros). A. Gago-Arias acknowledges the hospitality of Instituto de Investigación Sanitaria de Santiago during her stay at the institution. J. Pardo-Montero is funded by Instituto de Salud Carlos III (Miguel Servet II CPII17/00028, PI17/01428, and DTS17/00123 grants, FEDER cofunding). A. Gago-Arias is funded by Fondecyt (11170575 grant) and the EU H2020 program (Marie Sklodowska-Curie action TRIDOS, 839135). Ó. López-Pouso is supported by FEDER and Xunta de Galicia (GRC2013-014), and by the Spanish Ministry of Science, Innovation and Universities (MTM2017-86459-R).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

1.
Sperduto
PW
. 
A review of stereotactic radiosurgery in the management of brain metastases
.
Technol Cancer Res Treat
2003
;
2
:
105
110
.
2.
Chang
BK
,
Timmerman
RD
. 
Stereotactic body radiation therapy: a comprehensive review
.
Am J Clin Oncol
2007
;
30
:
637
44
.
3.
Lo
SS
,
Fakiris
AJ
,
Chang
EL
,
Mayr
NA
,
Wang
JZ
,
Papiez
L
, et al
Stereotactic body radiation therapy: a novel treatment modality
.
Nat Rev Clin Oncol
2010
:
7
;
44
54
.
4.
Kirkpatrick
JP
,
Meyer
JJ
,
Marks
LB
. 
The linear-quadratic model is inappropiate to model high dose per fraction effects in radiosurgery
.
Semin Radiat Oncol
2008
:
18
;
240
3
.
5.
Brenner
DJ
. 
The linear-quadratic model is an appropriate methodology for determining isoeffective doses at large doses per fraction
.
Semin Radiat Oncol
2008
:
18
;
234
9
.
6.
Brown
JM
,
Carlson
DJ
,
Brenner
DJ
. 
The tumor radiobiology of SRS and SBRT: are more than the 5 Rs involved ?
Int J Radiat Oncol Biol Phys
2014
;
88
:
254
62
.
7.
Song
CW
,
Park
H
,
Griffin
RJ
,
Levitt
SH
. 
Radiobiology of stereotactic radiosurgery and stereotactic body radiation therapy
. In:
Lewitt
SH
, et al
Technical basis of radiation therapy: practical clinical application
.
New York, NY
:
Springer Publishing Co
; 
2012
. p.
51
61
.
8.
Sperduto
PW
,
Song
CW
,
Kirkpatrick
JP
,
Glatstein
E
. 
A hypothesis: indirect cell death in the radiosurgery era
.
Int J Radiat Oncol Biol Phys
2015
:
91
:
11
13
.
9.
Song
CW
,
Kim
M
,
Cho
LC
,
Dusenbery
K
,
Sperduto
PW
. 
Radiobiological basis of SBRT and SRS
.
Int J Clin Oncol
2014
:
19
:
570
8
.
10.
Lee
Y
,
Auh
SL
,
Wang
Y
,
Burnette
B
,
Wang
Y
,
Meng
Y
, et al
Therapeutic effects of ablative radiation on local tumor require CD8+ T cells: changing strategies for cancer treatment
.
Blood
2009
:
114
:
589
95
.
11.
Filatenkov
A
,
Baker
J
,
Mueller
A
,
Kenkel
JA
,
Ahn
G-O
,
Dutt
S
, et al
Ablative tumor radiation can change the tumor immune cell microenvironment to induce durable complete remissions
.
Clin Cancer Res
2015
;
21
:
3727
39
.
12.
De La Maza
L
,
Wu
M
,
Wu
L
,
Yun
H
,
Zhao
Y
,
Cattral
M
, et al
In situ vaccination after accelerated hypofractionated radiation and surgery in a mesothelioma mouse model
.
Clin Cancer Res
2017
;
23
:
5502
13
.
13.
Wang
Y
,
Liu
Z
,
Yuan
H
,
Deng
W
,
Li
J
,
Huang
Y
, et al
The reciprocity between radiotherapy and cancer immunotherapy
.
Clin Cancer Res
2019
;
25
:
1709
17
.
14.
Herrera
FG
,
Bourhis
J
,
Coukos
G
. 
Radiotherapy combination opportunities leveraging immunity for the next oncology practice
.
CA Cancer J Clin
2017
;
67
:
65
85
.
15.
He
J
,
Yin
Y
,
Luster
TA
,
Watkins
L
,
Thorpe
PE
. 
Antiphosphatidylserine antibody combined with irradiation damages tumor blood vessels and induces tumor immunity in a rat model of glioblastoma
.
Clin Cancer Res
2009
;
15
:
6871
80
.
16.
Song
CW
,
Lee
YJ
,
Griffin
RJ
,
Park
I
,
Koonce
NA
,
Hui
S
, et al
Indirect tumor cell death after high-dose hypofractionated irradiation: implications for stereotactic body radiation therapy and stereotactic radiation surgery
.
Int J Radiat Oncol Biol Phys
2015
;
93
:
166
72
.
17.
García-Barros
M
,
Paris
F
,
Cordón-Cardo
C
,
Lyden
D
,
Rafii
S
,
Haimovitz-Friedman
A
, et al
Tumor response to radiotherapy regulated by endothelial cell apoptosis
.
Science
2003
;
300
:
1155
9
.
18.
Park
HJ
,
Griffin
RJ
,
Hui
S
,
Levitt
SH
,
Song
CW
. 
Radiation-induced vascular damage in tumors: implications of vascular damage in ablative hypofractionated radiotherapy (SBRT and SRS)
.
Radiat Res
2012
;
177
:
311
27
.
19.
Moding
EJ
,
Castle
KD
,
Perez
BA
,
Oh
P
,
Min
HD
,
Norris
H
, et al
Tumor cells, but not endothelial cells, mediate eradication of primary sarcomas by stereotactic body radiation therapy
.
Sci Transl Med
2015
;
7
:
278ra34
.
20.
Torok
JA
,
Oh
P
,
Castle
KD
,
Reinsvold
M
,
Ma
Y
,
Luo
L
, et al
Deletion of ATM in tumor but not endothelial cells improves radiation response in a primary mouse model of lung adenocarcinoma
.
Cancer Res
2018
;
79
:
773
82
.
21.
Blyth
BJ
,
Sykes
PJ
. 
Radiation-induced bystander effects: what are they, and how relevant are they to human radiation exposures?
Radiat Res
2011
;
176
:
139
57
.
22.
Serre
R
,
Barlesi
F
,
Muracciole
X
,
Barbolosi
D
. 
Immunologically effective dose: a practical model for immuno-radiotherapy
.
Oncotarget
2018
;
9
:
31812
19
.
23.
Serre
R
,
Benzekry
S
,
Padovani
L
,
Meille
C
,
André
N
,
Ciccolini
J
, et al
Mathematical modeling of cancer immunotherapy and its synergy with radiotherapy
.
Cancer Res
2016
;
76
:
4931
40
.
24.
Merrem
A
,
Bartzsch
S
,
Laissue
J
,
Oelfke
U
. 
Computational modelling of the cerebral cortical microvasculature: effect of x-ray microbeams versus broad beam irradiation
.
Phys Med Biol
2017
;
62
:
3902
22
.
25.
Paul-Gilloteaux
P
,
Potiron
V
,
Delpon
G
,
Supiot
S
,
Chiavassa
S
,
Paris
F
, et al
Optimizing radiotherapy protocols using computer automata to model tumour cell death as a function of oxygen diffusion processes
.
Sci Rep
2017
;
7
:
2280
.
26.
Daşu
A
,
Toma-Daşu
I
,
Karlsson
M
. 
Theoretical simulation of tumour oxygenation and results from acute and chronic hypoxia
.
Phys Med Biol
2003
;
48
:
2829
42
.
27.
Mönnich
D
,
Troost
EG
,
Kaanders
JH
,
Oyen
WJ
,
Alber
M
,
Thorwarth
D
. 
Modelling and simulation of [18F]fluoromisonidazole dynamics based on histology-derived microvessel maps
.
Phys Med Biol
2011
;
56
:
2045
57
.
28.
Espinoza
I
,
Peschke
P
,
Karger
CP
. 
A model to simulate the oxygen distribution in hypoxic tumors for different vascular architectures
.
Med Phys
2013
:
40
:
081703
29.
Konerding
MA
,
Malkusch
W
,
Klapthor
B
,
van Ackern
C
,
Fait
E
,
Hill
SA
, et al
Evidence for characteristic vascular patterns in solid tumours: quantitative studies using corrosion casts
.
Brit J Cancer
1999
:
80
:
724
32
.
30.
Konerding
MA
,
Fait
E
,
Gaumann
A
. 
3D microvascular architecture of pre-cancerous lesions and invasive carcinomas of the colon
.
Br J Cancer
2001
:
84
:
1354
62
.
31.
Tufto
I
,
Lyng
H
,
Rofstad
EK
. 
Vascular density in human melanoma xenografts: relationship to angiogenesis, perfusion and necrosis
.
Cancer Lett
1998
:
123
:
59
65
.
32.
Wouters
BG
,
Brown
JM
. 
Cells at intermediate oxygen levels can be more important than the "hypoxic fraction" in determining tumor response to fractionated radiotherapy
.
Radiat Res
1997
:
147
:
541
50
.
33.
Gago-Arias
A
,
Aguiar
P
,
Espinoza
I
,
Sánchez-Nieto
B
,
Pardo-Montero
J
. 
Modelling radiation-induced cell death and tumour re-oxygenation: local versus global and instant versus delayed cell death
.
Phys Med Biol
2016
:
61
:
1204
16
.
34.
Nordsmark
M
,
Høyer
M
,
Keller
J
,
Nielsen
OS
,
Jensen
OM
,
Overgaard
J
. 
The relationship between tumor oxygenation and cell proliferation in human soft tissue sarcomas
.
Int J Radiat Oncol Biol Phys
1996
;
35
:
701
8
.
35.
Jeong
J
,
Oh
JH
,
Sonke
JJ
,
Belderbos
J
,
Bradley
JD
,
Fontanella
AN
, et al
Modeling the cellular response of lung cancer to radiation therapy for a broad range of fractionation schedules
.
Clin Cancer Res
2017
;
23
:
5469
79
.
36.
Peña
LA
,
Fuks
Z
,
Kolesnick
RN
. 
Radiation-induced apoptosis of endothelial cells in the murine central nervous system: protection by fibroblast growth factor and sphingomyelinase deficiency
.
Cancer Res
2000
:
60
:
321
7
.
37.
Li
YQ
,
Chen
P
,
Haimovitz-Friedman
A
,
Reilly
RM
,
Wong
CS
. 
Endothelial apoptosis initiates acute blood brain disruption after ionizing radiation
.
Cancer Res
2003
:
63
:
5950
6
.
38.
Donker
M
,
Van Furth
WR
,
Mulder-Van Der Kracht
S
,
Hovinga
KE
,
Verhoeff
JJ
,
Stalpers
LJ
, et al
Negligible radiation protection of endothelial cells by vascular endothelial growth factor
.
Oncol Rep
2007
;
18
:
709
14
.
39.
Burton
AC
. 
On the physical equilibrium of small blood vessels
.
Am J Physiol
1951
:
164
:
319
29
.
40.
Sherman
TF
. 
On connecting large vessels to small. The meaning of Murray's law
.
J Gen Physiol
1981
:
78
:
431
53
.
41.
Milosevic
MF
,
Fyles
AW
,
Hill
RP
. 
The relationship between elevated interstitial fluid pressure and blood flow in tumors: a bioengineering analysis
.
Int J Radiat Oncol Biol Phys
1999
:
43
:
1111
23
.
42.
Steinbach
JP
,
Wolburg
H
,
Klump
A
,
Probst
H
,
Weller
M
. 
Hypoxia-induced cell death in human malignant glioma cells: energy deprivation promotes decoupling of mitochondrial cytochrome c release from caspase processing and necrotic cell death
.
Cell Death Differ
2003
:
10
:
823
32
.
43.
Ljungkvist
AS
,
Bussink
J
,
Kaanders
JH
, et al
Hypoxic cell turnover in different solid tumor lines
.
Int J Radiat Oncol Biol Phys
2005
:
62
:
1157
68
.
44.
Riquier
H
,
Wera
AC
,
Heuskin
AC
,
Feron
O
,
Lucas
S
,
Michiels
C
. 
Comparison of X-ray and alpha particle effects on a human cancer and endothelial cells: survival curves and gene expression profiles
.
Radiother Oncol
2013
;
106
:
397
403
.
45.
Lyubimova
NV
,
Coultas
PG
,
Yuen
K
,
Martin
RF
. 
In vivo radioprotection of mouse brain endothelial cells by Hoechst 33342
.
Br J Radiol
2001
;
74
:
77
82
.
46.
Harriss-Phillips
WM
,
Bezak
E
,
Yeoh
E
. 
The HYP-RT hypoxic tumour radiotherapy algorithm and accelerated repopulation dose per fraction study
.
Comput Math Methods Med
2012
;
2012
:
363564
.
47.
Ljungkvist
ASE
,
Bussink
J
,
Rijken
PF
,
Kaanders
JH
,
van der Kogel
AJ
,
Denekamp
J
. 
Vascular architecture, hypoxia, and proliferation in first generation xenografts of human head-and-neck squamous cell carcinomas
.
Int J Radiat Oncol Biol Phys
2002
:
54
:
215
28
.
48.
Li
W
,
Li
F
,
Huang
Q
,
Shen
J
,
Wolf
F
,
He
Y
, et al
Quantitative, noninvasive imaging of radiation-induced DNA double-strand breaks in vivo
.
Cancer Res
2011
;
71
:
4130
7
.
49.
Carlson
DJ
,
Stewart
RD
,
Semenenko
VA
. 
Effects of oxygen on intrinsic radiation sensitivity: a test of the relationship between aerobic and hypoxic linear-quadratic (LQ) model parameters
.
Med Phys
2006
:
33
:
3105
15
.
50.
Saltelli
A
,
Annoni
P
,
Azzini
I
,
Campolongo
F
,
Ratto
M
,
Tarantola
S
. 
Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index
.
Comput Phys Commun
2010
:
181
:
259
70
.