MATLAB®-based fitting method to evaluate survival fractions after multimodal treatment.

To easily analyse and visualize cell kill dynamics measured by survival fraction after single or combined treatments a MATLAB®-based application was developed. A statistical analysis with different options of visualisation of single and combined treatment effects can be performed in a few steps not requiring advanced knowledge of statistical programs.

Introduction

Clonogenic assays are used to quantify cellular response to anticancer agents including radiation [1]. The steepness of the initial slope of the survival curve correlates with clinical responsiveness [2], [3]. Therefore, clonogenic assays are frequently used tests in radiation research [4], [5]. However, the graphical and statistical evaluation of multimodal treatments is challenging because the mode of interaction of response modifiers and irradiation as well as the used radiation fraction size might limit the applicability of the LQ-model.

To analyse the effect of irradiation plus an additional modality also becomes tedious because a normalization of data points has to be performed for each treatment [6]. Employment of statistic software for evaluation is possible but programming of survival curve evaluation tools might become relatively complex and time consuming for a radiation oncologist or a radiobiologist as illustrated elsewhere [7].

Therefore, we developed a standalone GUI (graphical user interface)-application to analyse radiosensitivity measured by survival fraction after single and combined treatments in vitro. The LQ-Script was created in MATLAB®, which is a programming language available in many academic and research institutions providing reliable and proven algorithms and their implementation in GUIs [8]. The MATLAB®-based GUI application and their capabilities are demonstrated to give researchers of radiobiology an insight and easy access to relevant information on this script. Furthermore, this report can be utilized as an instruction, how to analyse combined treatments using survival fractions by this MATLAB®-based LQ-Script.

Methods

MATLAB® program

A MATLAB®-based program was developed to estimate α/β-parameters and to visualize approximated surviving fraction (SF) curves. System requirements were MATLAB® 2010 or newer. The linear quadratic (LQ)-Model [1], [9], [10], [11] was used to calculate cell kill kinetics after irradiation. The Eqs. (1), (2) describe the LQ regression with error term (±δ i.e. mean zero and homoscedastic variance σ2):

D is the radiation dose in Gy. On a log-linear plot of the survival curve, α is the cell kill per Gy of the linear component and β is the cell kill per Gy2 of the quadratic component. In assumption that in case of combined treatment the drug has it owns effect the following equation was applied:

Plating efficacy (PE) was calculated according to the following equation [7]:

Survival fraction (SF) was calculated according to the following equation

The script uses fit functions of MATLAB® including confidence intervals and prediction bounds. The algorithm is described elsewhere (http://uk.mathworks.com/help/curvefit/fit-postprocessing.html). The confidence interval and prediction bounds can be selected. Level of certainty is 0.95 (=95%). For users who have installed the “Statistics and Machine Learning Toolbox” additional options are available including standard errors for the estimated coefficients α, β, and SF0 along with t-statistics and p-values. The coefficient of determination (R2) is calculated as a goodness-of-fit measure and overall significance is measured by F-statistics. Both are given as output statistics. For advanced users, we integrated a feature allowing for plotting residuals of the fitted models. It is named “plot residuals”.

Data set

As example for combined treatment, we utilized a modified colorectal cancer cell line of HCT-116 which was cultivated as described elsewhere [12]. After graded radiation doses (0Gy, 2 Gy, 3.5 Gy or 5 Gy) with and without simultaneous in vitro cytotoxic drug (incubation with topoisomerase inhibitor type I) the colony forming assay was performed in delayed plating technique [13]. The plating efficiency (PE) and surviving fraction (SF) were calculated by the script as described above.

Script installation

The Supplement contains the Zip-File “Matlab_based_script.zip”. Please, unzip the file to a separate folder. The following five files are included:

LQFit.m – the MATLAB®-based LQ-Script

SampleData.mat – a MATLAB® sample data file

Main.Fig. and Main.m – graphical user interface of the script

InitLQFitStruct.m – a script facilitating data structure input

Results

A MATLAB® based script named “LQ-Script” was developed enabling visualization and statistical analysis of combined treatment effects. The instruction (development of data source structure, workflow and an example file of LQ-Script were attached (Supplement).

Algorithms

Two fitting methods are available. The first one is the classical linear fitting method, implemented as a build-in MATLAB® function named fit called with the parameter Method set to LinearLeastSquares. This method uses the linear least squares method optionally with a robust parameter (bisquare) to provide LQ-model approximation. The second one is the nonlinear least squares curve fitting method implemented as the same build-in MATLAB® function fit, called with parameter Method set to NonLinearLeastSquares.

Both methods use as default the unmodified MATLAB® implementation of the Levenberg–Marquardt algorithm [14], [15]. Optionally the Trust-Region algorithm could be used. The script sets the following options to run the nonlinear least squares method, Table 1.

Table 1

Initial conditions for the Levenberg–Marquardt algorithm. The initial guess of the parameter has to be provided to start running the nonlinear least squares method.

Parameters of the nonlinear least squares method	Item/value
Initial values for the coefficients	α:0
	β:0
	sf0:1
Algorithm	Levenberg-Marquardt
Maximum change in coefficients for finite difference gradients	0.1
Minimum change in coefficients for finite difference gradients	10e−8
Maximum number of evaluations of model allowed	600
Maximum number of iterations allowed for fit	400
Termination tolerance on model value	10e−6
Termination tolerance on coefficient values	10e−6

Curve visualization

A MATLAB®-based fit was approximated for an exemplary data set of combined irradiation and cytotoxic drug in vitro. The fitcurves were approximated und visualized on the basis of surviving fractions after irradiation. The number of data points which can be evaluated is not restricted.

In the “Main” window, the following evaluation tools can be selected (Fig. 1A): At first, the cell line(s) and drug concentrations are chosen for the evaluation by using “String” and a single mouse click (i.e. the sample “Some cell line” is chosen and the two drug concentrations 0 and 2.5 are selected). The next step is to define the regression type (linear/nonlinear, see above). Three different LQ-analyses are available:

Fig. 1

Curve visualization and statistics. A) “Main” window shows selectable evaluation tools. B) The calculated statistics is displayed in the “Command window”. C) In case of choosing “show additional statistics” the calculations are also displayed in the “Command window”, too.

The first option “LQ-Analysis normalized by the control” can be used to derive LQ-curves with normalization to control (i.e. 0 Gy) with and without a drug (Fig. 1A and A) to evaluate log-additive effects. The second (classical) option “LQ-Analysis normalized by SF0 [Gy]”. In this analysis the normalization to SF0 = 1 is performed as typically used for analysis of survival fractions and can be applied for evaluation of supra-log-additivity for instance if a response modifier is used (Fig. 1A and B). The third option “LQ-Analysis” equals LQ analysis without normalization to visualize the distribution of raw data (Fig. 1A and C).

The following display options are available and refer to colony counts. Options include “Plot mean value”, “Plot Confidence Interval” (CI), “Plot residuals”. The degree of certainty (usually 0.95) can be manually adapted for CI as well as for prediction bounds. “Plot measure points” and “Plot prediction bounds” (PB) are only available for “LQ-Analysis” (raw data). The robust regression can be selected either as “LAR” (least absolute residual method) or as Bisquare. The corresponding R2 is calculated and shown in the command window (Fig. 1B). The best fit (R near 1) should be selected. In addition, α, β, SF0 and goodness of fit are calculated and displayed. Using the box “show additional statistics” the F statistics, mean, standard error (SE), t-statistics, p-values are shown. Furthermore, confidence interval (CI), survival (%), Plating efficacy (PE in%) and CI of PE are calculated and displayed (Fig. 1C).

All permutations of fitcurves can be visualised (in our example irradiation alone, combined treatment at each radiation or chemotherapeutic dose). A 95%-confidence interval (CI) was chosen to graphically compare data, Fig. 2. Furthermore, the diagrams can be customized by a built-in Graph Editor with export options to different image file formats (e.g. .JPG or .PNG).

Fig. 2

Visualization of log-additive, supra-log-additive effects and raw data. A) Normalization to each control (i.e. with and without drug): The difference at 0 Gy between both curves demonstrates the log-additive effect of a drug and is indicated by an arrow. A significant log-additive effect is indicated by the dark blue area until the crossing (second arrow) of the confidence interval of both curves (red/green CIs). B) Normalization to 100 = 1 (classical analysis of radiosensitivity without a drug): A significant supra-log-additive effect is measured until crossing of both CIs of the radiation response curves with and without drug in a semi-logarithmic scale. In the example, the CI (red/green CIs) cross before the first applied radiation dose at ∼1.5 Gy (indicated by the blue area). Therefore, no supra-log-additive effect was measured since no radiation dose has been tested between 0 and 1.5 Gy. C) Raw data can also be visualized in a semi-logarithmic scale (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.).

Combined treatment i.e. combined irradiation and cytotoxic drug can be demonstrated in two different modes to measure additivity and supra-additivity. Additivity is defined as the sum of both single effects of two modalities i.e. addition of a drug to irradiation results in a parallel shift of the survival curve. Supra-additivity is defined as an effect being larger than the addition of single effects of two modalities [16]:

Normalization to each control (i.e. with and without drug): Survival fractions were normalized to each control (the “log”-additive effect of combined modalities was visualized in a semi-logarithmic scale). Heterogeneous response to multimodal treatment was visualized using the option to start the combined treatment curve at the level of cell kill related to the second modality (i.e. drug dose levels) at 0 Gy. A significant additive effect is indicated by the dark blue area until the crossing (arrow at ∼ 2.4 Gy) of the confidence intervals (CI) of both curves, Fig. 2A.

Normalization to 100=1 (classical analysis of radiosensitivity without a drug or to show supra-log-additive effects of a drug): The combined modality curve was normalized to 100=1 in order to visualize a potential radio-sensitizing (supra-log-additive) effect of the cytotoxic drug depending on concentration of chemotherapy. A potential supra-log-additive effect is measured until crossing of both CIs in a semi-logarithmic scale, Fig. 2B. In the example, the CIs (red/green) cross before the first applied radiation dose at ∼1.5 Gy (indicated by the blue area and the arrow). Therefore, no supra-log-additive effect was measured since no radiation dose has been tested between 0 and 1.5 Gy.

No normalization can be chosen to visualize raw data in a semi-logarithmic scale, Fig. 2C.

Taken together, log-additive and supra-log-additive effects (using different normalization methods, Fig. 2) can easily be graphically visualized and distinguished with this application.

Discussion

We developed a new software tool (LQ-Script) dedicated for visualizing combined treatment effects. All needed parameters are point-by-point requested from the user ensuring a comprehensive data visualization. Furthermore, a script for structure level definition is available to avoid challenging programming on the command line [7]. The user does not need to decide which parameters are needed, which coefficients should be used nor the analysis model. We demonstrated the functionality of the script. The resulting diagram enables a graphical comparison of curves in regard to statistical significance by evaluating confidence intervals and prediction bounds. Moreover, statistical features like F-statistics and goodness of fit were implemented. Supra-additivity as well as additivity of combined treatment can be evaluated by two modes of normalization. The α/β-value can be calculated allowing to annualize the effect of different fractionation schedules.

However, the possibility to present data on the basis of the LQ-algorithm should not lead to the assumption that combined treatment (i.e. radiation treatment and a response modifier) always follow the LQ-mode. The LQ-model can safely be applied to hyper-fractionated [17], [18] and conventional fractionated radiation schedules [19], [20]. There are some lines of evidence that the LQ-model can also cautiously be used for moderate and extreme hypofractionation including stereotactic body radiotherapy (SBRT) [21], [22]. However, criticism exists due to the not regarded impact of radioresistant subpopulations [23] and underestimation of isoeffective doses of normal tissue with large fraction sizes (>9 Gy) [24]. The applicability of the LQ-model can also be affected depending on mechanism of action of a response modifier and depending on culture conditions [25]. Depending on the choice of radiation dose range, another challenge could occur. In case of a strong radiosensitizing effect and radiation doses close to a potential threshold of the LQ-model, the threshold of the LQ-model could also be exceeded.

In summary, this MATLAB®-based script offers features enabling a graphical evaluation of surviving curves and permits to investigate the influence of radiation response modifiers. However, limitations of the LQ-model and modes of interaction have to be taken into account to carefully evaluate combined treatment effects. In addition, MATLAB@ is a software usually available at academic sites. However, biannual updates of MATLAB@ could lead optical changes of the LQ-Script.

Conclusion

An easy to use standalone MATLAB®-based GUI-application has been created visualizing single and combined treatment effects. The developed software resolves current limitations of other available templates of spreadsheet analysis programs like restrictions concerning total number of evaluated cells, combined treatment effects and comparison of different sample sizes. A complex data analysis can be conveniently performed by the MATLAB®-script. Graphical comparisons of different data sets become much easier. The wide range of options (presentation of all measure points and the mean, 95% CI, PB) facilitates data analysis. The clearly arranged, expandable input data structure can be used to perform some additional data processing. The software is available in the Supplement.

Conflict of interest

None.

Authors’ contributions

JM has made substantial contributions to conception, design, and discussion of the software script and writing of the manuscript. LM has developed the MATLAB®-based script and made substantial contributions to writing of the manuscript. PM has tested the MATLAB®-based script and made substantial contributions to evaluation of statistical tools and application of algorithms i.e. statistics. DZ was involved in revising the manuscript for important intellectual content and gave final approval of the version to be published. ACM has made substantial contributions to conception, design, and discussion of the software script and writing of the manuscript.

All authors read and approved the final manuscript.