Computational delineation and quantitative heterogeneity analysis of lung tumor on 18F-FDG PET for radiation dose-escalation

Quantitative measurement and analysis of tumor metabolic activities could provide a more optimal solution to personalized accurate dose painting. We collected PET images of 58 lung cancer patients, in which the tumor exhibits heterogeneous FDG uptake. We design an automated delineation and quantitative heterogeneity measurement of the lung tumor for dose-escalation. For tumor delineation, our algorithm firstly separates the tumor from its adjacent high-uptake tissues using 3D projection masks; then the tumor boundary is delineated with our stopping criterion of joint gradient and intensity affinities. For dose-escalation, tumor sub-volumes with low, moderate and high metabolic activities are extracted and measured. Based on our quantitative heterogeneity measurement, a sub-volume oriented dose-escalation plan is implemented in intensity modulated radiation therapy (IMRT) planning system. With respect to manual tumor delineations by two radiation oncologists, the paired t-test demonstrated our model outperformed the other computational methods in comparison (p < 0.05) and reduced the variability between inter-observers. Compared to standard uniform dose prescription, the dosimetry results demonstrated that the dose-escalation plan statistically boosted the dose delivered to high metabolic tumor sub-volumes (p < 0.05). Meanwhile, the doses received by organs-at-risk (OAR) including the heart, ipsilateral lung and contralateral lung were not statistically different (p > 0.05).

to deliver increased dose to high [ 18 F]FDG uptake regions while maintaining the dose to other parts of the heterogeneous tumor by dose painting. Intensive research has investigated the relations between tumor sub-volumes with different metabolic activities and the survival outcome 12,[14][15][16] . It has been proved that the localization and analysis of sub-volumes in metabolic tumor volume (MTV) are necessary for dose painting and precision personalized RT treatment planning 9 using [ 18 F]FDG-and [ 18 F]flortanidazole (HX4)-PET/Computed Tomography(CT) imaging [12][13][14][15][16] .
The usual approach to define the boundaries of MTV and sub-volumes is by manual delineation. Manual delineation, however, is time consuming, labour intensive, operator-dependent 17 and subject to inter-and intra-observer variations 18,19 . Most of the existing computational segmentation methods 17,20 including our previous work 21,22 focused on the overall tumor segmentation for accurate cancer diagnosis purposes. As the delineation of complete tumor volume together with the measurement and extraction of metabolic sub-volumes from [ 18 F]FDG-PET images are the basis of effective personalized treatment planning, and in order to better interpret the different levels of requests from diagnosis and treatment aspects, we propose an automated computational technique that delineates the target tumor and in the meanwhile, identifies and measures metabolic sub-volumes. A dose-escalation plan using the computational delineation and measurement results is performed by intensity modulated radiation therapy (IMRT) system. The technique is proposed with hope to provide second opinion for personalized precision RT treatment planning with radiation dose-escalation based on [ 18 F]FDG-PET.

Materials and Methods
Patient studies. 58 [ 18 F]FDG-PET/CT scans were retrospectively acquired in lung cancer patients. All patients were scanned on a Discovery LS PET-CT system (GE Healthcare, Milwaukee, WI, USA) after the injection of [ 18 F]FDG. The patient studies enrollment criteria include: patients who received 60Gy (standard dose) and concurrent chemotherapy but had local recurrence; tumor exhibited inhomogeneous intensity distributions or indistinct boundaries. Among the patient studies, there are 30 cases with tumor infiltrating chest wall, 21 cases with tumor adjacent to mediastinum and 7 isolated cases. The PET images were reconstructed using a matrix of 144 × 144 pixels with a size of 4 mm × 4 mm × 4 mm; the CT images were reconstructed using a matrix of 512 × 512 pixels with a size of 1.17 mm × 1.17 mm × 5 mm. The PET images were up-sampled to the size of CT to obtain a spatial correspondence between the volumes in a common space. Manual tumor delineations were drawn on the registered PET with CT images as reference by two radiation oncologists (referred to as Obs-1 and Obs-2) for automated segmentation accuracy evaluation.
All the data used in this study were approved by Medical Ethics Approval Committee at Shandong Cancer Hospital. All the patients consented to participate in this study and have signed a written informed consent document. The methods were performed in accordance with relevant guidelines and regulations.
Automated target tumor segmentation. The proposed tumor segmentation framework is given in Fig. 1. The first step is to separate the target tumor from surrounding structures by extracting a 3D masking surface (3D-MS). Secondly, a tumor customized boundary delineation method is proposed where the 3D-MS will shrink in a hill-climbing manner towards the tumor center until a joint stopping criterion is satisfied. Stage 1. The first stage is to initially extract a 3D-MS for tumor separation. As high FDG uptake tissues such as tumor appear as salient regions standing out from the background, given an input thorax volume and user defined seed for target tumor localization, max intensity projection (MIP) images are firstly generated along x, y and z coordinate directions. For the cases (as shown in Fig. 1(b)) where the tumor is attached to neighboring high uptake regions such as the heart, a saddle point where two salient regions attached is firstly located on each MIP image. Then the corresponding tangent line is generated to separate neighboring attached salient regions. The isocontour ( Fig. 1(d)) passing through the saddle point is used to obtain a region of interest (ROI) mask. Finally, the seeds and ROI masks obtained from each of the MIP images are projected back to 3D space to generate 3D-MS.

Stage 2.
Given the detected 3D-MS, the second stage is to shrink the surface towards the true tumor boundary. The aim of tumor segmentation can be considered as to assign each of the voxels inside 3D-MS a foreground/ tumor label l F or background label l B . To obtain the shrinking stopping criterion, we define the joint GM and intensity affinity of two adjacent voxels v i , v j as where the intensity-based affinity term 0 is a rotationally symmetric Gaussian low pass-filter where σ is the standard deviation, d i denotes the distance from the origin in axis. The filter is applied to GM image to reduce the effect of noise while preserving the salient boundaries. In our experiments, β is set as 60 (−60) when 0 β > (β < 0) and σ is set as 1.5. It is noted that the absolute value of β in Equation (1) has no effect on the segmentation results because what matters is the condition of β being positive or negative. With φ defined as Equation (1), the segmentation is implemented in an iterative searching process by Algorithm 1.

Segmentation accuracy evaluation.
To evaluate the segmentation accuracy of the proposed method, we calculate the spatial overlap and shape dissimilarity between the automated segmented tumor volume (TV auto ) and manually segmented tumor volume (TV manual ) by Dice similarity coefficient(DSC) and Hausdorff distance (HD) respectively. DSC and HD are defined as: where S denotes the surface of a volume, sup represents the least subset element and inf denotes the greatest subset element; d is the Euclidean distance between points i and j. The DSC value is 1 for a perfect segmentation and a low HD value indicates high segmentation accuracy. The performance of our method was compared with (1) a threshold of 40 % SUV max , referred to as RG40; (2) a threshold of 50 % SUV max , referred to as RG50; (3) 28 . All the methods are implemented with MATLAB R2017a on a PC with 3.50 GHZ Intel (R) Core (TM) i7-4770K CPU and 16.GB memory, running a 64-bit Windows operating system. RG40, RG50, FCM and TCD were implemented within a fixed bounding box of 60 × 60 × 20 voxels centred at the initial seed. The user seed was provided by a person with knowledge to identify a tumor.
Quantitative heterogeneity computing for dose-escalation. Heterogeneous sub-volume extraction and measurement. Metabolic sub-volumes are further segmented within our segmented tumor volume (TV PM ) for heterogeneity measurement and dose-escalation purposes. We partition tumor sub-volumes with low, moderate and high metabolic activities (referred to as TV low , TV mod and TV high ) using respective SUV ranges 29 of (0,25%] SUV max , (25%,50%]SUV max and (50%,100%]SUV max . According to the calculation grid voxel size of the mainly used Treatment Planning System (TPS) 30 , a sub-volume whose size is smaller than 5 × 5 × 1 voxels is merged to the range of its nearest sub-volume and not counted as one isolated sub-volume.
Dose-escalation based on FDG-PET driven dose painting. The segmented tumor and metabolic sub-volumes are tested on Eclipse (Varian Medical Systems, version13.5) treatment planning system for dose-escalation in FDG-PET driven dose painting. Two six-field inverse-planned IMRT treatment plans are designed. One plan (plan A) is with uniform dose prescription of the entire tumor volume, the other comparative plan (plan B) is with sub-volume oriented dose-escalation based on our extracted metabolic sub-volumes. These two plans are set as: • Plan A: a 10-mm margin is added to GTV manual to create the planning target volume (PTV planA   Treatment plan evaluation. To compare the two treatment plas A and B, we calculate the dose distribution and dose-volume histogram (DVH) of PTV and organ-at-risk (OAR). OAR is manually delineated by experts on the planning PET/CT which includes the heart, ipsilateral lung, contralateral lung and spine cord. Two heterogeneity indexes (HI), HI RTOG = I max /RI and HI D2,D98 31 = (D 2% − D 98% )/D p × 100, are calculated to evaluate the two plans, where I max = maximum isodose in the target, RI = reference isodose, D 2% and D 98% represent the near-maximal and near-minimal dose respectively. If HI RTOG ≤ 2, the treatment plan is considered to comply with the protocol. A lower HI D2,D98 value indicates a less heterogeneous target dose.

Results
Automated tumor segmentation performance. The inter-observer variation between obs-1 and obs-2 was 0.824 ± 0.079 with respect to DSC and 8.99 ± 6.67 mm with respect to HD. As shown by Table 1, the proposed method (PM) achieved the best tumor segmentation results with DSCs of 0.850 ± 0.044 and 0.847 ± 0.052, HDs of 7.47 ± 3.14(mm) and 7.38 ± 2.71(mm), followed by LSE and LARW. The improvement was statistically significant (p < 0.05). The computation time of the whole procedure by PM depends on the size of the lesion and was about 0.82 seconds for one slice. The time taken in initial tumor separation was about 0.68 second per case. The average computation time in tumor segmentation was 0.63 for one slice.
Two cases are shown in Fig. 2 where the tumor in Fig. 2(a) is with irregular margins and lies adjacent to the myocardium with similar FDG uptake. The local tumor extent is difficult to discern on PET. RG40, RG50, FCM and LARW failed in the tumor separation while LSE resulted in smaller tumor delineation and RW failed to delineate the whole tumor. PM was able to accurately separate the complete tumor from the surrounding structures with the DSC of 0.858 (Obs-1) and 0.831 (Obs-2), HD of 9.179 mm (Obs-1) and 9.356 mm (Obs-2) respectively. For the case in Fig. 2(b), the tumor has heterogeneous FDG uptake reflecting regions of necrosis. RG40, RG50, FCM excluded the tumor regions with low SUVs and failed to delineate the whole tumor. PM obtained  Fig. 2(b). PTV high , PTV mod and PTV low are shown by red, yellow and blue colors respectively. Heterogeneity measurement and dose-escalation analysis. The case in Fig. 2(b) with high heterogeneity is visually and quantitatively measured and analyzed as shown by Fig. 3. In this squamous cell carcinoma case, red, yellow and blue color maps are used to indicate the location and spatial relations of TV high , TV mod , TV low in the ranges of (50%, 100%], (25%, 50%] and %(0, 25%] SUV max . The TV low s located near the tumor boundary were small, however, the intensities were greater than 23.86 % SUV max which could be considered as an acceptable value for tumor definition in clinical practice 32 . In contrast, the SUV values of the other three larger TV low s located closer to the tumor centre were as low as 8.24 % SUV max . It can be learnt from the spatial SUV distributions that there were significant variations between the distributions of high, moderate and less active sub-volumes within this tumor. For the case analyzed in Fig. 3, its DVH and isodose distributions of treatment plans A and B are given in Fig. 4. The tumor was located in the left upper lung and had denser and higher dose distributions in plan B than comparative plan A. As pointed by the white arrows in Fig. 4(a), the isodose contour of 20 Gy in plan B is much closer to the target volume than plan A, which indicated a lower dose received by contralateral lung. A decrement of radiation field width can also be observed in axial view. As shown by DVH in Fig. 4(b), plan B decreased the volume fractions of OAR when the isodose value was between 10 Gy and 60 Gy.
The dosimetry results of seven patients (characteristics listed in Table 2) with respect to PTV and OAR are given in Tables 3 and 4. All these cases complied with the protocol (HI RTOG ≤ 2). The PTV homogeneity was higher in plan A. As shown by Table 3, plan B resulted in a significant boost to the maximal (D 2 ) and average (D mean ) doses received by PTV (p-values < 0.05). The doses received by OAR in plan B (shown by Table 4) were not statistically significantly different from that in plan A (p-values > 0.05).

Discussion
It has been proved that distinguishing tumor tissues with metabolic heterogeneity would be helpful for tumor customized radiation dose painting and treatment 13,33 . With the development of IMRT, the dose escalation approach called simultaneous integrated boost (SIB) is able to deliver tailored doses to specific tumor sub-volumes with different metabolic levels 29 . Thus the accurate target tumor delineation and metabolic sub-volume extraction are important and essential in dose painting by contours or numbers in order to achieve high-quality RT results 34 .
We proposed an efficient and objective computational technique to outline the complete target tumor volume with heterogeneity, localize and measure tumor sub-volumes of different uptake ranges. The proposed method reduced the variability between inter-observers. Our major findings were: (1) The automatically extracted masking surface from MIP images was able to separate the tumor from its surrounding high FDG uptake regions and also provided an improved initialization for tumor delineation. Using the MIP images, the shape and size  Table 3. Comparison between various metrics in plans A and B with respect to PTV. D mean is mean dose; V 60 (%) = PTV volume receiving at least 60Gy; the p-value is calculated using two-tailed paired t-test. of the masking surface and the tumor were coherent. The first stage ensured that a contour is returned which topologically contains the ROI so as to avoid the exclusion of inhomogeneous areas within the tumor. The initial surface incorporated the whole tumor volume and avoided the variability in initialization. For those even more challenging cases when the tumor region has overlapping with surrounding high update regions on all the three projection views of MIP images, the first step with three directions might have difficulty in separation. In this case, we would suggest to increase the number of projection directions in our separation mechanism. (2) The proposed hill-climbing stopping criterion achieved tumor-orientated boundary delineations which complied with manual definitions. This was because the shrinking of the initial masking surface towards the tumor avoided the over-segmentation caused by tumor heterogeneity. In addition, the joint affinity model of GM and intensity was able to obtain a balancing tumor-specified boundary identification. In our experiments over 58 cases, we found that the SUV values of our defined tumor boundaries were between 23.65 % and 37.76 % SUV max , which also complied with the findings that a value ranging from 15 % to 50 % SUV max could be used for tumor delineation 32 . The quantitative extraction and measurement of tumor heterogeneity provided insightful information of tumor metabolism which assisted inhomogeneous dose calculation and treatment response monitoring and assessment. Our heterogeneity quantitative measurement results were tested on IMRT planning system to perform a dosimetry comparison between a dose-escalation plan and a standard uniform dose prescription plan. The evaluation results of seven patient studies demonstrated that the dose-escalation plan statistically significantly boosted the dose delivered to the tumor sub-volumes with high metabolic activities. In the meanwhile, the doses received by OAR including the heart, ipsilateral lung and contralateral lung were not statistically and significantly increased. The results revealed that by using the sub-volume oriented dose-escalation plan, the dose delivered to the target tumor can be boosted while preserving the dose delivered to OAR.
With the proposed method, the accurate segmentation of target tumor volume and measurement of heterogeneous sub-volumes would be beneficial to tumor metabolism analysis, heterogeneous dose calculation and escalation in FDG-PET driven dose painting by contours or numbers so as to improve the local disease control rate (DCR). The proposed model may also contribute to image-centric prediction models such as survival prediction and treatment response prediction following treatment by chemotherapy, radiotherapy, and hypnotherapy. Given the PET images obtained at different scanning time intervals or treatment stages, the proposed technique would be able to quantitatively track the re-distribution of metabolic sub-volumes. For patients with NSCLC, molecularly targeted therapy using epidermal growth factor receptor (EGFR) inhibitors is well developed. Thus, future studies may include FDG PET/CT data before and after targeted therapy and then analyze the dynamic changes of MTV in segmented areas as well as texture parameters. By such, it could continuously monitor and evaluate the therapy responses and thereby could make contributions to personalized treatment plans.  Limitations. The dose prescriptions in treatment plan B were designed to evaluate the proposed computational technique in dose-escalation to tumor sub-volumes with varying metabolic activities. Although the plan complies with protocol, it is not guaranteed that this is a most optimized plan. As this is a retrospective study, the treatment response using this plan in clinical practices has not been evaluated in this work.

Conclusion
We proposed an automated target tumor delineation and quantitative tumor heterogeneity analysis approach.
The experimental results proved that our method achieved improved segmentation accuracy. The dose-escalation plan using our heterogeneity measurement boosted the tumor sub-volumes with varying metabolic activities while controlling the dose livered to the organs at risk. This technique would offer an efficient way to assist the precise radiation therapy procedures, particularly for the lung cancer patients who need personalized dose painting to achieve high-quality RT results and deal with the local recurrence. The treatment response using the proposed method in clinical practices should be evaluated before its wide application.