A Model of Indirect Cell Death Caused by Tumor Vascular Damage after High-Dose Radiotherapy

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.

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.

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).

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

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):

where g_max is the maximum oxygen consumption rate, and k is the oxygen pressure for half-maximum consumption rate. The variable g_rel(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.

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

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

α_ox, β_ox, OERα_m, and OERβ_m are the α-β values and the maximum enhancement ratios under fully aerobic conditions. K_m controls the slope of the curve.

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

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.

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

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:

In this equation, SF_e(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:

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:

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 r₀. The change in flow is,

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

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

If the pressure in the capillary, given by 40-ΔP_i (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, P_I, minus the wall pressure, P_W, which is the excess pressure that the capillary wall can hold (41).

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:

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

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 (u_th) 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):

In the above equation, SF_H^∞ represents the surviving fraction of a population of cells exposed to severe hypoxia during a time t_h (t_1/2^H = 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:

Time t' in the above equation is measured from the instant when cells enter in severe hypoxia. t_H^∞ 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 t_h (and SF_H^∞) 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:

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 t_e 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∼10⁴ s⁻¹ and b∼10⁻⁸ mm Hg m s. Values in the μm range were considered for r₀, 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 u_th. 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 t_1/2^H 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 t_H^∞ in our model to account for such effect. We have used t_H^∞ = 3t_1/2^H. We have considered oxygen-dependent radiosensitivities for tumor cells (32), and a half-life of doomed cells (t_1/2^d = 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 OER_m 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.

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, d_pimo, 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.

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.

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 T_max = 120 hours and dt = 2 hours.

Fitting SFs at t = 0 to the LQ model showed effective α and β values of 0.13 Gy⁻¹ and 0.002 Gy⁻², 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⁻¹, α_e = 0.009 Gy⁻¹, t_1/2^H = 20 hours, and proliferation rates of μ = 0.0058 h⁻¹ and λ = 0.0046 h⁻¹ 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.

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.

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 < u_th. 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.

The most sensitive parameters of the model are: r₀, α_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 r₀, which is expected due to the term r₀³ 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.