Isobio software: biological dose distribution and biological dose volume histogram from physical dose conversion using linear-quadratic-linear model.

PURPOSE: To develop an in-house software program that is able to calculate and generate the biological dose distribution and biological dose volume histogram by physical dose conversion using the linear-quadratic-linear (LQL) model.
MATERIAL AND METHODS: The Isobio software was developed using MATLAB version 2014b to calculate and generate the biological dose distribution and biological dose volume histograms. The physical dose from each voxel in treatment planning was extracted through Computational Environment for Radiotherapy Research (CERR), and the accuracy was verified by the differentiation between the dose volume histogram from CERR and the treatment planning system. An equivalent dose in 2 Gy fraction (EQD2) was calculated using biological effective dose (BED) based on the LQL model. The software calculation and the manual calculation were compared for EQD2 verification with pair t-test statistical analysis using IBM SPSS Statistics version 22 (64-bit).
RESULTS: Two and three-dimensional biological dose distribution and biological dose volume histogram were displayed correctly by the Isobio software. Different physical doses were found between CERR and treatment planning system (TPS) in Oncentra, with 3.33% in high-risk clinical target volume (HR-CTV) determined by D90%, 0.56% in the bladder, 1.74% in the rectum when determined by D2cc, and less than 1% in Pinnacle. The difference in the EQD2 between the software calculation and the manual calculation was not significantly different with 0.00% at p-values 0.820, 0.095, and 0.593 for external beam radiation therapy (EBRT) and 0.240, 0.320, and 0.849 for brachytherapy (BT) in HR-CTV, bladder, and rectum, respectively.
CONCLUSIONS: The Isobio software is a feasible tool to generate the biological dose distribution and biological dose volume histogram for treatment plan evaluation in both EBRT and BT.

Purpose

Radiation therapy plays an important role in cancer treatment, either by external beam radiation therapy (EBRT) or brachytherapy (BT), or by a combination of the two [1]. The variation in the biological effects of radiation on cells or tissues depends on radiobiological factors such as cellular sensitivity and tissue organization. It is necessary to consider the biological parameters for interpreting and weighting in radiobiological models in 3D volume [2]. Linear quadratic (LQ) is the basic radiobiological model widely used for biological dose determination, and it consists of two components: the linear part of irreparable cell death and the quadratic part of cell death from no repair or incorrect repair of repairable components. The relation of the LQ model is S = e –( + , where D is the total dose, α is the slope of the survival curve, which is the effect of irreparable cell death, and β is the slope of the survival curve, which is the effect of cell death in repairable components that are not repaired [3]. However, the LQ model still has some limitations as it does not include the overall treatment time factor and repopulation of tumor during treatment, which can be resolved by the linear-quadratic-linear (LQL) model [4].

The true biological dose can be determined by the biological effective dose (BED), which is derived from the radiobiological model by including the biological parameters in the calculation. Biological effective dose is a numerical measures of treatment intensity, which is not equal to any prescribed dose of fractionation, and it is difficult to relate it to radiation tolerance dose in clinical practice [5]. Therefore, the biological effect should be normalized to the conventional dose, 2 Gy per fraction, or the equivalent dose in 2 Gy fraction (EQD2) form. Moreover, EQD2 is useful for the determination of the new appropriate dose per fraction in unplanned gap situations where cell repair during days off is considered, in order to maintain effectiveness and correct incomplete repair between fractions.

Recently, biological models for plan optimization and/or evaluation have been introduced to predict treatment outcome. A commercial biological treatment planning system is now available for research and clinic use, but it is optional. Thus, in this study, we aimed to develop an in-house software to generate biological dose distribution in terms of EQD2 by physical dose conversion using the LQL model. This software will be useful for evaluating treatment planning by considering both target coverage for tumor control and normal tissue’s dose for decreasing complication probability, thereby, increasing the efficiency of treatment and improving the quality of life for patients.

Material and methods

Linear-quadratic-linear model

The LQL model resolved the LQ model limitation by including the overall time factor and the incomplete repair (or recovery of the normal tissue) of sublethal damage in multi-fraction per day in the calculation. An increase in the overall treatment time affects the radiation tolerance of the early reacting tissues by repopulation. The survival fraction in the LQL model is described in Equation 1 where G(δD)is the reduction in the survival due to interaction between lesions [6]. The LQL survival curve in the large dose region undergoes a change from the continuous bending of the LQ survival curve to the linear curve described in Equation 2 [6]. The two different survival curves are shown in Figure 1. The solid line is the LQ survival curve and the dashed line is the LQL survival curve.

Fig. 1

The linear quadratic (LQ) and the linear-quadraticlinear (LQL) survival curves

Biologically effective dose and equivalent dose in 2 Gy fraction

Biological effective dose calculation in the LQL model is separated for two specific organs by considering the event for tumor and the organs at risk. The dose per fraction (d) is considered in relation to the LQL threshold dose (d), by following Voyant et al. [4]. In the target volume, when d is greater than d in the n fraction, the BED calculation using the LQL model is shown in Equation 3, where T is the potential doubling time in day, is the LQL model parameter, and θ(T – T) is the Heaviside function. This equation is useful for tumor proliferation correction when the overall treatment time is longer than the kick-off time of the tumor cell (T). At low doses per fraction (d less than d), the BED calculation is applied from the standard BED equation, and the incomplete repair of damage (φ) for multi-fraction correction is included as shown in Equation 4. H is the LQ correct for multi-fraction in the m fraction per day to correct for incomplete repair of damage.

For organs at risk, only the term of lag dose by proliferation is modified by recovered dose (D) for normal tissue, where D is and T is not included in calculation. When dose per fraction is greater than d, the BED calculation is described by Equation 5, and Equation 6 for low doses per fraction.

EQD2 is defined as “the dose in 2 Gy fraction that is biologically equivalent to the total dose D given with a fraction size of d Gy” [5]. The biological effect of any dose per fraction will be normalized to be equivalent with the dose in 2 Gy. Therefore, the advantages of EQD2 are that it is more related to every day clinical practice, it is useful for comparing the efficiencies of different treatment schedules in clinical use, and determining new appropriate doses per fraction in an unplanned gap situation. EQD2 based on BED is shown in Equation 7.

Software development

The in-house Isobio software was developed by using MATLAB version 2014b for calculating and displaying biological dose distributions and biological dose volume histograms. The treatment data exported in the Digital Imaging and Communications in Medicine (DICOM) file from treatment planning system (TPS) were required to be used with the in-house software. The image, dose, and structure data were extracted through Computational Environment for Radiotherapy Research (CERR) to be inputted and were located in the MATLAB workspace for use in biological dose calculation in the EQD2 form. The EQD2 values in each voxel were directly converted from physical dose that was extracted from the MATLAB workspace using the LQL model with the BED base. The EQD2 distribution as well as the EQD2 volume histogram was displayed. This process is shown in Figure 2.

Fig. 2

The equivalent dose in 2 Gy fraction (EQD2) function in the Isobio software

Patient cases

Treatment planning data of cervical cancer patients treated with EBRT and BT at Chiang Mai University Hospital were used in this study. In the EBRT plan, the four-field box technique was used for 3D conformal radiation therapy (3DCRT) planning from the Pinnacle version 9.8 treatment planning system (Nucletron, an Elekta company, Elekta AB, Stockholm, Sweden) with a prescription dose of 50.0 Gy (2 Gy in 25 fractions). Oncentra TPS version 4.3 (Nucletron, an Elekta company, Elekta AB, Stockholm, Sweden) was used for BT planning with a prescription dose of 7 Gy at D90% for high-risk clinical target volume (HR-CTV).

Statistical analysis

Most of the data that were randomly selected were normal distribution, while some data were found to have slightly deviated when the normal distribution of the data was determined using the Shapiro-Wilk test. Therefore, the pair t-test statistical analysis was used for software verification to ensure the accuracy of the physical dose that was extracted through CERR, and the biological dose conversion by using IBM SPSS Statistics version 22 (64-bit).

Results

Software features

The features of the in-house Isobio software for the EQD2 function are shown in Figure 3. The EQD2 function page was divided into two components. First, the input of the treatment planning information, dose per fraction, number of fraction, fraction per day, and treatment time are the important parameters required to input for BED calculation. The number of structures aligned following CERR is also required because the operation is a link between the in-house software and CERR. Another component is the display function that could display 2D in three planes: transverse, coronal, and sagittal. This software can also display the iso-biological dose as well as the organ’s surface in 3D view. Moreover, the biological parameters used for calculation can be adjusted by the user. The iso-biological dose lines are displayed using different colors by interpolation dose in each voxel of each slice.

Fig. 3

The equivalent dose in 2 Gy fraction (EQD2) calculation and display function page in the Isobio software

Software verification

Physical dose verification

The accuracy of the physical dose that was extracted through CERR of the in-house software was verified using the difference in dose-volume histogram (DVH) between CERR and TPS. The comparison of DVH in the various organs of interest in Pinnacle is shown in Figure 4A, and that of Oncentra is presented in Figure 4B. Table 1 shows the difference between the DVH found in Oncentra and Pinnacle. There was a larger difference in Oncentra (3.33% at D90% of HR-CTV), and negligible difference (less than 1%) in Pinnacle.

Fig. 4

The differentiation of the dose volume histogram in each organ between treatment planning and Computational Environment for Radiotherapy Research (CERR): (A) Pinnacle treatment planning system (TPS) version 9.8 vs. CERR, and (B) Oncentra TPS version 4.3 vs. CERR

Table 1

Percentage differentiation of dose-volume histogram (DVH) between Computational Environment for Radiotherapy Research (CERR) and Pinnacle version 9.8 for external beam radiation therapy in D95% and D2%, and Oncentra version 4.3 for brachytherapy in D90% and D2cc from cervical cancer treatment plan

	Region of interest	Percentage dose difference
EBRT
D95%	CTVPTV	0.190.18
D2%	BladderRectumBowel	0.220.240.18
BT
D90%	HR-CTV	3.33
D2cc	BladderRectum	0.561.74

EBRT – external beam radiation therapy, BT – brachytherapy, CTV – clinical target volume, PTV – planning target volume

Biological dose verification

The EQD2 calculation of the in-house Isobio software was verified using the pair t-test statistics with confidence interval (CI) 99% between the output data from the software and the output data from the manual calculation, as shown in Table 2. The percentage differences in EQD2 of EBRT and EQD2 of BT in each organ of interest in cervical cancer treatment were separately analyzed. There was a 0.00% difference found in EBRT and BT, with no significant difference in the p-value> 0.01, at 99% CI. In addition, the difference in BED in EBRT and BT was considered as BED was an important parameter in EQD2 calculation. Similar to EQD2, the difference between the software calculation and the manual calculation for BED was also insignificant.

Table 2

Percentage differentiation of biological effective dose (BED) and equivalent dose in 2 Gy fraction (EQD2) verification from external beam radiation therapy (EBRT) and brachytherapy (BT) between software calculation and manual calculation in cervix, bladder, and rectum for cervical cancer treatment plan

Organ of interest			n	Cervical (HR-CTV and IR-CTV)	Bladder	Rectum
EBRT	BED	ManualxMB ± SE	100	64.16 ± 0.08	77.21 ± 0.31	72.94 ± 1.32
		SoftwarexSB ± SE	100	64.16 ± 0.08	77.21 ± 0.31	72.94 ± 1.32
		% dif.		0.00	0.00	0.00
		p-value		0.717	0.158	0.171
	EQD2	ManualxMQ ± SE	100	53.47 ± 0.06	53.45 ± 0.21	48.21 ± 0.87
		SoftwarexSQ ± SE	100	53.47 ± 0.06	53.45 ± 0.21	48.21 ± 0.87
		% dif.		0.00	0.00	0.00
		p-value		0.820	0.095	0.593
BT	BED	ManualxMB ± SE	100	132.55 ± 11.68	13.64 ± 1.22	14.57 ± 1.25
		Software xSB ± SE	100	132.55 ± 11.68	13.64 ± 1.22	14.57 ± 1.25
		% dif.		0.00	0.00	0.00
		p-value		0.185	0.319	0.066
	EQD2	ManualxMQ ± SE	100	110.46 ± 9.73	9.44 ± 0.84	9.63 ± 0.83
		SoftwarexSQ ± SE	100	110.46 ± 9.73	9.44 ± 0.84	9.63 ± 0.83
		% dif.		0.00	0.00	0.00
		p-value		0.240	0.320	0.849

xMB – the average value of manual BED calculation, xSB – the average value of software BED calculation, xMQ – the average value of manual EQD2 calculation, xSQ – the average value of software EQD2 calculation, SE – the standard error

Clinical results

Biological dose distribution and biological dose volume histogram

The biological dose was converted voxel by voxel from the physical dose. The distribution of the biological dose could represent the position of the tissue receiving the dose, which is already corrected for biological parameters that are displayed in 2D view, as shown in Figure 5, as well as 3D view of both EBRT and BT are shown in Figure 6. In addition, the biological doses were plotted against the volume of organs of interest, and shown as EQD2 volume histogram (EQD2VH). EQD2VH is a histogram relating the biological dose to the volume of organ of interest, as shown in Figure 7A for EBRT and Figure 7B for BT. The D95% of PTV and the D2% of the bladder and the rectum for EBRT, as well as the D90% of HR-CTV and the D2cc of the bladder and the rectum can be demonstrated by EQD2VH, as shown in Table 3.

Fig. 5

The 2D physical dose distribution of (A) external beam radiation therapy (EBRT) (B). Brachytherapy (BT) and biological dose distribution from physical dose conversion using the linear-quadratic-linear (LQL) model of (C) EBRT and (D) BT

Fig. 6

The 3D biological dose distribution from cervical cancer treatment in each organ of (A) external beam radiation therapy and (B) brachytherapy

Fig. 7

Dose-volume histogram (DVH) and equivalent dose in 2 Gy fraction volume histogram (EQD2VH) in each organ of cervical cancer treatment planning from (A) external beam radiation therapy and (B) brachytherapy

Table 3

Equivalent dose in 2 Gy fraction (EQD2) of region of interest in external beam radiation therapy and brachytherapy

	Region of interest	EQD2 (Gy)
EBRT
D95%	PTV	51.76
D2%	BladderRectum	57.3855.32
BT
D90%	HR-CTV	9.72
D2cc	BladderRectum	7.537.78

EQD2 – equivalent dose at 2 Gy, EBRT – external beam radiation therapy, BT – brachytherapy, CTV – clinical target volume, PTV – planning target volume

Discussion

The Isobio software is suitable for calculating and generating biological dose distribution and biological dose volume histograms; thus, it is useful in treatment plan evaluation. The difference between the physical dose volume histograms was about the same as obtained by Davenport DA [7], which showed about 3-5% difference between CERR and TPS, and resulted from the dose grid resolution and the calculation algorithm to fit DVH [8]. The difference was in the acceptable error range about research and clinical use. Adjusting the dose grid resolution in CERR to the same value of TPS could decrease this difference. Also, the structure’s volume calculation of TPSs and CERR affected the DVH accuracy, especially in the small volume [9]. Therefore, a larger difference between Oncentra and CERR was observed when D2cc of OARs was considered. The larger difference found in HR-CTV might be a result of the high dose gradient of BT in the region near the radiation source.

The values obtained from the EQD2 calculation with BED base using the LQL model between the software calculation and the manual calculation were not significantly different in both EBRT and BT. This equation corrected only the repopulation and the repair factors, and did not account for the redistribution and the reoxygenation correction factors. The overall treatment time factor was included in the EQD2 calculation. This parameter should be included in prolonged treatment time.

The biological dose distribution in 2D and 3D is useful with regard to consideration of tumor coverage, and over and under dose point in both target and OARs, as well as for volume effect determination. Volume effect depends on the dose and the volume of the received dose, because increasing the dose may also increase the severity of effect or the frequency of incidence, or both, in normal tissue [4].

The EQD2VH was obtain from the DVH computation in CERR. The structure contour could be approximated as small volumes associated with image. The value of the dose in every voxel inside the contour polygon was linearly interpolated for the center of the voxel. The algorithm can fail to compute when the contour polygon is bend within part of the voxel, and when there is a change in the longitudinal space or slice thickness. Likewise, the DVH computation in commercial TPSs uses an interpolated dose in each voxel and a calculation volume for DVH generation, which depends on the CT characteristics (e.g., slice thickness and pixel width), dose grid resolution, and TPS manufacturer, which influence the DVH uncertainty. The percentages uncertainty of DVH computation from various TPSs were found to be 5% and 2% of the volume expectation and the dose calculation, respectively, by comparing to 3D matrix trilinear (X, Y, and Z directions) interpolation in SWAN system [10]. The SWAN system can decrease the uncertainty from a longitudinal direction, which includes the dose in half image-slice regions outside the superior and the inferior of each volume. However, the interpolation uncertainty from the contour boundary remained in CERR and TPS as it depends on the voxel resolution, DVH computation algorithm, and the partial voxel in the boundary of the contour polygon [9, 10]. EQD2VH is suitable for treatment planning evaluation, determination, and selection of the best treatment plan when different doses per fraction are given [11]. EQD2VH gives different dose response curves in tumor and OARs due to the different biological parameters weighting. It is easier to discuss the probability of normal tissue complication when tumor control remains equal to other plans. In contrast, the DVH gives the same dose response curve regardless of changes in dose per fraction, and this makes selecting the best plan a hard decision. D2cc in the organ of interest from biological volume histograms should be carefully interpreted as the anatomical position of dose might not be the same in the different fractions.

The biological parameters in this software were modified from LQL_equiv that were collected from seven radiotherapy treatment centers in France [6]. Before using this software to predict the treatment outcome, the user should understand and check the biological parameters carefully by considering the biological effect end point. Organ contours should also be carefully and accurately determined as they are the most important step used to identify the biological parameter of each voxel in biological dose calculation. This in-house software was generated from MATLAB code; therefore, it is more convenient for the user to apply their own biological parameters.

It might be concluded that the biological effect base is superior in treatment plan evaluation to physical dose base as, in the former, tissue specificity and other biological parameters related to the radiobiological model were included.

Conclusions

The in-house Isobio software is suitable for calculating the biological dose as EQD2 using the LQL model with BED base from physical dose conversion and for generating the biological dose distribution in 2D and 3D as well as biological dose volume histograms for treatment plan evaluation. It is suitable for use in both EBRT, including advanced techniques such as stereotactic body radiation therapy (SBRT) and stereotactic radiosurgery (SRS), and BT. In the future, we aim to improve this software in order to achieve better clinical treatment outcome prediction, for example, tumor control probability (TCP) and normal tissue complication probability (NTCP).