2.1.

Structural Modeling and Metabolic Gradients

Here, we describe the algorithmic approach used to construct a semiparametric regression profile of the volumetric uptake data in terms of a reference ellipsoidal pattern. Hereafter, the $N$ observations $(\mathbf{x},Y)$ are described by their 3-D coordinates $\mathbf{x}\in \mathrm{\Omega }$ in the image domain $\mathrm{\Omega }\subset {\mathbb{R}}^3$ and measured uptake value $Y$. The spatial distribution of uptake observations is characterized in terms of its compliance to a rigid 3-D ellipsoidal model parametrized by $\theta =(\mu ,\mathrm{\Sigma })$ for shape $\mathrm{\Sigma }$ (the uptake data covariance matrix) and location $\mu$. Voxel location within this model can be expressed in terms of their radial position within the ellipsoid:

This provides the opportunity to represent the volumetric uptake information as a function of ellipsoidal radius, as shown in the profile plot in Fig. 1. A nonparametric (either isotonic or bitonic) regression $f$ of this profile $(u,Y)$ can be obtained given a set choice of ellipsoidal parameters $\theta$ in the form of a model:

Fig. 1

FDG-PET uptake profiles of two sarcoma studies: case A (49-year-old male with upper thigh soft tissue sarcoma, alive at last follow-up 4 years after baseline scan) presented with an active, homogeneous core (top) and case B (48-year-old male with pelvis soft tissue sarcoma, died 8 months after baseline scan) presented with a heterogeneous core with low activity (bottom). These features are visible from the transverse and coronal field-of-view images (left). Centre plots: corresponding uptake profile $[u(\widehat{\theta }),Y]$, initial stepwise monotonically decreasing fit $f[u(\widehat{\theta })]$ (inner navy curve), stepwise unimodal regression $\widehat{f}[u(\widehat{\theta })]$ (thinner, red curve), and smoothed bitonic uptake profile curve $\tilde{f}[u(\widehat{\theta })]$ (thicker, blue curve). Initial stepwise fit $f(u)$ is barely noticeable for case A but highlights the benefit of using a unimodal fit in case B. Rightmost plots summarize gradient information: top plot illustrates the gradient curve $g(u)$ Eq. (4), and uptake weighted gradients $g_Y(u)$, as functions of radial location $u$. Reference voxels are identified to illustrate decreasing and increasing uptake rates (a blue square and a red circle, respectively). The gradient lengths at those locations correspond to the absolute value of the normalized gradients $\tilde{g}$ Eq. (6). The plot for case B shows how the gradient measure has a negative sign for voxels located on the left hand side of the uptake profile mode (i.e., voxels closer to the core with a decreasing rate of uptake). Bottom plot: histogram of normalized gradients $\tilde{g}$ on a common binning scale. Broken vertical lines indicate the 25th and 95th quantiles of normalized gradients.

The original approach² consisted in fitting a stepwise isotonic decreasing nonparametric least squares regression function $f$ to the uptake data $Y$, following the idealized representation of gradually decreasing ellipsoidal level curves at increasing ellipsoidal radii. A measure of lack-of-fit of this template structure was used directly to quantify heterogeneity (this is described in more details further). Here, we used this approach to quantify heterogeneity but introduced a second, bitonic (i.e., unimodal) regression for $f$ in order to perform further assessment of the structure of the 3-D uptake distribution. This unimodal level profile $f$ provides a more reasonable description of a tumor that might have a central necrotic core (see Fig. 1). (This principle also applies to other tumor types for which distributions have an overall decreasing radial uptake pattern, as it has the ability to differentiate sharper concentrations of uptake.) An algorithm to implement the unimodal fit is described in Appendix A. Given the fitted nonparametric unimodal uptake profile, it is now possible to evaluate associated gradients $-f^{\prime }[u(\theta )]$ at any point within the tumor. These signed uptake curve gradients can be used as a variable for tumor characterization (and prognostic assessment, as discussed hereafter). As the unimodal regression curve $\widehat{f}$ obtained for Eq. (2) is a stepwise function, a smoothing spline is subsequently applied to it in order to obtain a final unimodal continuous uptake profile function $\tilde{f}$, i.e.,

This creates the opportunity to define a variable for tumor characterization in the form of the signed metabolic gradient:

2.2.

Structural Variables and Interpretation

A number of structural variables can be extracted from the modeling approach described in Sec. 2.1 for tumor characterization. As in previous works from the group,^2,3 intratumoral heterogeneity may be defined as a measure of lack-of-fit of the spatial model to the observed uptake distribution. Two versions of such a heterogeneity variable can be constructed as follows using Eq. (3):

Various gradient summaries may be obtained from Eq. (4) to further describe intratumoral status and activity. Normalized gradients may be derived from Eq. (4) to conform to a universal scale by

As each voxel $i$ receives a gradient value $g_i=g[u_i(\widehat{\theta })]$, the sample of $N$ signed gradients may be analyzed, for example, in terms of its first quartile or its 95th percentile, to capture summarizing features of the rate of metabolic change within the VoI. In this view, taking the first quartile of (normalized) gradients provides a single-valued evaluation at the lower end of the scale of metabolic rates: one would expect a lower first quartile for a tumor with decreasing metabolic activity, for example, resulting from necrosis. Likewise, the 95th percentile of (normalized) gradients provides an analysis of higher rates, where a higher value would likely correspond to a more rapidly changing tumor metabolism. Figure 1 illustrates this idea.

The sample of metabolic gradients may also be weighted by voxel uptake into

2.3.

Textural Quantitation

First-order, second-order, and regional features commonly found in other contributions^1,6,22–24 were computed as per definitions provided by the Image Biomarker Standardization Initiative²⁵ (version 1.5). Second-order features are defined as moments and functions of the grayscale level co-occurrence matrix (GLCM), which is made of the frequencies $P_{ij}$ of adjacency of any two uptake values $(i,j)$ in a given direction.¹² The output texture features are computed and averaged over all 13 directions in the volume. We use the normalized GLCM, i.e., $\sum _{i,j}{P}_{ij}=1$, where matrix $P$ is symmetric. Regional features were derived from a gray level size zone matrix (GLSZM), which evaluates the numbers and sizes of contiguous homogeneous regions of equal (discretized) gray level.^22,24 Abundant descriptions of the above-mentioned features can be found in the literature cited above. The list of features considered here is by no means comprehensive but includes the main first- and second-order features found in the literature on PET-derived radiomics. The two regional features were selected on the basis of their relevance in other works,^6,22,24 their relative algorithmic simplicity, and to introduce higher-order evaluations in our analysis.

Other variables commonly found in the literature in the context of radiomics include assessment of shape²⁶ or morphology.^24,25 We also included in our analyses an evaluation of volume asphericity²⁶ proportionally to the ratio of the segmented boundary surface $S$ and its volume $V$ by $S/{(36V\pi )}^{1/3}-1$. A comparison of morphological ellipsoidal features²⁵ with our model-derived features was aptly suggested by a reviewer. On this basis, five morphological descriptors for ellipsoidal characteristics²⁵ were included that rely on volume eigenvalues ${\lambda }_{\text{major}}\ge {\lambda }_{\text{minor}}\ge {\lambda }_{\text{least}}$ derived by principal component analysis (PCA). These features assess major, minor, and least axis lengths (defined, respectively, as $4\sqrt{{\lambda }_{\text{major}}},4\sqrt{{\lambda }_{\text{minor}}}$, and $4\sqrt{{\lambda }_{\text{least}}}$) as well as volume elongation and flatness ($\sqrt{{\lambda }_{\text{minor}}/{\lambda }_{\text{major}}}$ and $\sqrt{{\lambda }_{\text{least}}/{\lambda }_{\text{major}}}$, respectively).

3.1.

Dataset and Analysis Framework

The dataset of primary sarcoma tumors was acquired at the University of Washington in Seattle, United States, between August 1993 and January 2003, after patients were diagnosed by biopsy. After exclusion of unsuitable cases due to lack of complete data availability, the cohort ($N=197$, 88 deaths observed) is composed of 130 soft tissue, 51 bone, and 16 cartilage sarcomas, in patients aged between 17 and 86 years of age (median 45), of which 86 females and 111 males, with 99 high-grade, 66 intermediate, and 32 low-grade tumors.

All patients underwent standard [18F] FDG-PET on a GE Advance PET scanner before either neoadjuvant chemotherapy or surgical resection. Patients received an intravenous injection of [18F] FDG (259-370 MBq) after fasting for at least 12 h and rested for between 45 and 60 min before emission and attenuation scan acquisition of the tumor fields of view (FoVs). All imaging data were reconstructed using filtered backprojection (FBP), rendering 3-D images using a Hanning filter after scatter correction, resulting in a reconstructed resolution of 10 mm.²⁷ Output image sets had voxel size of $4.30\mathrm{mm}\times 4.30\mathrm{mm}$ in the transverse plane and slice thicknesses of 4.25 mm. Details on the imaging protocol were reported in previous reports.^2,28 Raw counts were scaled into SUV with respect to activity in injected dose per unit weight of the patient ($\mathrm{kBq}/\mathrm{g}$).

Primary tumors were identified by a radiologist. Crude input ellipsoidal VoIs were drawn around the entire FDG-PET tumor volume in AMIDE²⁹ for volumetric (3-D) statistical analysis. Our model-based quantitation approach is robust to the choice of input VoI, and in fact, the analysis can be performed on the crude bounding box (or bounding ellipsoid),^2,3 which could include a (relatively small) proportion of background voxels. However, texture features do require volume segmentation, which was performed in three different ways in our initial analyses. Here, we present results obtained for a fixed-threshold segmentation similarly, e.g., to Yan et al.,³⁰ where the threshold value is set for each study based on the subsample of the lower 15% of uptake values (so as to include background and healthy tissue activity only). For a given study, the segmentation threshold is defined as the mean subsample value plus three standard deviations of this subsample. This approach yielded satisfactory delineation throughout the cohort. Alternative segmentation techniques considered consisted in using either the convex hull of this threshold-segmented volume or the volume segmented obtained by localized fitting of a tubular representation.³ Results for these alternative approaches are left out here, as the analysis of their impact is outside the scope of this paper and may be considered in future works. The general conclusions of this paper remained unchanged by the choice of segmentation strategy.

Given the near-homogeneous voxel dimensions, no interpolation was performed prior to VoI resegmentation for texture analyses. Thus, both structural and textural groups of features were obtained from the same segmented VoI. Uptake values $Y$ were requantized into $Q=32$ gray levels $Y^Q$ by fixed bin number transformation:

Experimentation indicated that although texture quantitation changed with the choice of $Q$, this did not meaningfully impact general conclusions from our survival analyses. We have not considered the impact of alternative quantization techniques; this question is beyond the scope of this paper.

A total of 44 variables were considered and may be identified in three frames as follows: (i) routine clinical variables [tumor grade, clinical volume, clinical tumor subtype, patient age, patient sex, maximum standardized uptake value (SUVmax), mean uptake value (SUVmean) and total lesion glycolysis (TLG) were collected for this cohort]; (ii) structural features including heterogeneity using both ${\mathcal{H}}_0$ and ${\mathcal{H}}_1$, and the raw, normalized, and uptake-weighted gradients (Sec. 2.2); and (iii) a set of image summaries including morphologic and texture features (Sec. 2.3).

The primary endpoint for all survival analyses was overall survival status at last follow-up. The sarcoma cohort comprises tumors of three different grades and subtypes described above. In particular, it includes a high number of high grade soft tissue sarcomas (STS). For this reason, we excluded a total of 15 patients, who were lost to follow-up within the first 2 years for survival analyses, a common horizon in survival studies of STS.³¹ This threshold seems a reasonable compromise as median survival times can drop to 11 months for higher stages, irrespective of grade, according to the TNM classification of malignant tumours (TNM).^31,32

Univariate Cox survival analyses performed for all of the variables revealed the following were significant, at least at the 5% significance level: (clinical) tumor grade, subtype, tumor volume, SUVmax, SUVmean, and TLG; (structural) ${\mathcal{H}}_0$, ${\mathcal{H}}_1$, ${\tilde{g}}_{[0.25]}$, $g_{Y,[0.95]}$; (textural) asphericity, ${\mathrm{PCA}}_{\mathrm{major}}$, ${\mathrm{PCA}}_{\mathrm{minor}}$, ${\mathrm{PCA}}_{\mathrm{least}}$, ${\mathrm{mean}}_{\mathrm{HIST}}$, ${\text{skewness}}_{\mathrm{HIST}}$, ${\text{kurtosis}}_{\mathrm{HIST}}$, ${\text{median}}_{\mathrm{HIST}}$, ${\text{entropy}}_{\mathrm{HIST}}$, ${\text{uniformity}}_{\mathrm{HIST}}$, $\min .{\mathrm{gradient}}_{\mathrm{HIST}}$, ${\text{entropy}}_{\mathrm{GLCM}}$, ${\text{dissimilarity}}_{\mathrm{GLCM}}$, ${\text{homogeneity}}_{\mathrm{GLCM}}$, ${\text{contrast}}_{\mathrm{GLCM}}$, ${\text{autocorrelation}}_{\mathrm{GLCM}}$, and ${\max .\mathrm{probability}}_{\mathrm{GLCM}}$.

3.2.

Feature Space Exploration

3.2.1.

Correlation

The correlation matrix of the set of variables provided in Fig. 2 demonstrates the strong separation between structural and textural variables, with a few noticeable exceptions: coefficient of variation of the histogram (${\mathrm{CoV}}_{\mathrm{HIST}}$) and max. ${\text{probability}}_{\mathrm{GLCM}}$ aligned with $g_{[0.95]}$ and ${\tilde{g}}_{[0.95]}$, and ${\text{entropy}}_{\mathrm{GLCM}}$ and ${\text{homogeneity}}_{\mathrm{GLCM}}$ with ${\mathcal{H}}_0$.

Fig. 2

Correlation matrix of the continuous variables used in the survival analyses. All correlations within $(-0.65,+0.65)$ have been reduced to the predominant gray level in the image. Gray levels darker and brighter than that level indicates, respectively, negative and positive correlations greater than 65%.

3.3.

Feature Selection

The above feature space exploration indicates that structural and radiomic variables tend to span different areas of the information space. This suggests that (i) structural variables derived from spatial modeling of the volumetric FDG uptake distribution have prognostic potential in the analysis of sarcomas and (ii) they remain useful when used in combination with texture features commonly used in radiomics. The availability of many quantitative variables for the characterization of tumor metabolism raises the question of final selection of prognostic variables. Multivariate prognostic models may be defined directly from available features (for instance, using only structural features, only radiomic features, or a combination of both), or on the basis of a feature selection approach. We considered several typical alternatives for this purpose, exploring two avenues: using principal components analysis and clustering; and using typical machine learning classifiers. These methodologies are described hereafter.

4.1.

Feature Selection from Outcome-Based Classifiers

4.1.1.

Outcome-based classification performance

Figure 3 shows the cross-validated ROCs obtained by averaging ROCs across all $M=1000$ cross-validation samples for all classifiers (except the single-layer neural network, which did not perform better than the two-layer network and was left out of further analyses). The corresponding cross-validated AUCs indicate that the LASSO, forward-stepwise logistic regression, and random forest models yield close, acceptable performance, and outperform the other classifiers. Interestingly, different classifiers tend to select very different feature sets, as we see next.

Fig. 3

Receiver operating characteristic curves for the classifiers for outcome-guided feature selection (SVM: support vector machines; F/B-logistic: forward- and backward stepwise logistic regression; NN: neural networks; RF: random forest). AUCs are specified in increasing order (inset). Random forests, forward-stepwise logistic regression, and the LASSO yielded effective discrimination rates.

4.1.2.

Feature selection rates per classifier

The top plot of Fig. 4 summarizes cross-validated variable selection rates for each classifier. This method-specific summary highlights features that were selected for over 50% of the $M=1000$ cross-validation samples–a darker box indicating a higher selection rate. Rates are shown for “smaller” ($k=6$) and “larger” ($k=16$) models in order to emphasize how feature selection evolved with $k$. Surprisingly, the random forests and the neural networks did not elect grade consistently as one of the most predominant features, unlike the more traditional methods (stepwise logistic regression, LASSO, and simulated annealing). Some structural variables, and in particular, ${\mathcal{H}}_0$ and ${\mathcal{H}}_1$, are selected rather consistently. Random forests and neural networks noticeably favor a selection of structural features over grade (the former may capture grade-related aspects of the tumor characteristics²). For radiomic features, the selection process is less systematic, in that variables are elected either less often or in different combinations for varying cross-validation samples.

Fig. 4

(a) Feature selection rates per classifier. F/B-log: forward- and backward stepwise logistic regression; SA-k: simulated annealing; RF-k: random forest; NN-k: neural networks; where “-k” denotes using the $k$ most important features on average from that model (hence considering 6-variable and 16-variable models). Selection rates below 50% were coded in white; rates above this threshold are represented with a proportionally darker color. The two vertical lines are supplied as visual guides to identify the routine (left-most), structural (central), and textural (right-most) subsets of variables. (b) Feature selection rates overall, i.e., combined across all $M=1000$ cross-validated experiments, for the same subset of six classifiers comprising of both logistic classifiers, the LASSO, and simulated annealing, random forest and neural network classifiers with $k=16$. Features are ranked in decreasing order of popularity and color-coded according to feature frame.

Table 2 presents the multivariate prognostic analyses based on the most popular models obtained from five selection schemes: forward-stepwise logistic regression selection, simulated annealing, the LASSO, random forests, and neural networks. For each method, the $k=15$ most frequent variables were identified. Of these subsets, those variables selected for over 50% of Monte Carlo repetitions were used together into a multivariate model. The table indicates that clinical and structural features are selected consistently by all methods and that at least one gradient summary is among these variables deemed most important by the selection scheme.

Table 2

Cox proportional hazard model analyses of multivariate prognostic models comprising of only clinical and structural variables, only clinical and radiomics variables, combining both, or obtained from forward logistic selection, simulated annealing, the LASSO, random forests and neural networks with k=15 variables. In each model, the variables selected by the corresponding technique for over 50% of the M=1000 cross-validation sets were used. This cut-off selection rate was chosen arbitrarily but different values yield to similar conclusions. This table indicates both the most frequently selected features for each model, and contributions that are significant at the 5% significance level (in bold).

	Model [HR (p)]
Variable	Structural	Radiomics	Combined	F-logistic	SA	LASSO	RF	NN
Grade (int)	0.46 (0.0020)	0.44 (0.0009)	0.48 (0.0034)	$-\mathbf{0.47}$(0.0029)	0.46 (0.0032)	0.47 (0.0061)	—	—
Grade (high)	0.21 (0.0035)	0.20 (0.0022)	0.21 (0.0030)	$-\mathbf{0.21}$(0.0030)	0.19 (0.0022)	0.21 (0.0047)	—	—
Sex (male)	—	—	—	—	1.45 (0.1116)	—	—	1.34 (0.2177)
Subtype (cart)	—	—	—	—	—	0.67 (0.5541)	—	—
Subtype (STS)	—	—	—	—	—	1.14 (0.6119)	—	—
SUVmean	—	—	—	—	—	1.24 (.0259)	1.63 (0.0105)	1.24 (0.2011)
SUVmax	1.13 (0.2658)	1.56 (0.0007)	1.47 (0.0062)	—	—	—	0.81 (0.5959)	—
TLG	—	—	—	—	—	0.99 (0.9012)	1.09 (0.6597)	1.00 (0.9806)
${\mathcal{H}}_0$	—	—	—	1.42 (0.0001)	1.23 (0.0871)	1.24 (0.0818)	1.54 (0.0166)	0.96 (0.8053)
${\mathcal{H}}_1$	1.89 (${<\mathbf{10}}^{-\mathbf{5}}$)	—	1.78 (0.0017)	—	1.35 (0.0217)	1.36 (0.0139)	1.07 (0.7337)	2.32 (0.0004)
${\tilde{g}}_{[0.25]}$	—	—	—	0.83 (0.0438)	0.68 (0.0026)	0.83 (0.0744)	0.93 (0.5722)	—
$g_{Y,[0.25]}$	—	—	—	—	1.41 (0.0077)	—	—	1.13 (0.4705)
${\tilde{g}}_{[0.95]}$	1.46 (0.0101)	—	1.99 (0.0009)	—	—	—	—	1.94 (0.0072)
$g_{Y,[0.95]}$	—	—	—	—	—	—	0.88 (0.6935)	—
${\mathrm{PCA}}_{\text{minor}}$	—	—	—	—	—	—	1.02 (0.9327)	—
${\mathrm{PCA}}_{\mathrm{elongation}}$	—	—	—	—	—	—	0.96 (0.7048)	—
Max. grad.${\mathrm{Gray}}_{\mathrm{HIST}}$	—	—	—	—	—	—	—	1.16 (0.1955)
${\mathrm{CoV}}_{\mathrm{HIST}}$	—	0.42 (0.0002)	0.44 (0.0053)	—	—	—	—	—
${\mathrm{Skewness}}_{\mathrm{HIST}}$	—	—	—	—	—	—	—	1.00 (0.9961)
${\mathrm{Entropy}}_{\mathrm{HIST}}$	—	0.03 (0.0003)	0.12 (0.0443)	—	—	—	0.45 (0.2353)	—
${\mathrm{Contrast}}_{\mathrm{GLCM}}$	—	—	—	—	0.85 (0.1851)	—	—	—
${\mathrm{Entropy}}_{\mathrm{GLCM}}$	—	84.59 (0.0019)	15.56 (0.0603)	—	—	—	1.83 (0.3316)	—
${\mathrm{Homogeneity}}_{\mathrm{GLCM}}$	—	6.35 (0.0074)	2.79 (0.1274)	—	—	—	—	—
${\text{Intensity}}_{\mathrm{GLSZM}}$	—	—	—	—	0.95 (0.6064)	—	—	—
Concordance	0.72	0.72	0.74	0.70	0.73	0.73	0.69	0.70

There was more variability in the selection of textural features (partly due to the textural frame being larger than the other two feature frames). These lower selection rates do not indicate that this quantitation strategy is less effective for prognosis (Fig. 5, described further, indicates the contrary), but rather that several groupings of these features seem to hold comparable prognostic potential. In other words, different subgroups are in turn important for different cross-validation samples of the cohort. This suggests that the information captured by these conventional radiomic features may be more versatile.

Fig. 5

Rates of selection of either clinical-only models, structural-only models, radiomics-only models, or models that combine (at least) structural and radiomic features, out of all cross-validation samples and from the logistic models, the LASSO and six-variable simulated annealing, random forest and neural network classifiers considered.

4.2.

Outcomes Analyses

4.2.1.

Outcomes analyses from conventional models

Table 2 includes structural-only and radiomics-only prognostic models to explicitly demonstrate the prognostic potential of either feature frame in sarcoma. These two models were selected as they include the same clinical variables, have optimal AIC, and achieve high concordance compared to other structural-only and radiomics-only models. This analysis establishes the significance of texture analysis for prognostic modeling of overall survival (Vallières et al.²⁰ showed the predictive potential of texture analysis in a joint PET-MRI predictive model for the presence of lung metastasis in sarcoma). Similarly, structural variables ${\mathcal{H}}_1$ and ${\tilde{g}}_{[0.95]}$ are prognostic when used instead of the radiomic features.

Figure 6 shows Kaplan–Meier analyses obtained from various multivariate risk models. Each analysis stratifies the sarcoma cohort into lower- and higher-risk groups on the basis of maximum log-rank test statistic so as to maximize the dichotomic separation. In all cases, separation was statistically significant ($p<{10}^{-7}$). The left-hand-side panel shows output from three traditional models: the structural model of Table 2 (solid lines), the prognostic model comprising of the PCA components of Sec. 3.2.2 (dashed lines), and the $k$-means centroid model Eq. (8), using the 15 cluster centroids as risk predictors (solid lines). Although risk stratification varies with the model, risk separation remains comparable.

Fig. 6

Kaplan–Meier survival curve estimates for dichotomic (low- versus high-) risk group segregation on the basis of various feature selections. For each model, the corresponding two Kaplan–Meier curves were obtained by maximizing the risk group separation in terms of highest log-rank test statistic via grid search. (a) Analyses using models derived from PCA model, $k$-means centroids Eq. (8), and AIC-based selection of structural variables (out of ${\mathcal{H}}_0$, ${\mathcal{H}}_1$, and the gradient summaries). (b) Similar Kaplan–Meier analyses for models obtained from machine learning techniques (LASSO, random forests, and neural networks). All models, although different, yield comparable survival analysis performances and significant risk stratification (log-rank test $p<{10}^{-7}$).

In the multivariate Cox proportional hazard analysis for the $k$-means centroid model Eq. (8) used in Fig. 6, significant cluster centroids were those of clusters ${\mathcal{C}}_1$, ${\mathcal{C}}_7$, ${\mathcal{C}}_8$, and ${\mathcal{C}}_9$ ($p<{10}^{-4}$, $p=0.0369$, $p=0.0129$, and $p=0.0096$, respectively), which are all summaries of routine clinical and structural assessment (cf. Table 1). In this model, the grouping of a relatively small number of structural features captures critical prognostic information for the cohort, unlike the summaries of conventional radiomic features.

5.1.

Summary of Findings

Building on the results of recent spatial modeling approaches for the analysis of sarcoma,^2,3 this paper presents a variable, namely a measure of metabolic gradient, for further characterization of the metabolic activity in sarcoma tumors based on FDG-PET imaging data. A previous uptake profiling approach² has been adapted to allow for a coherent calculation of gradients of the metabolic uptake profile obtained from a spatial model of the tracer uptake distribution. The profile and gradients are evaluated at each voxel and can therefore be mapped onto the imaged volume of interest for visualization. This sample of gradients can also be summarized into a single value (for example, its 95th percentile), to provide opportunities for multivariate prognostic evaluation along with other features derived from the structural representation of the VoI. Here, we used this information in three forms by considering the raw, uptake-normalized, and uptake-weighted gradients. Each of these summaries has its own practical interpretation, but overall they all evaluate local rates of change in glycolysis.

This paper also demonstrates that the proposed model-based structural variables (heterogeneity and metabolic gradients) complement characteristics captured by radiomic features, in that each quantitation methodology captures distinct parts of the information space. Ellipsoidal morphologic features were especially of interest given they are based on the same reference shape as is our model. In previous analyses (e.g., by O’Sullivan et al.⁵⁴), morphology was not found to be an independent risk factor for sarcoma prognosis. Here, we found no correlation between our metrics and those morphologic features, which rather capture volumetric aspects. Similarly, no relationships were found between metabolic gradient summaries and conventional histogram gradient metrics.

Different mechanisms were considered for feature selection, including PCA and machine learning techniques, which are particularly suited to the PET imaging analysis context, where the number of covariates is likely to be (much) greater than the size of the cohort. In total, 37 schemes were considered, including the LASSO, random forests, and neural network models. It was seen that feature selection varied with the technique considered. Overall, structural variables assessing heterogeneity and metabolic gradients were rated among the most important prognostic variables across all methods, along with routine variables, such as tumor grade and SUVmean. The importance of these features may vary with the disease and imaging modality but appears to be relevant in the case of FDG-PET imaging of sarcoma.

In summary, the contribution of this paper is threefold. (i) It defines a promising approach for evaluating the tumor metabolic gradient based on FDG-PET imaging data and validates its statistical significance in prognostic modeling of overall patient survival. Through this process, the paper also demonstrates the (separate) prognostic utility of typical texture features in sarcoma. (ii) It demonstrates that descriptors derived from our spatial model of PET tracer uptake capture aspects of patient information that complement those described by morphologic and textural features. (iii) It illustrates that structural and conventional radiomic features can be successfully combined for multivariate prognostic modeling from a range of statistical learning techniques.

5.2.

Discussion

Structural variables allow for both quantitative and qualitative assessment due to interpretability of their underlying model and could be used in various ways in clinical settings. For example, image mapping of the metabolic rates could be used to guide understanding of the metabolic process, e.g., in predicting areas in the volume more likely to see an increase in glycolysis at short term. More generally, both structural and textural features are typically used to measure (some form of) intratumoral heterogeneity, and the literature often understands or summarizes these quantitations without explicit distinction in many contributions.^1,2,5,17 Some descriptions of these approaches acknowledge subtle differences in what they capture, by referring to, e.g., spatial or textural heterogeneity.^1,2,5 It is usually accepted that PET imaging allows one to assess macroscopic heterogeneity, which is not clearly linked to microscopic heterogeneity in the literature. Recent works suggest that only a small number of texture features may actually correlate with microscopic assessment (notably some work on mice⁵⁵), but no report clearly establishes a correlation with other macroscopic assessments. On the other hand, structural analyses, such as the ellipsoidal-based template approach considered here, have links to histologic grading of sarcomas² (spatial modeling of the tracer uptake distribution also has the ability to capture other characteristics of the disease). These two methodologies may in fact capture very different aspects of the tumor histologic profile. Simple examples can highlight that texture features assess tumor characteristics that are not structural in nature; a toy example presented in Appendix B illustrates this.

Incorporation of subtype into the analysis was not found significant in any of the prognostic models present in this work–this is consistent with the result reported by Eary et al.² The observed complementarity between structural and radiomic variables may merit more detailed examination of sarcoma subtypes information. Further exploration of this issue may be warranted.

The structural methodology proposed here is segmentation-insensitive by nature, as technically it does not rely critically on the whole sample of observations. For example, lower extreme values can be removed or down-weigthed for model-fitting without impacting the output analysis. This is not true for textural analyses, where any processing of the input sample will directly affect the GLCM and other critical features that the quantitation may rely on.

The clinical dataset analyzed here is attractive by its size (sarcoma having relatively low prevalence) and the breadth of patient follow-up information. It was, however, acquired using an old-generation tomograph using FBP reconstruction. Many publications have highlighted the significant and varying influence of image reconstruction settings on texture features in 18F-FDG PET^30,56–58 and the question also arises for our model-based quantitation. A simulation study presented in Appendix C illustrates the relative stability of model-based quantitations with respect to the choice of PET reconstruction technique applied to a 2-D simulated phantom. It shows that for both FBP and maximum likelihood (ML) reconstruction techniques, error characteristics for ${\mathcal{H}}_0$, ${\mathcal{H}}_1$, and ${\tilde{g}}_{[0.25]}$ are comparable and with decreasing variance as count increases. Our modeling approach was also successfully carried out on a clinical PET/CT nonsmall cell lung dataset acquired between 2012 and 2015 on a GE Discovery VCT with OSEM reconstruction. Exploratory analyses demonstrated a similar correlation structure between our structural variables and a set of textural features comparable to the one used here.^59,60

Final feature selection based on machine learning classifiers suggests a less stable role found for radiomic features. These could either be represented significantly or on the contrary, given much lower importance depending on the method used. This also depends on the classification reference, i.e., what patient outcome or disease information is used to guide the selection problem, if no histological information is available. It seems unlikely that a consensus may be reached as to what feature selection method should be used when implementing radiomics analyses. These considerations may again change when larger sets of texture features are included or when censored time-to-event data is taken into account;^38,61 this will be the scope of future investigations.