A Biomathematical Model of Tumor Response to Radioimmunotherapy with αpDL1 and αcTLA4

1545-5963 (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TCBB.2022.3174454, IEEE/ACM Transactions on Computational Biology and Bioinformatics 1 A Biomathematical Model of Tumor Response to Radioimmunotherapy with αPDL1 and αCTLA4 Isabel Gonz´ alez-Crespo, Antonio G´ omez-Caama˜ no, ´ Oscar L´ opez Pouso, John D. Fenwick, and Juan Pardo-Montero Abstract—There is evidence of synergy between radiotherapy and immunotherapy. Radiotherapy can increase liberation of tumor antigens, causing activation of antitumor T-cells. This effect can be boosted with immunotherapy. Radioimmunotherapy has potential to increase tumor control rates. Biomathematical models of response to radioimmunotherapy may help on understanding of the mechanisms affecting response, and assist clinicians on the design of optimal treatment strategies. In this work we present a biomathematical model of tumor response to radioimmunotherapy. The model uses the linear-quadratic response of tumor cells to radiation (or variation of it), and builds on previous developments to include the radiation-induced immune effect. We have focused this study on the combined effect of radiotherapy and αPDL1/αCTLA4 therapies. The model can fit preclinical data of volume dynamics and control obtained with different dose fractionations and αPDL1/αCTLA4. A biomathematical study of optimal combination strategies suggests that a good understanding of the involved biological delays, the biokinetics of the immunotherapy drug, and the interplay between them, may be of paramount importance to design optimal radioimmunotherapy schedules. Biomathematical models like the one we present can help to interpret experimental data on the synergy between radiotherapy and immunotherapy, and to assist in the design of more effective treatments. Index Terms—Radioimmunotherapy, radiotherapy, αPDL1, αCTLA4, biomathematical modeling. F 1 INTRODUCTION CANCER immunotherapy (IT) is a therapeutic strategy against cancer that aims at boosting and exploiting the natural immune response to control and cure tumors [1], [2]. Checkpoint inhibitors are a type of IT which is used for treatment of several cancers, including melanoma, prostate, NSCLC and leukemia [3], [4]. These inhibitors block different checkpoint proteins, like CTLA-4 and PD- 1/PD-L1, which are well known suppressors of the immune response against tumors. Inhibitors of these proteins have shown promising results in preclinical experiments, and there are several monoclonal antibodies against PD-1/PD- L1 approved to treat different types of cancer [5]. However, efficacy of IT as cancer treatment is still limited, but for particular cases. For example, ipilimumab has shown an improvement on the survival of melanoma patients, but response rates are low, in the 10-15% range [6]. •I. Gonz´alez-Crespo is with the Group of Medical Physics and Biomathematics, Instituto de Investigaci´on Sanitaria de Santiago, and the Department of Applied Mathematics, Universidade de Santiago de Compostela, Spain. •A. G´omez-Caama˜no is with the Department of Radiation Oncology, Clinical University Hospital of Santiago, Spain. •O. L´opez Pouso is with the Group of Medical Physics and Biomathematics, Instituto de Investigaci´on Sanitaria de Santiago, and the Department of Applied Mathematics, Universidade de Santiago de Compostela, Spain. •J.D. Fenwick is with the Department of Molecular and Clinical Cancer Medicine, Institute of Translational Medicine, University of Liverpool, Liverpool, United Kingdom. •J. Pardo-Montero is with the Group of Medical Physics and Biomathematics, Instituto de Investigaci´on Sanitaria de Santiago, Spain. E-mail: [email protected] Many preclinical studies have shown that the combination of radiotherapy (RT) and IT, in particular inhibitors of CTLA-4, and PD-1/PD-L1, is significantly more effective than RT and IT alone [7], [8], [9], [10], [11]. The dominant biological cell killing mechanism behind the effect of RT is the generation of double strand breaks in the DNA by ionizing particles [12]. However, there is evidence that radiation can trigger other cell killing mechanisms, particularly when delivered at high-doses per fraction, which may be important for the synergy with immunotherapy. High-doses of radiation can damage the tumor vascular system, eventually triggering cell death [13], and also lead to an immune response against surviving tumor cells [10], [14], [15], increasing the likelihood of tumor control. The mechanisms behind these induced immune effects seem to be related to the increased liberation of tumor antigens, which cause the activation of antitumor T-cells, and the modification of tumor microenvironment, killing immune down-regulators like Treg and MDSC cells, and facilitating T-cell infiltration in the tumor [7], [16], [17]. Despite the potential of IT and RT, how the best combination of both therapies can be achieved is still a matter of study. In this regard, validated biomathematical models of RT+IT (in silico tumor models) would be very useful. Biomathematical models, based on solid experimental and clinical data are of high importance in order to interpret results, and may also assist on the design of optimal therapeutic strategies, potentially guiding clinicians in the selection of optimal treatments. Modeling the response of tumors to IT has been addressed in the biomathematical literature, following both phenomenological and systems biology approaches [18], [19], [20], [21], [22], [23], [24]. On the other hand, modeling the synergistic combination of IT Authorized licensed use limited to: UNIVERSIDADE DE SANTIAGO. Downloaded on January 16,2023 at 09:26:31 UTC from IEEE Xplore. Restrictions apply. 1545-5963 (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TCBB.2022.3174454, IEEE/ACM Transactions on Computational Biology and Bioinformatics 3 model can be shown in a simple way as: dynamics C=proliferation −radiation death −immune death dynamics ˆ A=natural release +RT-mediated release −natural elimination −T-cell activation dynamics Ta=activation/infiltration −radiation death −immune death −natural elimination The model is built following a mechanistic approach rather than a purely phenomenological approach. It includes many simplifications, though, in order to present a tractable problem. Among the most relevant simplifications, we should cite: i) We do not include spatial coordinates in our multi-compartmental model; ii) The process of immune death is overly simplified. In particular, we do not include other cell types that participate in this process, either favoring immunity or acting as suppressors, like natural killers or T-regs; iii) We mostly rely on the linear-quadratic model (LQ) to account for radiation cell death, but departures from the LQ model at high-doses are studied and discussed; iv) The process of T-cell infiltration in the tumor is not included in the model. 2.2 Direct Cell Death and Kinetics Dose delivery is modeled as instantaneous (it seems like a good approximation, as in typical fractionated treatments dose delivery takes minutes, while the typical times of our model are days). Radiation tumor cell death is typically modeled with the linear-quadratic (LQ) model: log SF =−αd −βd2(1) where SF is the surviving fraction of a population of cells after being irradiated to a radiation dose, d, and αand βare the LQ linear and quadratic parameters. We will fit data with different fractionation schedules. It is well known that cell death can depart from the standard LQ-model, especially at high doses per fraction. Several effects contribute to this, including re-oxygenation [31], saturation [32], or vascular damage [13]. In order to investigate possible departures from the LQ model, we will check whether other models provide better fits. In particular, we will investigate a simple ad hoc modification of the quadratic term of the LQ-model [33]: β→β01 + c√d(2) where cis a free parameter and dis the dose per fraction. We will also use the Linear-Quadratic-Linear (LQL) model [34], which obtains the surviving fraction as: log SF =−αd −2βxd −1 + e−xd x2(3) where xis an extra parameter that modulates the slope of the curve. We will use (1), (2) and (3) to model radiation-induced tumor cell death. T-cell death will be described using the LQ-model. Cells fatally damaged by radiation (doomed) do not die instantly, but follow a given kinetics, generally a pause (mitotic delay) followed by a progressive death as lethally damaged cells enter mitosis and suffer mitotic catastrophe. We will model this by considering a mitotic delay followed by an exponential death [35]. Viable tumor cells can proliferate, and we describe this by using the logistic formalism [36]. On the other hand, while doomed cells may carry some proliferative capacity (abortive divisions [37]), it should be limited and does not contribute to the long-term cell population. Therefore, we will ignore it. From the above considerations we can write the following equation for viable tumor cells: dC dt (t) = λ1C(t)[1 −λ2Ctot(t)] −KC(t)(4) where λ1and λ2are constants and Ctot(t) = C(t) + Cd(t) is the total number of tumor cells at time t.KC(t)is an impulse term accounting for the effect of the radiation dose: KC(t) = (1 −SF C(d(t)))C(t)X i δ(t−ti)(5) where {ti}is the vector of radiation delivery times, {di} are the doses delivered at the times {ti}, and SF Cis the surviving fraction of tumor cells given by (1), (2) or (3). δ(x)is the Dirac delta function. We consider that each radiation fraction creates new doomed cells, but does not interfere with the radiation kinetics of existing ones. Therefore, we can split the compartment of doomed cells into ncompartments created by ndose fractions: dCd,i dt (t) = KC(t)−φω(¯ ti)Cd,i(t)(6) Cd(t) = X i Cd,i (7) Here, Cd,i(t)denotes doomed cells created by the radiation dose fraction di,tiis the delivery time of that fraction, and ¯ ti=t−ti. Notice that Cd,i(t)is defined as zero for t < ti. The parameter φis the death rate, and ωmodels the mitotic delay and progressive incorporation of damaged cells to cell death kinetics after a radiation fraction: ω(¯ t) =        0,for ¯ t≤τd1 ¯ t−τd1 τd2−τd1 ,for τd1<¯ t≤τd1 1,for ¯ t>τd2 (8) Radiation also kills T-cells present in the tumor as: dTa dt (t) = −KT(t)(9) where KT(t)is the impulse term for T-cells, which has the same form of (5), but with SF Tas the surviving fraction of T-cells. It is assumed that radiation-damaged T-cells die instantly, and so there are no kinetic terms associated to such process. 2.3 Antigen Release and T-cell Activation Antigens are considered to be released both naturally (rate proportional to the number of tumor cells, both viable and Authorized licensed use limited to: UNIVERSIDADE DE SANTIAGO. Downloaded on January 16,2023 at 09:26:31 UTC from IEEE Xplore. Restrictions apply. 1545-5963 (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TCBB.2022.3174454, IEEE/ACM Transactions on Computational Biology and Bioinformatics 4 doomed), and during radiation-induced cell death (proportional to the rate of cell death). We also include a term describing natural elimination. A biological delay, τ1, between antigen release and T-cell activation is included (which may be interpreted as the time that APCs take to collect antigens and carry them to the activation sites, Fig. 1(A)): dˆ A dt (t) = ρCtot(t−τ1) + ψφω(t−τ1)Cd(t−τ1) −σˆ A(t)(10) The activation of T-cells against tumor cells is modeled through four bilinear equations which describe the generation of activated T-cells ( ˆ Ta) or blocked T-cells ( ˆ Tb) (through the CTLA-4 receptor) from a pool of blank T-cells ( ˆ T): dˆ A dt (t) = −aˆ A(t)ˆ T(t)−bˆ A(t)ˆ T(t)(11) dˆ T dt (t) = −aˆ A(t)ˆ T(t)−bˆ A(t)ˆ T(t) + h(12) dˆ Ta dt (t) = aˆ A(t)ˆ T(t)(13) dˆ Tb dt (t) = bˆ A(t)ˆ T(t)(14) The constants aand b(r= 1 + b/a) describe the affinities for activation/inactivation, respectively [25]. The pool of blank T-cells starts from ˆ T(0) = T0, which is assumed to be the carrying capacity of T-cells. It can be depleted due to activation/inactivation, and in that situation it can renew at constant rate (due to maturation of new T-cells), h. A constraint is imposed to avoid the T-cell compartment from exceeding the carrying capacity: ˆ T(t)≤T0. Active T-cells, ˆ Ta, migrate and infiltrate in the tumor (with a biological delay τ2, which phenomenologically models the time needed by active T-cells to act on tumor cells) where they become part of the compartment Ta(note the hat notation): dTa dt (t) = dˆ Ta dt (t−τ2) = aˆ A(t−τ2)ˆ T(t−τ2)(15) We have tested the hypothesis that vascular damage at high radiation doses [13] may reduce the effectiveness of radioimmunotherapy by limiting the infiltration of T-cells in the tumor. Therefore, we include a dose and time dependent T-cell infiltrating parameter to account for vascular damage and recovery. Inspired by [13] and [38], we consider critical vascular damage for doses beyond 15 Gy, and a progressive recovery of vascular function as, f(t) = min{0.05t, 1}(16) where the time post-irradiation, t, is measured in days. This term represents the fraction of active T-cells reaching the tumor, and multiplies (15). 2.4 T-cell Mediated Tumor Cell Death Interaction between active T-cells and tumor cells results in the partial depletion of both. Following the work of de Pillis et al. [18], we model this interaction with a bilinear term for the compartment Ta, in addition to an exponential natural elimination: dTa dt (t) = −ιTa(t)Ctot(t)−ηTa(t)(17) On the other hand, T-cell mediated tumor cell death is modeled with the following term [18]: dC dt (t) = −p(Ta(t)/Ctot(t))q s+ (Ta(t)/Ctot(t))qC(t)(18) The same expression holds for Cd(t). Note that these terms are coupled to the earlier equations. 2.5 The Effect of αPDL1 and αCTLA4 The concentration biokinetics (in arbitrary units) of αPDL1 (p1) and αCTLA4 (c4) is modeled as an instantaneous source term at injection times ({tp1}and {tc4}, respectively) and a continuous exponential elimination: dc4 dt (t) = ic4(t)δ(t−{tc4})−νc4(t)(19) dp1 dt (t) = ip1(t)δ(t−{tp1})−µp1(t)(20) The precise pharmacokinetic modeling of these drugs is beyond the scope of the present article. However, this seems a good approximation for the kinetics of αCTLA4 as Selby et al. [39] investigated different αCTLA4 drugs biokinetics, finding that they follow linear forms like (19). Although Deng et al. [40] investigated the biokinetics of αPDL1, finding a more complex non-linear behaviour. Rather than considering a complex kinetics model for the characterization of the effect of αCTLA4 on the deinhibition of T-cells, we model it as a simple dependence on the parameter bin (11), (12) and (14) on c4, similarly to [25]: b→b 1 + c4(t)(21) On the other hand, the effect of αPDL1 on the immunedeath of tumor cells is modeled by introducing a dependence on the parameter pin (18) as: p→p(1 + p1(t)) (22) 2.6 The Complete Model The assembled model is the following for C,Cd,ˆ Aand Ta: dC dt (t) = λ1C(t) (1 −λ2Ctot(t)) −KC(t) −p(1 + p1(t)) (Ta(t)/Ctot(t))q s+ (Ta(t)/Ctot(t))qC(t) dCd,i dt (t) = KC(t)−φω(¯ ti)Cd,i(t) −p(1 + p1(t)) (Ta(t)/Ctot(t))q s+ (Ta(t)/Ctot(t))qCd,i(t) dTa dt (t) = −KT(t) + aˆ A(t−τ2)ˆ T(t−τ2)(23) −ιTa(t)Ctot(t)−ηTa(t) dˆ A dt (t) = ρCtot(t−τ1) + ψφ X i ω(¯ ti−τ1)Cd,i(t−τ1) −σˆ A(t)−aˆ A(t)ˆ T(t)−b 1 + c4(t)ˆ A(t)ˆ T(t) In addition, (12)-(14) control the activation of T-cells against tumor cells, and (19, 20) control the biokinetics of αPDL1 and αCTLA4. Authorized licensed use limited to: UNIVERSIDADE DE SANTIAGO. Downloaded on January 16,2023 at 09:26:31 UTC from IEEE Xplore. Restrictions apply. 1545-5963 (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TCBB.2022.3174454, IEEE/ACM Transactions on Computational Biology and Bioinformatics 5 2.7 Modeling Tumor Control Probability: Markov model We have employed the clonogenic cell hypothesis [41] to obtain tumor control probablities (TCP) from our model. It states that in order to control the tumor, all cells with proliferative capacity, which we identify with the compartment C, need to be eliminated. As defined in the previous section, the model is continuous and deterministic. In order to calculate TCPs we need a discrete model (numbers of cells) and stochasticity. Therefore, for low numbers of cells (C < 1000 cells), the model is converted to a Markov birth/death stochastic process [42] by interpreting terms in the differential equations as birth/death probabilities. In a simulation, the tumor is considered controlled if Creaches 0. In addition to the stochasticity of the Markov model, to obtain populational TCPs we also implemented random perturbations of the model parameters to simulate the heterogeneity of a population. For the population (a given number of simulations), TCP is computed as: TCP =number of controls number of simulations (24) More details about the TCP calculation and implementation are provided in the Supplementary Material. 2.8 Experimental Data and Model Fitting In [8] the authors studied the response of tumors in mice to RT (different fractionations) and αCTLA4, either as monotherapies or in combination. Tumor cells (TSA breast carcinoma cells) were planted on the side of mice and let grow for 12 days, when they reached a volume of ∼32 mm3. Treatments started at that time, and evolution of tumor volumes were monitored every 3 days. They studied tumor response to different combinations of RT+IT. In particular: i) no treatment; ii) RT alone, 20 Gy single-fraction; iii) RT alone, 3 fractions of 8 Gy (days 12, 13 and 14); iv) RT alone, 5 fractions of 6 Gy (days 12, 13, 14, 15, and 16); v) IT alone, delivered in 3 fractions (days 14, 17, 20); vi) combined RT+IT, 20 Gy + 3 fractions of IT (ii+v); vii) combined RT+IT, (6 Gy×5) + 3 fractions of IT (iv+v); viii) combined RT+IT, (8 Gy×3) + 3 fractions of IT (iii+v); ix) combined RT+IT, (8 Gy×3) + 3 fractions of IT at days 12, 15, 18; x) combined RT+IT, (8 Gy×3) + 3 fractions of IT at days 16, 18, 20. The study also reported the fraction of animals where tumor control was achieved (no evidence of tumor at the time of euthanasia). In [9] the authors presented responses of tumors to radiotherapy and αPDL1. Cancer cells (TUBO breast carcinoma cells) were implanted in mice, and tumors were allowed to grow for 14 days, when they had volumes around 120 mm3. At that time, treatments started and tumor volumes were monitored up to day 35. The different experimental arms were: i) no treatment; ii) single dose of 12 Gy; iii) four fractions of αPDL1 at days 14, 17, 20 and 23; iv) 12 Gy + αPDL1 (ii+iii). In order to fit the reported evolution of (populationaveraged) tumor volumes, we used the continuous model. Firstly, we let the modeled tumors to freely grow until they reached the relevant pre-treatment volumes reported in [8] and [9]. That time was defined as reference (day 0), and treatment times are defined relative to it. Notice that the time to reach those volumes differs from the experimental results, as we start with different numbers of cells than those experimentally injected and our model does not aim to describe the process of tumor growth, which may be dominated by different mechanisms than tumor response. Tumor volumes in our model are computed by considering the populations of both tumor cells (viable and doomed) and T-cells located in the tumor: Vmodel(t) = Ctot(t)VC+Ta(t)VT(25) where VCand VTare the volumes of individual tumor cells and T-cells, respectively. A simulated annealing method [43] was implemented to find best fitting parameters. The objective function to be minimized is the weighted sum of square differences between model and experimental values: F=X curves X points (Vmodel −Vexp)2 u2(26) where Vexp and Vmodel are the experimental and model results, and uare the experimental uncertainties. The optimization method has been applied to all response curves of each study at once (i.e. the sum above runs over different time points and different combinations of radiotherapy and immunotherapy), to avoid different best-fitting parameters for each curve. 2.9 Evaluation of goodness-of-fit We have used the Akaike Information Criterion (AIC) with sample size correction [44] to compare best fits obtained with different direct damage terms (1), (2) and (3). This methodology ranks models according to the likelihood of the fit, L, and number of free parameters of the model, k: AIC =−2 log (L)+2k+2k(k+ 1) N−k−1(27) where Nis the number of experimental data points, and Lis the maximum likelihood (of the best fit), calculated assuming a normal distribution of the experimental points and uncertainties. The model with the lowest AIC is considered the best model. 2.10 Biologically Effective Dose The biologically effective dose (BED) [45] was used to design different radiobiologically iso-effective fractionations. The BED of a schedule delivering a total dose Din fractions of dose dis given by: BED =D1 + d α/β (28) where α/β is the ratio of the linear and quadratic terms in the LQ model. 2.11 Qualitative Sensitivity Analysis We have performed a qualitative local parametric sensitivity analysis. In order to do so, we have evaluated the sensitivity of the cost function to parameter perturbations around bestfitting values as: Si=|F(x+∆xi)−F(x)|(29) Authorized licensed use limited to: UNIVERSIDADE DE SANTIAGO. Downloaded on January 16,2023 at 09:26:31 UTC from IEEE Xplore. Restrictions apply. 1545-5963 (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TCBB.2022.3174454, IEEE/ACM Transactions on Computational Biology and Bioinformatics 6 where Fis the cost function (26), xthe set of best-fitting parameters and ∆xi= (0, ..., 0,0.1xi,0, ..., 0) (setting a 10% perturbation with respect to the best-fitting value xi). 2.12 Implementation and Parameters The model was implemented using different functions in Matlab (The Matworks, Natick, MA). The model is solved by employing an explicit Euler method [46], with a time step of 0.05 days (details about the behavior of the method are available in the Supplementary Material). The main functions and the data used for model fitting are available from the Dataverse repository [47]. Not all model parameters were free during the fit to experimental data. Cell volumes in (25) were set to VC= 10−6mm3and VT= 2 ×10−7mm3[48]. For the radiosensitivity of T-cells we have fixed the LQ parameters to αT= 0.1and βT=αT/10. The biokinetic elimination rates of αPDL1 and αCTLA4 were set to µ= 0.5days−1 and ν= 0.1days−1respectively, from fits to data reported in [39], [40]. The parameters characterizing the mitotic delay of radiation-damaged cells were set to τd1= 1 days, τd2= 1.5days, which is in the range of reported mitotic delays. The parameter r, relative to the activation/inactivation rate of T-cells is set to 5, as in [25]. On the other hand, values of best-fitting parameters were constrained to qualitative reasonable bounding intervals when deemed necessary, in order to avoid unphysical/unreasonable values (for example, α > 0 and β > 0 in the LQ-model, or λ > 0 for tumor proliferation). The bounding intervals are shown in Supplementary Table 4. 3 RESULTS 3.1 The Model Can Fit Pre-clinical Data of Tumor Response to Combined Therapies of Radiation and αPDL1/αCTLA4 In Fig. 2 we report best fits of our model to volume dynamics data presented in [8], when employing the modified LQ- model for direct radiation death (2). In this dataset, when no treatment is delivered, tumor volumes grow exponentially, αCTLA4 alone has no significant effect on tumor response, RT alone causes a moderate tumor response, and the combination of RT and αCTLA4 leads to an important tumor response, achieving tumor control in some cases (mostly for 8 Gy×3). Our model reproduce these progression patterns, as shown in Fig. 2. Best-fitting parameters are presented in Supplementary Table 1. We highlight the most relevant parameters associated with proliferation (λ1≃0.14 day−1), radiation damage (αC≃0.02 Gy−1,βC≃0.007 Gy−2, c≃ −0.2 Gy−1/2) and the immune effect on tumor cells (p≃24.4 day−1) and T-cells (ι≃2×10−8day−1). Data fitting point to a decrease in relative radiosensitivity with increasing dose (negative parameter cin (2)). Such behavior can be described by the LQL model. Therefore, we have performed the same fit with the LQL model instead (Fig. 3), obtaining similar results. Best-fitting parameters are reported in Supplementary Table 2, including the most relevant parameters associated with proliferation (λ1≃0.14 day−1), radiation damage (αC≃0.04 Gy−1, βC≃0.017 Gy−2,x≃8.4 Gy−1) and the immune effect on tumor cells (p≃23.3 day−1) and T-cells (ι≃ 2×10−8day−1). The best-fitting value of the cost function is F= 29.56 and F= 42.18 for the modified LQ and LQL model, respectively. In Fig. 4 we report best fits of our model to data presented in [9], which shows tumor responses to radiotherapy and αPDL1. This dataset presents similar patterns of tumor response: αPDL1 alone has not significant effect on tumor response, RT alone causes a moderate tumor response, and the combination of RT+αPDL1 presents synergy and leads to an important tumor response. To fit the data of Fig. 4 we have kept fixed most of the best-fitting parameters obtained when fitting Fig. 2: only parameters related to the dose of αPDL1, tumor cell proliferation (λ1≃0.12 day−1) and tumor cell radiosensitivity (αC≃0.03 Gy−1) were allowed to vary. Because this dataset only includes one dose per fraction, we have used the LQ model with αC/βC= 10 Gy (i.e. only αCis a free parameter). While there are differences in the clones and tumors that could justify using different host-related and tumor-related parameters, we think that imposing such constraints on the optimization poses a serious test to our model, and avoids reaching good fits by over-fitting. Bestfitting parameters are reported in Supplementary Table 1. 3.2 Vascular Damage May Limit the Effectiveness of Radioimmunotherapy Large radiation doses can seriously damage tumor vasculature, which might limit the infiltration of active T-cells in the tumor. This might also explain the poorer results obtained with the 20 Gy single-fraction irradiation in [8]. In order to test the hypothesis that vascular damage may affect the effectiveness of radioimmunotherapy, we removed the dependence of the tumor cells β-term on the radiation dose (see Section 2.2), and we included a dose and time dependent T-cell infiltrating parameter in our model (16), to account for vascular damage and recovery. Inspired by [13], [38], we consider critical vascular damage for irradiation above 15 Gy, followed by a progressive recovery of vascular function. This factor represents the fraction of active T-cells reaching the tumor. In Fig. 5 we show best fits of our model to data presented in [8]. The model provides a good fit to the experimental values. The goodness of the fit is slightly better than those reported in Fig. 2 and Fig. 3: AIC = 436.24 (k= 19, F= 25.33,N= 50), versus AIC = 446.02 (k= 20, F= 29.59) and AIC = 458.75 (k= 20,F= 42.18) obtained with the modified LQ and LQL models, respectively. Best-fitting model parameters are shown in Supplementary Table 3, including the most relevant parameters associated with proliferation (λ1≃0.14 day−1), radiation damage (αC≃0.02 Gy−1,βC≃0.002 Gy−2) and the immune effect on tumor cells (p≃24.9 day−1) and T-cells (ι≃6×10−9day−1). In Table 1 we rank the sensitivity of the cost function to model parameters. The model is most sensitive to parameters describing tumor cell proliferation, radiosensitivity and immune-mediated tumor cell killing. Volumes presented in previous figures include tumor cells and T-cells. Certainly, we do not want the model to reproduce tumor volumes by including low fractions of tumor Authorized licensed use limited to: UNIVERSIDADE DE SANTIAGO. Downloaded on January 16,2023 at 09:26:31 UTC from IEEE Xplore. Restrictions apply.