Chemoradiotherapy of locally-advanced non-small cell lung cancer: analysis of radiation dose-response, chemotherapy and survival-limiting toxicity effects

Summary: 2-year overall survival in 68 trial arms was best described by a model that accounted for classical radiobiological factors, chemotherapy effects and survival-limiting toxicity. The fitted  /  ratio was 4.0 Gy and repopulation negated 0.38 Gy/day. Modelled survival peaked at 80 (stage IIIA) and 87 Gy (IIIB) for radiotherapy and sequential chemoradiotherapy delivered in 2 Gy fractions over 40 days, and 67 (IIIA) and 73 Gy (IIIB) for concurrent chemoradiotherapy, before falling at higher doses.

according to the model fit, if critical normal structures in which survival-limiting toxicities arise can be identified and selectively spared.

Introduction
Overall survival (OS) following radiotherapy (RT) for locally-advanced non-small cell lung cancer (LA-NSCLC) remains disappointing.While some dose-escalation trials have achieved promising results, [1] outcomes have been inconsistent, with median OS in better performing studies being 25-28 months.[2,3] Surprisingly, the RTOG-0617 phase III study of concurrent chemoradiotherapy (CRT) found a sub-unity median survival ratio (MSR) of 0.71 for 74 Gy in 37 daily fractions versus 60 Gy in 30 fractions.[2] Although hyperfractionated and moderately hypofractionated schedules have been trialled, [1,[3][4][5][6] the standard-of-care remains 60-66 Gy in 30-33 fractions with some UK centres preferring 55 Gy in 20 fractions.Improved survival has recently been demonstrated using radioimmunotherapy [7] but the optimal RT schedule has yet to be determined.
A need therefore exists to reconcile apparently inconsistent trial results, to improve outcome prediction for modified RT of NSCLC.Partridge et al analyzed 2-year diseasefree survival (DFS), demonstrating a dose-response using a probit tumour control probability (TCP) model.[8] Differences in dose-per-fraction and schedule duration were accounted for by standardizing to equivalent doses in 2 Gy fractions delivered over a fixed treatment duration (EQD2T), using the time-corrected linear-quadratic formalism.Values of /,  and TK parameters describing fractionation sensitivity, accelerated tumor proliferation and onset delay are not well established for NSCLC, [9][10][11][12], and therefore Partridge et al chose 10 Gy, 0.6 Gy/day and 21 days, similar to values obtained from analyses of outcomes for head-and-neck squamous cell carcinoma (HNSCC).[13] A recent meta-analysis of randomized studies reported an MSR of 1.13 (p=0.002) for higher radiation dose arms in RT-only and sequential CRT trials, but 0.83 (p = 0.02) for higher dose arms in concurrent CRT trials [6], demonstrating the need for dose-response models to describe differential chemotherapy effects and possible reductions in survival at high radiation doses.Here, we analyze 2-year OS rates (OS2yr) compiled for ~7000 NSCLC patients from published phase I-III trials of RT alone and CRT.Starting from the Partridge model, we freely fit /,  and Tk before parsimoniously extending the model to describe sequential and concurrent chemotherapy effects, and survival reductions at high doses due to radiation toxicity.[14] 2. Materials and methods

Data
Phase I-III trial data were identified from PubMed, ScienceDirect and Google Scholar searches for the MeSH term 'NSCLC radiotherapy dose-escalation'.Citation-following yielded further data.The dataset was limited to studies published after 1995 describing outcomes for >20 patients per trial arm.Reports missing unambiguous details of dose, fractionation, treatment duration, chemotherapy scheduling, stage-mix and OS2yr data were excluded.Strongly hypofractionated schedules used to treat predominantly earlystage disease (dose-per-fraction >4.0 Gy) were also excluded, substantially weighting the dataset towards higher-stage disease.
Reported prescribed doses and doses-per-fraction were increased by 5% for North American trials not employing lung tissue heterogeneity corrections.[15] Doses prescribed in the RTOG-0617 study were also raised by 5% despite heterogeneity corrections having been applied, since prescription was to 95% of the planning target volume, generating physical isocentre doses similar to those of isocentrically-prescribed treatments planned without heterogeneity corrections.[16]

Fitted dose-response models
Three models have been fitted to OS2yr data.The first accounts for dose, dose-per-fraction and schedule duration.The second contains two additional parameters describing systemic and local chemotherapy effects.The third contains two further parameters describing survival reductions at high doses, and one accounting for improvements in survival over time.

Model 1: Standard probit-EQD2
Model 1, based on that of Partridge et al, [8] describes TCP varying sigmoidally with EQD2T delivered to the tumor, EQD2Ttum where  is the cumulative normal distribution, EQD2T tum,50 is the tumor EQD2T required for 50% control, and m defines the dose-response gradient.EQD2T tum was calculated as D and d being prescribed dose and dose-per-fraction, and T treatment duration.[8] Following Partridge's approach, cohort-specific OS was calculated as a weighted sum of modelled TCPs for each disease stage where fi is the fraction of patients in the cohort with stage i (I, II, IIIA or IIIB) disease, and EQD2T tum,50 (  ) the EQD2Ttum,50 for that stage.Initially, /tum,  and Tk values were fixed at the levels chosen by Partridge, and the five parameters m and EQD2T tum,50 ( I, II, IIIA, IIIB ) were fitted to achieve the best description of the data ('Partridge model').Subsequently all eight model parameters were fitted ('Model 1').

Model 2: Adding chemotherapy effects
Survival is improved by chemotherapy-driven reductions in the distant failure-rate.[17,18].
Model 2 describes this effect via an additional factor where OS2yr-max values <100% account for patients who died of distant metastases or unrelated causes despite achieving loco-regional control.For RT alone, OS 2yr-max RT-only was fixed at 85%, based on a reported distant failure-rate of 15% post-surgery for patients who did not receive chemotherapy.[19] For CRT treatments, OS 2yr-max CRT was fitted to obtain the best match of the model to the data, a value >85% describing reduced distant failurerelated mortality due to systemic effects of chemotherapy.
Survival is longer following concurrent than sequential CRT.[20,21] This is accounted for by scaling the EQD2Ttum doses of concurrent CRT treatments by an additional fitted parameter RS cCRT >1 to describe possible tumor radiosensitization by concurrent chemotherapy Model 2 has ten fitted parameters.Its residuals revealed an unfitted trend for OS2yr to rise with study publication year (see Supplementary Information), presumably reflecting advances in treatment and staging over time unrelated to CRT scheduling, and concurring with a reported trend for significantly longer survival in more recent studies [6].Survival was therefore adjusted further via Y being the number of years before 2016 a study was published, and R a fitted parameter describing survival benefit with time.
2.2.2.Model 3: Reduced survival at higher doses Model 2 cannot describe falling survival at high doses.However, the concept of a therapeutic window implies that OS must rise, plateau and eventually fall as RT dose increases.To accommodate this, we extended the model to where the modelled survival-limiting toxicity-rate, SLT, is given by and increases sigmoidally with the normal tissue EQD2 calculated using / = 3 Gy.The fitted quantities mNT and EQD250,NT define the toxicity dose-response gradient and the EQD2NT at which survival is halved by toxicity, and bring the number of fitted parameters in Model 3 to 13.

Statistical methods
Models were fitted to OS2yr data using the maximum-likelihood method, via the mle2 package within the 'R' language (v3.4.0).[22] Significances of improvements in model-fit were determined using likelihood-ratio testing.[23] Relative model performance in future data was estimated using the Akaike Information Criterion (AIC), which trades off goodness-of-fit against a penalty term that rises with the number of fitted parameters to account for possible overfitting of the dataset.[24] Leave-one-out cross-validation was performed to check the AIC findings, calculating weighted sums-of-squared residuals.
For each model, asymptotic confidence intervals (CIs) were calculated on all fitted parameters.For the model judged best on AIC and cross-validation, profile-likelihood CIs were also determined.[22] 3. Results
Figure 1a shows observed OS2yr plotted against prescribed physical dose.Although the data appear highly dispersed, local regression (LOESS, smoothing = 0.7) shows survival increasing with dose before plateauing at ~80 Gy.Physical dose and dose-per fraction are plotted against RT duration in Figure 1b and Supplementary Figure S1.Higher doses were typically given using longer schedules (Spearman's rho = 0.62, p <0.001), but dose-perfraction and schedule duration were not significantly correlated (p = 0.66).

Data fits
Model fits to 2-year survival data are detailed in Table 2.The data were described significantly better when classical radiobiological parameters /,  and TK were fitted rather than fixed at the levels of the Partridge model (p<10 -5 ).The fit was further improved by additional terms describing chemotherapy effects (p<10 - 19 ), longer survival in later studies (p<10 -8 , Supplementary Table S1 and Figure S2), and reduced survival at high dose-levels (p<10 -6 ).Model 3 included all these factors and described the data best (Figure 2).
Values of OS2yr predicted by Model 3 are plotted in Figure 3a and rise with observed OS2yr, with notably less data dispersion than in Figure 1.A plot of predicted and observed OS2yr versus EQD2Ttum shows no trend in residuals (Figure 3b).On AIC and cross-validation measures Model 3 substantially out-performed the other models (Table 2), indicating that its success was not due to over-fitting of this specific dataset.
For Model 2b a higher / value of 7.4 Gy (95%CI: 2.7, 12.4 Gy) was obtained, but Model 3 described the dataset substantially better.In fits of Models 1 and 2, tumor repopulation ran faster at 0.71 and 0.97 Gy/day starting at day 28, but again Model 3 described the data significantly better.The Model 3 fit described a maximum dose-response gradient of 1.25% (95% CI: 1.13, 1.42%) gain in OS2yr per 1% increase in EQD2Ttum, in the absence of survival-limiting toxicity.

Dose-response curves
Dose-response curves for OS2yr described by Model 3 are plotted in Figure 4 for RT alone, sequential and concurrent CRT.For ease of comparison with a reference schedule delivering 60 Gy in 30 fractions over 6 weeks (40 days), OS2yr is plotted against dose delivered in 2 Gy fractions over 40 days, EQD2(40d)tum, where EQD2(40d)tum = EQD2Ttum +  (40-TK) = EQD2Ttum + 10.6 Gy for Model 3 (10) The plots show predicted OS2yr peaking at 52% and 38% respectively for stages IIIA and IIIB NSCLC treated using RT alone, and at 59% and 42% for CRT treatments.For RT alone and sequential CRT, peak survival is reached at EQD2(40d)tum doses of 80 and 87 Gy for stages IIIA and IIIB respectively.For concurrent CRT, peak survival requires 67 and 73 Gy for IIIA and IIIB disease.Some of these dose-levels, 87 Gy for example, would exceed tolerance for serious toxicities [8].At higher dose-levels the model fit describes falling survival, gains in tumor control being more than offset by survival-limiting toxicity.

Discussion
NSCLC OS2yr data was described best by Model 3 (Table 2, Figure 2), which includes terms accounting for sequential and concurrent CRT, and toxicity-related reductions in survival at high doses.Although survival is known to be improved by CRT, [20] and is suspected to be limited by toxicity at high doses, [14,65] to our knowledge this is the first time these factors have been included in comprehensive dose-response modelling for NSCLC.
The 4.0 Gy (95%CI: 3.0, 5.4 Gy) / ratio of the Model 3 fit is notably lower than the 10 Gy value often assumed for NSCLC [10], and allowed the model to describe the data significantly better (p <.002, likelihood-ratio).A similar value of 4.9 Gy (95%CI: 3.0, 6.8 Gy) was obtained when the radiobiological parameters of the standard Partridge model were freely-fitted (Model 1, Table 2), indicating that the result is robust and not an artefact of modelling CRT and toxicity-limited survival.We are unaware of such a low / ratio being previously reported for LA-NSCLC, but the 4.0 Gy / value is consistent with an observation of Partridge et al. [8] that hypofractionated schedules appeared to overperform in plots of DFS versus EQD2Ttum calculated for / = 10 Gy, and concurs with / estimates of 3.9 and 2.8 Gy obtained from fits to stage I NSCLC data [66,67].The fitted 0.38 Gy/day loss of EQD2 to repopulation is similar to an early estimate of 0.45 Gy/day made for NSCLC by Koukourakis et al [68], but is lower than estimates for HNSCC (~0.6 Gy/day starting 3-5 weeks into RT [13]), and depends strongly on the introduction of the toxicity term into Model 3. Fits of the simpler but less well-performing Models 1, 2 and 2b described higher loss-rates of 0.6-1.0Gy/day typically starting at day 28 of RT, a consequence of these models using rapid repopulation rather than survival-limiting toxicity to fit the observed plateau in OS2yr at high doses, albeit less well.
The systemic effect of chemotherapy compared to RT alone was described in the Model 3 fit by a 10% rise in OS2yr-max, the overall survival achievable by a treatment if 100% of primary tumours were cured and no survival-limiting toxicity occurred.And the local effect of concurrent chemotherapy was described by a modelled 23% increase in effective tumor EQD2, greater than estimated previously [9] but consistent with the 0.86 hazard ratio reported by O'Rourke et al for concurrent rather than sequential CRT [20].
The model fit was significantly improved by the toxicity term, which was introduced to account for survival reductions at high dose-levels, and described survival-limiting toxicity rates of 4%, 11% and 22% at 50, 66 and 80 Gy EQD2NT respectively following RT alone or sequential CRT, and 8%, 23% and 42% following concurrent CRT.These modelled rates suggest that survival in dose-escalation trials may have been reduced by toxicity, and are in line with data from a cause-of-death analysis [69], which reported an ~15% rate of nonmalignant cardiac mortality at 2 years post-treatment for a curatively treated stage III cohort.Research is being carried out to identify mediastinal substructures in which survival-limiting toxicities arise, with a focus on heart [14], offering the prospect of improved outcomes if doses to these structures can be limited without compromising tumor coverage.
Whilst Model 3 was designed to parsimoniously account for several classical and novel dose-response factors, it nevertheless requires fitting of 13 parameters.This relative complexity is justified by the model's superior AIC and cross-validation scores, compared to those of simpler models.Further terms might provide a more complete description of the dataset, or might simply overfit it.We have investigated stage-specific survival-limiting toxicity terms, motivated by the observation that critical structures are less likely to receive high doses during RT of early-stage disease: resulting improvements in AIC score were marginal (Supplementary Information).
In summary, we have fitted a dataset of OS2yr rates reported for many RT and CRT schedules using a model that combined a classical description of radiation dose-response with novel terms accounting for chemotherapy effects and survival reductions at high doses.The model fit had an / ratio of 4.0 Gy and described a rate of EQD2 loss due to tumor repopulation of 0.38 Gy/day, implying that moderate acceleration achieved via hypofractionation would produce useful survival gains.
The fit predicts maximal OS2yr rates of 52% (stage IIIA) and 38% (IIIB) for RT alone when given in 2 Gy fractions over 6 weeks, and 59% (IIIA) and 42% (IIIB) for CRT, these peak rates being achieved at prescribed doses of 80 Gy (IIIA) and 87 Gy (IIIB) for RT alone and sequential CRT, and 67 Gy (IIIA) and 73 Gy (IIIB) for concurrent CRT.According to the fit, 10-20% further improvements in OS2yr might be achieved if normal tissues in which survival-limiting toxicities arise could be identified and spared without compromising tumor dose coverage.Table S2.Effect of fitting the fraction f(SI) of stage I patients in whom critical normal structures receive the prescribed dose-level and survival-limiting toxicities can occur (set to 100% for all stages in Model 3, and fitted for stage I only in the modified model).
Tumor repopulation negates 0.38 Gy EQD2 per day of treatment extension beyond day 12 of RT according to the Model 3 fit, making RT acceleration worth ~3 Gy per week of schedule shortening.Since the fitted tumour / value of 4.0 Gy lies close to the generic late toxicity / ratio of 3 Gy, such acceleration might be achieved efficiently via moderate hypofractionation (dose-per-fraction >2 Gy) including simultaneous boost techniques.For example, 56 Gy in 20 fractions over 26 days offers a modelled 5% advantage in OS2yr for stage IIIB NSCLC compared to 63 Gy in 30 fractions over 40 days, the two schedules being equivalent for late normal tissue damage (/ = 3 Gy).

Figure 1 .
Figure 1.a) Observed OS2yr vs. prescribed physical dose for the analyzed trial arms.LOESS regression (solid line) indicates dose-response in the region 40-80 Gy. b) Physical dose vs. RT schedule duration (RT only -open, sequential chemo-RT -striped, concurrent CRT -closed; larger symbols indicate better OS2yr).

Figure 2 .
Figure 2. Calibration plots of predicted versus observed OS2yr across the dataset, for fits of Models 1-3.Patients per trial arm are represented by areas of plotted points.Weighted least-square fits to the data (solid lines) are shown along with the line-of-identity.The gradients of the plots increase with improving fit quality.

Figure 4 .
Figure 4. (a-c) OS2yr dose-response curves calculated from the Model 3 fit for (a) RT alone, (b) sequential CRT, (c) concurrent CRT.Curves for stage I, II, IIIA and IIIB NSCLC are plotted as dashed grey, dashed black, solid grey and solid black lines.The peaks of the IIIA curves are picked out.(d) Schematic showing composition of OS as a product of survival unlimited by toxicity (OSno-tox) and (100% -survival-limiting toxicity percentage rate).

Figure S1 .
Figure S1.Dose-per-fraction plotted against RT duration for the trial arms analyzed.Patient numbers-per-trial arm are represented by areas of plotted points.

Figure S2 :
Figure S2: Model 2 fit residuals (differences between observed and fitted OS2yr) versus year of publication.Symbol sizes reflect patient numbers per trial arm.The solid line is a weighted least squares linear fit to the data.
*1174 patients were split 50/50% stage I versus II; the split for a further 37 patients was not published and 50/50% was assumed.*