Combined morphologic and metabolic pipeline for Positron emission tomography/computed tomography based radiotherapy response evaluation in locally advanced pancreatic adenocarcinoma.

BACKGROUND AND
PURPOSE: Adaptive radiation planning for pancreatic adenocarcinoma (PA) relies on accurate treatment response assessment, while traditional response evaluation criteria inefficiently characterize tumors with complex morphological features or intrinsically low metabolism. To better assess treatment response of PA, we quantify and compare regional morphological and metabolic features of the 3D pre- and post-radiation therapy (RT) tumor models.
MATERIALS AND METHODS: Thirty-one PA patients with pre and post-RT Positron emission tomography/computed tomography (PET/CT) scans were evaluated. 3D meshes of pre- and post-RT tumors were generated and registered to establish vertex-wise correspondence. To assess tumor response, Mahalanobis distances (M dist∣Fusion) between pre- and post-RT tumor surfaces with anatomic and metabolic fused vectors were calculated for each patient. M dist∣Fusion was evaluated by overall survival (OS) prediction and survival risk classification. As a comparison, the same analyses were conducted on traditional imaging/physiological predictors, and distances measurements based on metabolic and morphological features only.
RESULTS: Among all the imaging/physiological parameters, M dist∣Fusion was shown to be the best predictor of OS (HR = 0.52, p = 0.008), while other parameters failed to reach significance. Moreover, M dist∣Fusion outperformed traditional morphologic and metabolic measurements in patient risk stratification, either alone (HR = 11.51, p < 0.001) or combined with age (HR = 9.04, p < 0.001).
CONCLUSIONS: We introduced a PET/CT-based novel morphologic and metabolic pipeline for response evaluation in locally advanced PA. The fused M dist∣Fusion outperformed traditional morphologic, metabolic, and physiological measurements in OS prediction and risk stratification. The novel fusion model may serve as a new imaging-marker to more accurately characterize the heterogeneous tumor RT response.

Introduction

Pancreatic adenocarcinoma (PA) has one of the worst outcomes of all solid tumors with a 5-year survival rate of 7% [1], [2]. At the time of diagnosis, only 10% − 15% of patients are eligible for resection, which is the only therapy offering potentially prolonged survival [1]. A major contraindication for surgery is vascular involvement by primary tumor, which precludes a margin-negative resection. As a result, approximately 40% are considered borderline resectable or locally advanced (BRPA, LAPA). Following aggressive multiagent chemotherapy and radiotherapy, only <30% of patients are downstaged [3]. However, an accurate post-radiotherapy assessment with computed tomography (CT) remains challenging and may not accurately reflect tumor response [4].

To assess tumor treatment response, size-based criteria, such as Response evaluation criteria in solid tumors (RECIST) and World Health Organization (WHO), have been widely adopted in current clinical practice [5], [6]. While easy to implement and generalize, these uni-dimensional measurements inaccurately characterize tumors with complex geometry (Fig. 1 (A)) and are insensitive to tumor heterogeneity, which is particularly important for dense stromal tumors such as PA [7]. For example, in a study investigating the association between RECIST response and overall survival (OS), similar median OS durations were observed in patients who underwent resection and had no radiographic evidence of response to neoadjuvant therapy and patients who did respond [4]. In addition to size-based criteria, F Fluorodeoxyglucose (FDG) positron emission tomography (PET) has been commonly used to quantify tumor metabolic activities [8], [9], [10]. However, basic PET values, such as the maximum or mean FDG uptake, is inadequate to describe the heterogeneous tumor metabolism (Fig. 1 (B), P5).

Fig. 1

A: Representative tumor short and long axes on axial CT images and the corresponding 3D constructions for a 64-year-old female (P1), a 76-year-old male (P2) and a 77-year-old female pancreatic adenocarcinoma (PA) patient (P3). All 3 patients show similar values of tumor short and long axes [5], [6], while their 3D tumor constructions revealed varied shape features. B: Representative PET/CT fusion images of PA patients with varied SUV distributions within tumors. Boundaries of tumors are circled in red. P4, Axial (left column), coronal (middle column), and sagittal (right column) views of a 77-year-old female patient with low SUV PA tumor. P5, a 61 years old male patient with heterogeneously distributed SUV within PA tumor. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Moreover, not all tumors present with increased FDG uptake compared to surrounding normal parenchyma [8] and comorbidities such as hyperglycemia and diabetes, which are common in PA, can disrupt tissue FDG uptake leading to false-negative results [11]. Therefore, traditional assessments based on current standard of care functional (SUVmax, SUVmean, etc.) and structural imaging (RECIST, WHO criteria) are inadequate to evaluate treatment response in PA. Integrated anatomic and metabolic assessment may provide a more comprehensive and accurate indication of tumor response [10]. In this study, we integrate morphologic and metabolic metrics of LAPA and BRPA using CT/PET to assess treatment response and clinical outcomes using multivariate distances between the resulted 3D pre- and post- RT tumor surfaces. Our aims are: first, to develop a combined 3D anatomic and metabolic pipeline; second, to test the efficacy of this pipeline as a quantitative survival predictor for patients with LAPA or BDPA.

Materials and methods

Subjects and data

With Institutional Review Board approval, 31 patients with biopsy-confirmed PA were retrospectively selected from our institutional database between 2011 and 2017. All patients had LAPA/BRPA treated with radiotherapy and received both pre and post-RT PET with diagnostic CT scans. The post-RT PET/CT was taken 6–8 weeks after radiotherapy.

Patients were instructed to fast for a minimum of 6 h and then undergone skull to thigh imaging with a PET/CT scanner. Three scanners were included in this study: one Discovery 710 (64-slice PET/CT, General Electric Healthcare, Boston, Massachusetts, U.S), one Biograph (64-slice TruePoint PET/CT, Siemens Healthcare, Erlangen, Germany) and one Gemini (64-slice PET/CT, Philips Healthcare, Amsterdam, Netherlands). All PET/CT scans were obtained approximately 50 mins after injection of F-FDG (5 MBq/kg). All the acquired PET images were constructed using CT-based attenuation correction. Quality assurance was routinely performed to maintain consistent calibration on all PET/CT. In addition, all patients had carbohydrate antigen 19-9 (CA 19-9) values – the most extensively studied and validated physiological biomarker in PA prognosis [12] - measured pre-RT, and twenty-nine of them had CA 19-9 collected post-RT.

Data preprocessing

On attenuation corrected PET scans, standardized uptake values (SUV) were calculated and normalized to the patient’s body weight as:where Cimg(T) is the radioactivity concentration measured at time T, DA denotes administrated dose [MBq] at the time of injection, and BW represents patient’s body weight [kg] on the day of imaging. For simplicity, SUVBW will be shorted as SUV in the remaining text. To transform the PET indices into the CT space, SUV maps were up-sampled to the resolutions of respective CT images and used as the main source of metabolic inputs in later steps.

For each patient, boundaries of baseline and residual primary pancreatic gross tumor volumes (GTV) were manually defined by a single experienced gastrointestinal radiation oncologist on pre- and post-RT diagnostic CT images. To assist tumor definition, the contrast CT scans were registered to the corresponding PET/CT using a commercial B-spline based deformable image registration software (Velocity AITM Varian, Palo Alto, CA). The intra-observer reproducibility was evaluated with the same radiation oncologist contouring the tumor twice for four patients at two different time points eight weeks apart. The mean Jaccard index (intersection over union) was used to compare intra-observer variability. Pre- and post-RT CT images were rigidly registered to align the tumor orientation. The resultant transformation matrices were propagated to the corresponding PET images and tumor segmentation. This step enables following regional tumor surface matching. The propagated pre- and post-RT tumor segmentations (GTVpre and GTVpost) were binarized and assessed using the pipeline illustrated in Fig. 2.

Fig. 2

Illustration of our pipeline. Processing steps represented in numbers are: a) 3D construction; b) Surface-based registration; c) PET indices sampling; d) Anatomic parameter (detJ) calculation; e) Anatomic parameter (Radius) calculation; f) Multivariate surface vertex cloud assembling; g) Mahalanobis distance calculation on multivariate surface vertex ouds.

Surface measurement generation

As shown in Fig. 2, a 3D surface mesh was generated from GTV segmentation using the matching cube algorithm [13]. To obtain a regional correlation of pre- and post-RT tumor surfaces, GTVpre and GTVpost were matched using surface-based constrained harmonic registration [14], [15]. As such, shape changes were tracked within the structure-of-interest, excluding the influence of surrounding parenchyma [16]. For each surface vertex, a multivariate vector was formed using both the anatomic and metabolic components.

The anatomic component consisted of two geometrical features based on the same surface model: the normal direction, represented by the midline to surface radius , and the surface area, represented by the Jacobian of the transformation (J) from surface registration [16], [17]. Two Jacobian derived metrics were used: the determinant of Jacobian (detJ) and the log-Euclidean projected deformation tensor (log). detJ indicates the extent of warping needed to deform a surface to match its target, and log further encapsulates directional information in the deformation [15], [16].

To characterize the surface vertices using tumor metabolism, the midlines of all the 3D tumor volumes were first determined. In CT-aligned SUV maps, SUV measurements within the 3D tumor mesh were then sampled as surface vertices along the corresponding radius to midline, as illustrated in [19]. Specifically, using [18], [19]:and

where X , M, P denote the (x, y, z) coordinates of a vertex in the surface, the corresponding point of the vertex in the midline, and a voxel within the 3D representations, respectively. A pre-defined distance (D = 1.2) were chosen to ensure necessary voxel assignment while minimizing the overlap with neighboring vertices.

Statistical analysis

Paired Mahalanobis distance (Mdist), a dissimilarity measurement of observations in multivariate space, was used to assess tumor response. Specifically, for each tumor, given p-dimensional pre-RT surface observation , and the corresponding post-RT observation , the sample covariance matrix was calculated as:where n is the number of surface vertices, denotes the difference between and , and represents the arithmetic mean vector of over the surface vertices.

For each patient, Mdist of pre- and post- RT surfaces was obtained using:

Surface variables were analyzed accordingly: 1) metabolic features only - the surface projected SUVmax and SUVmean; 2) morphological features only – radius and log; 3) combined metabolic and morphological features – , detJ, and SUVmax. While log adds directional information in surface warping, it also increases the noise when combined with other metabolic and morphological features [16]. As a result, we chose detJ over log in the combined metrics. Mdists were calculated based on these three combinations and denoted as Mdist|SUV, Mdist|Shape, and Mdist|Fusion, respectively.

The prognostic values of the three Mdist metrics were evaluated using Cox proportional hazards regression with OS. OS was defined as the time span from the time of diagnosis to the time of patient death or last encounter. Correlation analyses were performed between OS and baseline (B) measurements of serum CA19-9, GTV, global SUVmax, global SUVmean, as well as post-RT versus pre-RT ratios (R) of the above four measurements. The enrolled subjects in our study were gender-matched (52% male and 48% female) but had a wide range of age. As a result, we also investigated the prognostic value of age as an OS predictor.

To evaluate the classification capacity of each risk predictor, all patients were stratified into two risk groups based on the three Mdist metrics. The low versus high-risk populations were determined according to 2 methods, M1 and M2. M1 first sorted each of the above metrics and then grouped the lower half subjects as a low-risk group (LG) and the upper half subjects as high-risk group (HG). M2 used k-means clustering algorithm to automatically separate LG and HG based on multivariate data, consisting of the tested metric and patient’s age at diagnosis. The inclusion of age in the analysis was again due to the wide range of age presented in this study and the significant association between age and OS, as revealed by our data. Log-rank tests were performed between LG and HG stratified patients using both M1 and M2, respectively. As a comparison, log-rank statistics of Mdist derived metrics were also compared with widely used imaging metrics (post- to pre-RT ratios of global SUVmax, global SUVmean, and GTV). Additionally, because the enrolled subjects were heterogeneous in terms of gender, resectability (BRPA versus LAPA), RT methods (SBRT versus CRT), and post-treatment resection (unresected versus resected), survival differences between subgroups were also examined by log-rank tests.

Results

Patient demographics and disease characteristics are listed in Table 1. The age of the enrolled subjects ranged from 45 to 84 years with a mean of 69.2 ± 9.5 (standard deviation) years. The mean Jaccard index of intra-observer GTV contour variability was 0.89. Results of univariate Cox regression analyses of all predictors are shown in Table 2. Among all the imaging parameters, Mdist|Fusion - the extent of tumor response considering both the morphologic and metabolic image information, was the best predictor of OS (p = 0.008) with a hazard ratio (HR) of 0.52 (95% CI, 0.19–0.78). Post-RT versus pre-RT ratio of SUVmax was associated with OS with borderline significance (HR = 3.89, 95% CI, 0.94–16.02, p = 0.060). All other imaging metrics (baseline or ratios of global SUVmean, GTV, as well as Mdist metrics that based on morphologic or metabolic information alone) were not associated with OS. Additionally, patient’s age at diagnosis was significantly associated with OS, (HR = 1.14, 95% CI, 1.04–1.25, p = 0.004), while traditional physiological biomarker CA19-9 did not present significant association with OS in our small cohort of patients.

Table 1

Patient characteristics.

Characteristics	Patients (%) (N = 31)
Gender
Male	16 (52%)
Female	15 (48%)
Age at Baseline Scan
Median	69
Mean	69.2
Std	9.5
Range	45–84
Clinical Tumor Stage
Stage I	1 (3%)
Stage II	1 (3%)
Stage III	7 (23%)
Stage IV	21 (68%)
Unknown	1 (3%)
NCCN1Resectability
Locally Advanced	23 (74%)
Borderline Resectable	8 (26%)
Primary Tumor Location
Uncinate, head, and neck	23 (74%)
Body and tail	4 (13%)
Combined	4 (13%)
Scan and Treatment Time Interval
Pretreatment (w)	4.54 (3.76)
Posttreatment (w)	7.14 (4.28)
RT Type2
SBRT	6 (19%)
CRT	25 (81%)
Resected3
Yes	11 (35%)
No	20 (65%)

National Comprehensive Cancer network.

RT types include stereotactic body radiation therapy (SBRT) and chemoradiotherapy (CRT).

This refers to the resection after RT.

Table 2

Results of univariate Cox proportional hazards regression between OS and 12 predictors. Four traditional prognostic measurements (serum CA19-9, global SUV, global SUV, and GTV) are categorized as baseline measurement (B), and post-RT versus pre-RT ratio (R). Surface measurements based on metabolic features only, anatomic features only and combined anatomic and metabolic features are indicated as M|, M|, M|, respectively. Statistical significance are marked as * for p-values < 0.05, and ** for p-values < 0.01.

	Mean	Std	HR/exp(β) (95% CI)	p values
Age	68.20	9.52	1.14 (1.04–1.25)	0.004**
CA19-9 (B) (k)	2.62	11.8	0.99 (0.99–1.00)	>0.1
CA19-9 (R) (n = 29)	0.68	0.66	1.63 (0.70–3.77)	>0.1
SUVmax (B)	6.00	2.37	0.90 (0.72–1.13)	>0.1
SUVmax (R)	0.81	0.39	3.89 (0.94–16.02)	0.060
SUVmean (B)	2.41	0.63	0.81 (0.32–2.04)	>0.1
SUVmean (R)	0.91	0.31	1.75 (0.45–6.75)	>0.1
GTV (B) (cc)	51.71	35.60	1.00 (1.00–1.00)	>0.1
GTV (R)	0.90	0.38	3.01 (0.57–15.71)	>0.1
Mdist|SUV	3.25	4.87	0.39 (0.24–1.13)	>0.1
Mdist|Shape	0.70	0.77	1.22 (0.69–2.14)	>0.1
Mdist|Fusion	2.68	4.62	0.52 (0.19–0.78)	0.008**

The capability of each measurement to classify patients is more evidently reflected in risk stratification analysis. Patients presented no survival differences between genders (female vs. male, p > 0.1), resectability (BRPA vs. LAPA, p > 0.1), treatment methods (SBRT vs. CRT, p > 0.1), and post-treatment surgery (unresected vs. resected, p > 0.1). Kaplan-Meier plots of risk groups stratified using M1 and M2 are shown in Fig. 3, respectively. Results of log-rank tests comparing each pair of the survival curves can be found in Table 3. Evidently, from Fig. 3, the survival risk of patients was best stratified by Mdist|Fusion, either alone or combined with patients’ age. Using sorted levels of each parameter (M1), Mdist|Fusion successfully classified an HG group with a mean of 0.57, and LG with a mean of 4.65, the former has a hazard ratio almost twelve times higher than the latter (p < 0.001). Mdist|SUV, the surface based distance measurement based on metabolic features only, also classified two groups of patients with significantly different risks (p = 0.021), while all other parameters failed in classifying survival risks using M1. Using unsupervised clustering (M2), a combination of Mdist|Fusion, and age significantly partitioned two risk groups (p < 0.001), with the HG having a nine-fold hazard risk than the LG. Mdist calculated on metabolic features alone (Mdist|SUV), and post- to pre-RT ratio of SUVmax also significantly classify HG and LG, with respective p values equal to 0.015, and 0.008. Other combinations were unable to stratify the patient cohort.

Fig. 3

Kaplan-Meier plots of OS for 31 patients stratified by post- versus pre-RT ratios of Global SUVmax, GTV, and the new M|, using M1, the evenly dividing method (a-c) and M2, the k-means clustering method (d-f). Censored observations are marked by black circles.

Table 3

Results of log-rank tests for patient risk groups stratified by 6 imaging based metrics, using M1 and M2. Surface measurements based on metabolic features only, anatomic features only and combined anatomic and metabolic features are indicated as M|, M|, M|, respectively. Log-rank statistics of groups stratified using M1 and M2 are displayed in white and gray rows, respectively. Significance are marked as * for p-values < 0.05, and ** for p-values < 0.01.

	LG		HG		HR (95% CI)	p values
	N	Mean (Std)	N	Mean (Std)
SUVmax (R)	16	0.51 (0.19)	15	1.12 (0.29)	3.11 (1.06–9.11)	0.073
	20	0.58 (0.22)	11	1.22 (0.27)	6.97 (1.94–25.07)	0.008**
SUVmean (R)	16	0.68 (0.12)	15	1.16 (0.27)	2.78 (0.93–8.27)	>0.1
	18	0.71 (0.14)	13	1.19 (0.27)	2.40 (0.80–7.24)	>0.1
GTV (R)	16	0.63 (0.23)	15	1.19 (0.28)	2.78 (0.88–8.78)	>0.1
	13	0.57 (0.22)	18	1.14 (0.28)	1.47 (0.50–4.34)	>0.1
Mdist|SUV	16	5.33 (6.13)	15	1.03 (0.44)	4.30 (1.42–13.02)	0.021*
	15	5.50 (6.32)	16	1.15 (0.54)	4.89 (1.56–15.29)	0.015*
Mdist|Shape	16	1.18 (0.81)	15	0.17 (0.10)	1.90 (0.66–5.52)	>0.1
	20	0.24 (0.15)	11	1.52 (0.76)	1.04 (0.35–3.08)	>0.1
Mdist|Fusion	16	4.65 (5.82)	15	0.57 (0.35)	11.51 (3.48–38.02)	<0.001**
	13	5.37 (6.27)	18	0.73 (0.50)	9.04 (2.89–28.33)	<0.001**

Discussion

Accurate response assessment during and after radiotherapy for LAPA and BRPA is critical to clinical decision making as it may help identify radioresistant versus radiosensitive tumors and influence the decision for definitive resection. Although single dimension tumor measurements using abdominal CT remains the primary criteria for assessment of treatment response, a number of studies have shown the synergistic role of PET to abdominal CT in response evaluation. As with other neoplasms, FDG uptake was found to significantly differentiate responders and non-responders and predict survival in patients with PA [8], [20], [21]. While the FDG uptake promisingly depicts tumor malignancy, rudimentary analyses, typically global SUVmax, and SUVmean, have limited accuracies in assessing the biologically complex PA tumors (Fig. 1). As demonstrated by genomic sequencing analysis, multiple tissue samples taken from the same tumor displayed a wide variety of cytogenetic differences [22], [23]. Therefore, quantitative imaging methods describing the intra-tumor heterogeneity are needed.

In the past decade, fused PET/CT has become the standard imaging modality localizing FDG-intense foci. At the same time, more sophisticated PET/CT analysis algorithms have been developed to describe intra-tumor spatial heterogeneity. In particular, Tixier and colleagues calculated textural features, such as homogeneity and entropy, on baseline esophageal tumors and successfully predicted response to RT in 41 patients [24]. Yue and colleagues further compared pre- and post-RT intra-tumor PET textural features in 26 PA patients, and identified 3 borderline-significant prognostic variables, including homogeneity, variance and cluster tendency [25]. Instead of focusing only on FDG-avid regions, textural features describe the distribution of FDG uptake in the entire tumor region, which improved the characterization of tumor metabolic heterogeneity. However, in these PET textural studies, tumor volumes were only served as a binary mask for textural feature extraction, without contributing morphological information into the analyses. New alternatives that integrate both anatomic and metabolic information of tumor may provide a more complete picture of tumor response.

Tumors exhibit heterogeneity not only in space but also across time. In one hand, some PA cells grow faster than others, for instance, tumor cells in close proximity to pancreatic stellate cells have accelerated proliferation [26]. On the other, certain PA tissues, such as tumorigenic cancer cells defined by the CD133 expression, are more treatment-resistant than others [27]. Collectively, targeting the tumor microenvironment variations and probing into spatial-temporal coincidence in longitudinal imaging data may provide valuable insight of tumor treatment response [25], [28]. In this study, we developed a 3D tumor surface based pipeline combing both the metabolic and morphological features to quantify tumor variations between pre- and post-RT. The tumor surface model characterizes the 3D tumor with regional spatiotemporal features. Comparing to univariate shape descriptors over the entire GTV, the use of 3D surface model facilitates comparisons of regional morphological features on locally aligned surface vertices, thus are sensitive to small shape alterations from tumor responses [13], [14]. By projecting the SUV onto the tumor surface, regional metabolic features were integrated into the 3D model, and the intra-tumor metabolic heterogeneity was largely preserved in the response evaluation.

Our results showed that the combined morphological and metabolic metric, Mdist|Fusion, significantly predicted OS and stratified risks for our PA patients, and outperformed traditional CA19-9, global SUV, and GTV based measurements. Furthermore, using the same 3D surface model, the morphological and metabolic fused metric Mdist|Fusion were superior to morphological or metabolic features alone (Mdist|Shape, and Mdist|SUV) in OS predicting and risk stratification. This, in turn, indicates the more accurate characterization of tumor response using the combined model. In addition to the PA patients, the presented pipeline may be applicable to other solid tumors for response evaluation.

There are several limitations in this study. First, in this pilot study, the enrolled participants were retrospectively solicited from our radiation oncology practice. The small sample size limits the power of statistical tests, and preclude the possibility to reliably remove the influence of all confounding factors. The comparisons between new multimodal based metric to traditional imaging metrics were made in the same condition and thus less biased by the heterogeneity of patient characteristics (tumor stage, radiotherapy fraction, censorship, etc.). However, the use of multiple PET/CT scanners may have led to some inconsistency in metabolic metric calculations despite our best effort to calibrate these scanners to be consistent. A calibration comparison study on 11 regularly maintained PET/CT scanners observed good reproducibility and accuracy of all dose calibrators and reported no major deviation across different PET calibration systems [29]. Nevertheless, many quantitative descriptors draw from GTVs, for instance, textural features, have been reported to have limited robustness against multi-scanners [30]. Future larger cohorts preferably scanned in the same PET/CT scanner with high-resolution PET protocols, are needed to validate our method and findings. Additionally, segmentation accuracy is paramount to quantitative imaging analysis. A previous study investigating the reproducibility of PET radiomic features against GTV delineation has reported increased robustness of derived radiomic descriptors when using semi-automated segmentation over manual slice-by-slice delineation [31]. Same concerns also apply to the presented method, as it is sensitive to tumor volume segmentation errors that are prominent with less conspicuous tumors. With recent developments of deep learning, several automated tumor segmentation studies have shown promising results [32], [33]. While accurate pancreatic tumor delineation leads to ongoing research, semi-automated or automated segmentation tools may be integrated into our pipeline to improve consistency and accuracy.

In conclusion, based on PET/CT, we introduced a novel combined morphologic and metabolic pipeline for treatment assessment of locally advanced or borderline resectable PA. The fused Mdist outperformed traditional morphologic, metabolic, and physiological measurements in OS prediction and risk stratification. The presented fused model shows the promise to be a new imaging-marker capable to improve characterization of tumor response and personalized PA care. Future large-scale studies will be needed to validate the capacity of the combined metric generated here.

Conflicts of interest

All authors confirm that there are no known conflicts of interest associated with this publication.

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.