Challenges in radiobiological modeling: can we decide between LQ and LQ-L models based on reviewed clinical NSCLC treatment outcome data?

Aim To study the dose-response of stage I non-small-cell lung cancer (NSCLC) in terms of long-term local tumor control (LC) after conventional and hypofractionated photon radiotherapy, modeled with the linear-quadratic (LQ) and linear-quadratic-linear (LQ-L) approaches and to estimate the clinical α/β ratio within the LQ frame. Material and methods We identified studies of curative radiotherapy as single treatment through MedLine search reporting 3-year LC as primary outcome of interest. Logistic models coupled with the biologically effective dose (BED) at isocenter and PTV edge according to both the LQ and LQ-L models with α/β = 10 Gy were fitted. Additionally, α/β was estimated from direct LQ fits. Results Thirty one studies were included reporting outcome of 2319 patients. The LQ-L fit yielded a significant value of 11.0 ± 5.2 Gy for the dose threshold (Dt) for BED10 at the isocenter. The LQ and LQ-L fits did not differ substantially. Concerning the estimation of α/β, the value obtained from the direct LQ fit for the complete fractionation range was 3.9 [68 % CI: 2.2–9.0] Gy (p > 0.05). Conclusion Both LQ and LQ-L fits can model local tumor control after conventionally and hypofractionated irradiation and are robust methods for predicting clinical effects. The observed dose-effect for local control in NSCLC is weaker at high doses due to data dispersion. For BED10 values of 100–150 Gy in ≥3 fractions, the differences in isoeffects predicted by both models can be neglected. Electronic supplementary material The online version of this article (doi:10.1186/s13014-016-0643-5) contains supplementary material, which is available to authorized users.


Introduction
The linear-quadratic (LQ) model was developed to describe experimental survival curves of both normal and tumor cells after irradiation. The LQ model fits the cell surviving fraction through a second-order polynomial on the dose per fraction, with coefficients α and β. The ratio between both coefficients describes the repair capacity of the cells and thus sensitivity to fractionation [1,2].
The LQ model provides an accurate description of fractionation effects at doses between 1 and 8-10 Gy per fraction [3]. Essentially, this formalism enables isoeffect calculations in current clinical practice, defining the relationships between the biological irradiation effect and key parameters such as dose per fraction, total number of fractions and treatment time. Advancements in this model led to the two most extended, complementary approaches for isoeffect calculation: the biologically effective dose (BED) and the equivalent dose in 2 Gy per fraction (EQD2) [4].
Current treatment of choice for stage I non-small cell lung cancer (NSCLC) is surgical tumor extraction. Since it became technically feasible, radiotherapy has been used as an alternative treatment method in inoperable cases. Early approaches used 3D-conformal conventionally fractionated techniques, whereas today, stereotactic body radiotherapy (SBRT) allows for highly precise delivery of radiation, thus enabling hypofractionation to deliver ablative radiation doses in 1-5 fractions [5]. Therefore, SBRT evolved to be the current treatment of choice for early-stage NSCLC in medically inoperable patients and in patients who do not consent to surgery. The high precision of dose delivery facilitates normal tissue sparing, even allowing for dose escalation to potentially improve local control. Despite a growing pool of clinical outcome data, the optimal total dose and fractionation scheme to reach the intended biological effects in terms of both local tumor control and side effects are still under debate [6].
Hypofractionation requires reliable isoeffect calculations. Thus, the relevance of the question has been renewed, whether the description of radiobiological effects based on the LQ formalism is appropriate for hypofractionated treatments [7][8][9].
If no α/β ratio estimation is available for a specific tumor entity, a generic value of 10 Gy is used for BED calculations, although the precision achievable with such a standard α/β ratio is assumed to be lower. Numerous attempts to calculate clinical α/β values have been made, using available clinical outcome data [10]. Such estimations are specially needed if the α/β of a specific tumor entity is suspected to be lower than 10 Gy. In this case a modified fractionation scheme, which reduces tumor cell recovery between fractions, could increase the therapeutic ratio as is the case e.g. in prostate carcinoma or breast cancer.
Many recent studies aim at outcome review and modeling of the dose response relationship of NSCLC [11][12][13][14][15]. However, studies attempting to estimate the α/ β ratio for NSCLC are scarce [11,15,16]. It is subject of current debate if the improved outcomes of hypofractionated SBRT are a consequence of an α/β ratio lower than 10 Gy, or even lower than the α/β value of the surrounding normal tissue, which could add a radiobiological rationale to the use of hypofractionation. Alternatively, the improvement could be caused by a reduced repopulation in a shorter overall treatment time.
In addition, in particle radiotherapy, a currently emerging field in radiation oncology, radiobiological considerations are of importance. For proton radiotherapy hypofractionated concepts are aimed for partially as motion management strategy [17][18][19][20], so that isoeffect calculations are essential. Apart from isoeffect calculations current treatment planning strategies for light ion therapy also require the attribution of radiobiological properties to both tumor and normal tissues. Specifically, in scanned-beam carbon ion therapy, radiosensitivity is characterized through α/β values obtained from photon irradiation experiments in vitro [21] in one of the mathematical models describing the enhanced biological effect in the Bragg peak, the so called local effect model (LEM), which is implemented in commercial treatment planning systems.
In such situation, clinical long time follow up data is the most valid source of data for modeling approaches, which is however inherently limited by inhomogeneity of treatment parameters and treatment techniques evolving over time. We focused our analysis on the comparison of the LQ-L versus the LQ model, since most of the mathematical model corrections to the LQ model proposed need additional input parameters [16,[22][23][24], which are not available for the specific clinical situation.
Therefore, this work aims at: 1) investigating the dose-response of NSCLC tumor control data from conventionally fractionated (CF) and stereotactic, hypofractionated radiotherapy treatments (HF), based on a review of published long-term outcome results, 2) evaluating the validity of the LQ and LQ-L models for both conventional and SBRT treatments, 3) and obtaining an estimation of the clinical α/β ratio of NSCLC.

Study design
We identified inclusion criteria, search strategy, outcome measures of interest and indispensable treatment parameters for the study. The analysis keeps standards of the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA) statements [25].

Selection criteria
PubMed was searched (July 2015) without restrictions on the publication date. Abstracts of conference proceedings were excluded and language was restricted to English. Repeated publications based on the same cohort were excluded as were outcome reviews in order to avoid duplicity of cohorts. The selection of studies based on the following criteria was made by two independent researchers. Study cohorts were eligible only if the following criteria were fulfilled: 1) Patients with stage I NSCLC (cT1/2, cN0, cM0) with either central or peripheral tumor location. 2) Treatment with photon radiotherapy with curative intent either 2D or 3D-conformal radiotherapy or SBRT as single modality treatment. Treatment could be delivered with CyberKnife, GammaKnife or linacbased without restrictions on fractionation schemes (e.g. normofractionated, accelerated, hyper-or hypofractionated schedules, single fraction irradiation), provided complete information on number of fractions, dose per fraction, total absorbed dose and overall treatment time was available. 3) Reported outcome of interest (actuarial 3-year local control estimations according to Kaplan-Meier or other methods, in which death was a censored event) with a median follow-up time in each cohort ≥17 months. 4) From an original study (i.a. prospective randomized controlled trial or prospective or retrospective observational study or case series). 5) The reported cohorts must include at least 25 patients. 6) Patients could be judged to be either operable or inoperable if the clinical decision was made in favor of radiotherapy. 7) A small amount of heterogeneity in reported dose and patient parameters for a given outcome data point was accepted (e.g. a marginal group of patients treated with a deviating radiation scheme, or a small proportion of tumors included in the report of a stage I population with clinical staging other than I).

Data extraction
The main endpoint of this review was local tumor control (LC) at 3 years. If this value was not explicitly provided in the text, it was extracted from the Kaplan-Meier diagrams. The fractionation concept (dose per fraction, number of fractions and total treatment time), and the planning technique were extracted together with further treatment-and patient-related parameters. The range of variation in reported cohort parameters for each local tumor control data point was qualitatively assessed and highly heterogeneous cohorts were excluded. Reported mean or preferably median dose values were used to describe the outcome of each specific patient cohort. All prescription doses were translated to doses at the isocenter and at PTV edge, calculated according to the information provided in each publication (more details in the Additional file 1). Mathematical modeling was performed with both, doses at the isocenter and at the PTV edge. The treatments were classified to be hypofractionated (HF) if 1-10 fractions were delivered with doses per fraction at the isocenter above 6 Gy. Treatments were classified as conventionally fractionated (CF) with a broader definition than in the clinical convention, based on the validity limits of the LQ model according to our current knowledge, namely treatment delivered in more than 10 fractions with fraction doses at isocenter ranging between 1.2 and 6 Gy.

Data analysis and mathematical models
Model parameters were fitted with nonlinear least square optimization methods and confidence intervals were calculated with likelihood profiling. A logistic relationship between tumor control probability (TCP) and the biological effective dose (BED) was assumed, according to the parameterization described in Okunieff et al. [26] and Bentzen et al. [27]. BED was based on the LQ model, calculated from the number of fractions and the dose at the isocenter, taking into account neither repopulation nor hypoxia, according to Eq. 1: where TCD 50 is the dose necessary to obtain a local tumor control of 50 % and k is a parameter with dose units that is used to calculate the normalized slope, γ 50 . This parameter quantifies the change in the expected TCP when a 1 % change in dose occurs, evaluated at the dose level of the TCD 50 , and represents the maximal slope of the dose-response relationship. It can be calculated from k and TCD 50 with the expression [27,28]: The same logistic model was implemented with an alternative BED definition, including a transition from linear-quadratic dependence for the cell survival to purely linear beyond a certain dose level, the dose threshold D t , as described in [22]: where n is the number of irradiation fractions, d is the fraction dose, and D t is the threshold for the fraction dose.
The LQ and LQ-L models were fitted to the joint dataset, and to the CF and HF subsets separately, for BED doses calculated both at the isocenter and the PTV edge. First, the linear-quadratic (LQ) model was applied with α/β fixed to 10 Gy, as it is universally accepted for conventional fractionation. The alternative BED definition derived from the LQ-L model was also applied with α/β equal to 10 Gy, to test if the inclusion of a dose threshold D t would improve the previous fit. Additionally, a study to tentatively estimate the α/β ratio from these clinical data was carried out.
In summary, the following fits were calculated: 1) LQ model with α/β ratio fixed to 10 Gy on the full dataset, and on the HF and CF datasets separately. 2) LQ-L model with α/β ratio fixed to 10 Gy on the full and on the HF datasets. 3) LQ model with free γ 50 , TCD 50 and α/β ratio on the full and the CF datasets.
Finally, in order to compare with the results of the recent analysis of Chi et al. [15] the Spearman's correlation between 3-year LC and BED calculated with different fixed α/β ratios was also investigated on the full and on the HF datasets.
The quality of the fits was assessed with checks of the residuals for normality. Models were ranked according to the Akaike Information Criterion (AIC), and for the LQ and LQ-L, being nested models, maximum likelihood ratio tests were made. All fitted and calculated values are reported together with their 68 % confidence intervals (CI), whenever possible. Statistical significance was assumed for p values < 0.05. Data handling, statistical analysis, model fitting and graphing were done with the software package R, version 2.15.0 [29].

Selected patient cohorts and description of the studies
In total, 31 studies were identified, which fulfilled the selection criteria, Of those, 8 studies report outcomes after conventionally fractionated treatments of a total of 344 patients [30][31][32][33][34][35][36][37] and 23 studies including 1975 patients reporting on hypofractionated irradiations . A total of 34 local control -schedule data points, with doses per fraction ranging from 30 to 1.2 Gy, applied in 1-58 fractions, were collected (see Tables 1 and 2, and Additional file 1 for details of the publication search).
Of all reported tumors, 63.6 % were confirmed to be stage T1, 36.4 % T2 (  ] years in the CF and HF groups respectively. Patients, who received conventionally fractionated RT were treated in the time period from 1976 to 2010, whereas patients treated with hypofractionated regimes were irradiated in the time period from 1996 to 2012. In the CF cohort only in one study PET-CT was performed for staging in 6 out of 31 patients (Bogart et al. [36]), whereas for many of the HF cohorts PET was a routine procedure; for many of the most recent studies PET-staging was even an inclusion criterion in the retrospective series.
In the 8 series of the CF group, generally a margin of 1-1.5 cm was added around the gross tumor volume (GTV), which was in some cases estimated from port films if no planning computer tomography (CT) scan was available. In the HF series, most frequently no GTV-to-CTV (clinical target volume) margins were added, except in 5 out of 23 series. Internal target volume (ITV) concepts were applied in 13/23 studies, based  BED 10 biologically effective dose with α/β = 10 Gy, PTV planning target volume, D total dose, d dose per fraction, T total treatment time, LC local control, ns not specified, dens inhom corr density inhomogeneity correction, PB convol pencil beam convolution, CS convolution superposition either on addition of GTV from 3D-CT scans in expiration/inspiration and free breathing in 4 cases, on slow CT scans in 5 cases, and on 4D-CT scans in 4 cases. Nine out of 23 studies did not apply any ITV concept. The most frequently used CTV-to-PTV (planning target volume) margins were 0.5 cm in axial and 1 cm in cranio-caudal directions. A minimal margin of 0.2 cm was added in one patient cohort treated with the Cyber-Knife where tumor tracking was used to correct for intrafractional target motion. In 2 out of 23 series the PTV margin definition was patient-specific. Different dose reporting concepts were found throughout the selected references. Only 5 CF reports out of 8 explicitly mentioned that the dose was prescribed to the isocenter. When not specified, prescription to the isocenter was assumed. In the case of the hypofractionated SBRT data, 9 references reported the prescribed dose to isocenter and 11 to the isodose line encompassing the PTV, which ranged from the 50 to the 100 % isodose, most frequently to the 80 % isodose line. Only one of the SBRT cohorts was treated with IMRT, and in this case the dose was prescribed to the 95 % isodose line enclosing the PTV. More information can be found in the Additional file 1.
Median [min -max] applied BED at the isocenter, calculated with an α/β ratio of 10 Gy (BED 10  Median value of the documented median follow-up times were 28.5 and 27 months for the CF and HF groups, respectively. Very few of these publications state explicitly the number of patients at risk at each follow-up time point (four cohorts). In only 5 of the selected studies an estimation of the 95 % CI of the calculated actuarial local control rates was reported and one publication presents the standard error. Therefore, no information could be collected about the precision of the estimated actuarial local control rates, apart from the cohort size and the median follow-up time.

3-and 5-year clinical outcomes for stage I NSCLC
The median value [range] for the 3-year LC was 86 [range: 55-100] % for the HF, and 63 [range: 50-91] %, for the CF series, respectively. Twenty-two out of 26 HF 3-year LC data points lie above the 80 % level, while all except two of the CF datasets lie below. All except one HF treatments delivered a BED 10 of 100 Gy or higher at the isocenter but none of the CF treatments reached a BED 10 of 100 Gy.
Spearman's correlation coefficients between dosimetric parameters (total absorbed dose and BED 10 ) and the different clinical outcomes were calculated. For the complete dataset and the total BED 10 at the isocenter, a significant Spearman's correlation of 0.716 for the local control with BED 10 was found. For the BED 10 evaluated at the PTV edge, a significant correlation of 0.638 was found between LC and BED 10. The total absorbed dose at either dose point however did not correlate with LC.
Modeling local control versus BED Linear-quadratic (LQ) model with α/β ratios fixed to 10 Gy (two-parameter fit) All 3-year LC data points for both the HF and CF treatments were fitted to a logistic model coupled with BED 10 at both isocenter and PTV edge. For the isocenter, the TCD 50 Fig. 1a, and all fit parameter values are summarized in Table 5. The p value was found to be < 0.05 for both TCD 50 and k simultaneously only for the model based on the full dataset.
The model for the complete dataset based on BED 10 at PTV edge yielded a TCD 50   These models are represented in Fig. 2a.
Linear-quadratic-linear (LQ-L) model with α/β fixed to 10 Gy (three-parameter fit) The alternative BED definition including a linear portion in the dose-effect beyond a certain fraction dose was tested, according to Eq. 3, with an α/β ratio fixed to 10 Gy. Fits were performed either based on the full dataset or on the HF subset alone. For the models based on BED 10 at the isocenter, the fit based on the full dataset generated a D t value of 11.0 Only two data points were found to be below this D t , which was too few to obtain a reliable estimate. Thus, the fit parameters did not produce significant p values. This fit is shown in Fig. 2b, together with the LQ fit for comparison.
When BED 10 doses at PTV edge were used, similar D t values were found, namely, 12.4 Gy for the complete dataset, and 9.9 Gy for the HF dataset (this fit is shown in Fig. 2b). Likelihood ratio tests showed no difference between LQ-L and LQ fits, independently of which BED 10 doses were used, at isocenter or PTV edge (see Table 5).

Correlation of local control with BED
To complete the modeling study, the correlation of the 3-year LC with BED under different assumptions for α/β equal to 5, 8.6, 10, 15 and 20 Gy was investigated (see Table 4). For the complete dataset and BED 10 values at isocenter, a Spearman's correlation of 0.716 (p < 0.0001) with BED 10 was found. For all other α/β ratios, correlation values increased marginally from BED 5 (r = 0.706) to a maximum at BED 10 and decreased again for BED 20 (r = 0.706), in all cases being significant. The Spearman's correlations based on the BED values at the PTV edge decreased with growing α/β values for the complete dataset, from 0.680 to 0.510, all of them being significant and consistently lower than the respective values for the BED 10 at the isocenter (see Table 4).
In contrast, for the HF subset, the correlation of the LC with the same series of BED α/β at isocenter increased minimally from 0.575 to 0.618, and also the values for BED α/β at PTV edge, from 0.531 to 0.601, all of them being statistically significant.   (Table 5).

Differences in the prediction of isoeffects
In order to translate the impact of the different BED model assumptions on clinical treatment schedules we calculated as an example the doses per fraction which would be necessary to reach selected BED 10 levels at the isocenter (Fig. 4). Under the assumptions of the LQ model with α/β values of 8, 10, or 15 Gy and the LQ-L model at the isocenter with α/β equal to 10 Gy and a D t of 11.0 Gy, in order to reach 100 Gy (BED), the maximum differences among models for one and two fractions are 10.5 and 4.5 Gy, respectively. Maximal differences remain below 3.3 Gy for treatments delivered in 3 or more fractions. For 200 Gy (BED), discrepancies between models in fraction size increase to 30.1 and 10.5 Gy for one and two fractions, respectively and remain below 5.6 Gy for 3 fractions and more.

Discussion
Review of clinical outcome data after radiotherapy treatment represents the only possibility to gather long-term information from large numbers of patients, which could serve as basis for statistical analysis for radiobiological modeling. However, this task presents a number of challenges since these datasets are intrinsically heterogeneous. Variability among radiotherapy centers applies to aspects such as target volume definition, dose prescription, planning concepts and delivery techniques with different precision levels. Additional effects hindering    precise dose-effect modeling are: the use of multiple and sometimes less appropriate dose calculation algorithms especially in the case of lung tumors (e.g. with limited heterogeneity correction), which may lead to misestimation of the absorbed dose. Additionally, outcome of different tumor stages, dose levels and fractionation schemes are frequently reported together. Furthermore, including historical cohorts implies dealing with changes over time in the standard diagnostic and therapeutic procedures, e.g. for staging, recurrence assessment and radiotherapy image guidance.
To counteract this variability, we applied strict inclusion criteria to the selected publications. The analysis is based on 3-year local control in order to depict the dose-response in mature outcome data, while maintaining a sufficient number of data points. To our knowledge this is one of the largest, most homogeneous patient collectives among similar studies. Through the combination of conventionally fractionated and hypofractionated data, a broad range of doses and fractionation schemes is covered, which is a further requirement in order to achieve conclusive modeling results.
Dosimetric heterogeneity in the PTV can be very pronounced in dose distributions for stereotactic treatments, reaching dose differences between isocenter and PTV edge of up to 50 %. It is not possible to know a priori, which reported dosimetric parameter will describe best the dose-effect relationship: the dose at the PTV edge or the dose at isocenter. Therefore, we calculated models based on both, isocenter and minimum target doses. The TCD 50 doses estimated for the isocenter doses were in general higher than the ones from the models at PTV edge.
Variability in the estimated doses at PTV edge could arise for instance, from variations in the CTV and PTV margin definitions among institutions, or uncertainties in the dose calculation methods, which in the case of outdated, less accurate dose calculation algorithms for the lung, would produce large dose mis-estimation and underdosage at the PTV edge.
We calculated the Spearman's correlations between outcome parameters and the BED 10 evaluated at both dose specification points and observed that the correlations with BED 10 at the PTV edge were without exception lower than the corresponding correlations based on BED 10 at the isocenter. This could be interpreted as an indication of the isocenter doses being more robust than the doses at the PTV edge for retrospective modeling studies, in agreement with previous studies [61].
The BED 10 fits based on the CF data and the complete dataset differ in both the TCD 50 and k values. This can be explained by the fact that the information required defining the slope of the dose-response curve and the TCD 50 , which together shape the sigmoid region of the logistic function, is provided by the CF data, where generally a lower BED was applied. This explains why TCD 50 and k are not consistently determined across models, which use different input data: CF + HF, or HF/ CF alone. We represented the range of the model input data and the extrapolation regions in all figures, to stress that special care must be taken predicting doses in the extrapolation region.
For comparison, Fig. 1 includes the previously published dose-response models of Martel et al. [28] and Guckenberger et al. [61]. Martel et al. found a TCD 50 , which is higher as compared to our findings as they  Logistic fits for LC versus BED α/β at isocenter, with constraint to approach the sigmoidal curves to the coordinate origin, a with α/β fixed to 10 Gy for the CF and HF datasets separately, and b with three free parameters: α/β, TCD 50 and k, for the combined dataset with conventionally fractionated and hypofractionated treatments. Logistic fits without constraint are provided for comparison compiled stage I/II data whereas our dataset included stage I tumors only. The values for the slope γ 50 that we obtained for the complete dataset and the CF subset are lower than the value in Martel et al. [28] and further values in the literature. For instance, Stuschke et al. [11] reported a value of 1.5, and Okunieff et al. [26] found a value of 1.6 on average.
Our fit based on the HF data alone shows a shallower curve than the fits including CF data, since the outcomes after hypofractionation are well above the inflection point of the sigmoidal curve. Our HF fit results for both BED 10 (isocenter) and BED 10 (PTV edge) are graphically similar to the findings by Guckenberger et al. [61], which were based on a multicentric compilation of individual patient SBRT data.
The LQ-L model was proposed to account for the experimental observation that curves of the (log cell survival) versus radiation dose often show a more straightened portion at doses beyond 10 Gy than predicted by the LQ model. This effect is in the majority of cases due to heterogeneity of the radiosensitivity distribution of the cells. This is partly due to differences in the stage of cells in the proliferation cycle, but can also be due to partially hypoxic conditions. Resistant cells tend to survive even at larger doses, which causes the survival curve to become less steep than predicted by the LQ model. This explanation implies that the straightening in the curve is not caused by a fundamental mechanism but by a simple to explain heterogeneity in the distribution of sensitivity. Both models converge at dose of less than about 6 Gy.
Applying state-of-the-art fitting methods to compare the performance of the LQ versus the LQ-L models for different fractionation schemes, we did not find significant differences. Therefore, it was not possible to decide, which model better predicts clinical NSCLC outcome data. For BED corresponding to doses per fraction below 11.0 Gy, the BED points for both concepts overlap, whereas above D t there is a contraction in the BED LQ-L values, which has no consequence in the fit itself, since it takes place in the region where the TCP approaches 100 %. Clinical consequences of using one model or another are only relevant for highly hypofractionated schedules aiming at delivering BED values well above 100 Gy.
The LQ-L model was also fitted to the HF datasets alone, although we could not obtain a 68 % CI for the TCD 50 , and the p values of all three parameters were > 0.05. Fig. 2b clearly demonstrates the similarity of our LQ-L fit for the HF data subset to the LQ fit previously presented, therefore TCP predictions will be similar with either model. Although our LQ-L fit based on HF data did not yield significant estimates, the results suggest a D t estimate in the same range of magnitude of 10 Gy. This fit was tested previously also on hypofractionated data alone by Guckenberger et al. [61]. Their dataset had a median dose per fraction at isocenter of 20.8 Gy with range  Gy. This group found a D t value of 22 Gy with a broad 68 % confidence interval,  Gy, whose addition to the model did not improve the prediction power.
Estimation of the α/β ratio by fitting three parameters simultaneously on a clinical dataset presenting high dispersion -as is the case of the current work -is challenging. Our LQ model with free α/β did not yield significant values for α/β, nor for TCD 50 and the slope k, which appears to be too shallow after visual inspection. We concluded that there is no indication for larger α/β values than 10 Gy if the complete range of fractionations is considered. The opposite trend (α/β > 10 Gy) was found for the HF dataset as was also the case in Chi et al. [15], although no significant p values could be obtained in this case, neither.
There are few works aiming to the estimation of a clinical α/β for radiotherapy of NSCLC [10,11,15]. Thames et al. [10] published an extremely high α/β value for lung tumors but these authors did not obtain a reliable confidence interval and so, their calculations must be regarded with caution. A similar work was also carried out by Stuschke et al. [11], who found an α/β value of 8.2 Gy. They used a fit similar to ours, but set a constraint to force the model to approach the axes origin by adding a point with low BED and null TCP, with a high fit weight. We also tested this approach (full model information in the Additional file 1), adding a data point at 0 Gy (BED α/β ) and 0 % TCP. We observed that this constraint influenced the TCD 50 value to a small extent only, but could have a strong effect on the slope of the curve, and also on the α/β value, for instance, 3.9 [2.2-9.0] versus 12.6 [10.5-15.0] Gy for the complete dataset and BED 10 (isocenter). The LQ-based fits with α/β of 10 Gy for the CF data alone did not vary much with and without constraint. In contrast, if the fit was based on the HF dataset alone, adding this constraint on the original fit had a pronounced effect on the steepness and TCD 50 of the curve, which in the constrained fits approached the curves based on CF data (Fig. 4a). For this reason, we think it is preferable not to set a constraint to the model, which largely influences the estimates for k and α/β, and also their standard errors. It seems reasonable to accept that fits done on the HF dataset alone will not reproduce a realistic fall-off in LC at low doses nor a plausible, clinical TCD 50 , since the input dose-response data are well above that region in this specific case.
The dispersion in the collected data points was high. Due to this fact, in the hypofractionation range the doseeffect relationship appears to be weaker than in other reports [13,61]. Specifically, for the fraction of HF data above 100 Gy (BED 10 ), the Spearman's correlation between LC and BED 10 is low. It can be speculated that in the region of high tumor control probability and highly hypofractionated treatments the relative contribution of non-radiobiological factors to the treatment effect is larger. Such factors are, for instance, subjectivity in target delineation and geographical miss, among others.
The logistic fit was applied also in the study at high dose regions since it is widely used for dose-effect description. However, other functional dependencies of LC with BED might be also appropriate in the high BED range. Concerning the BED concept, it seems likely that more than one approach can fit equivalently well on an inherently noisy dataset like this, for instance the LQ applying a higher α/β ratio, or the LQ-L model with α/β of 10 Gy, even if these models have contrasting radiobiological implications.