2.2.1.

Image segmentation

Image segmentation is a fundamental process in automated image analysis, enabling the precise delineation of a tumor, its subregions, and the surrounding infiltrated anatomy. For example, segmentation of a glioblastoma, the most malignant brain tumor, can delineate the enhancing tumor (ET) and nonenhancing tumor (NET) (possibly necrotic) parts of the tumor core, its surrounding edematous (ED) tissue, and the distant normal-appearing tissue regions.^15,16 Such delineations enable the extraction of features specifically from each tumor subregion, and surrounding anatomy, thus allowing for more accurate quantification of the tumor’s entirety and spatiotemporal changes. It is important to note that the term delineation refers to regions based on their radiographic appearance and may differ from the actual tumor delineations.

The CaPTk suite offers several segmentation modules, ranging from general purpose, user-guided segmentation, to specialized segmentation methods tuned to the specific characteristics of certain tumors and organs. A representative example of general purpose segmentation is ITK-SNAP,¹⁷ a well-established interactive tool, which has been integrated into the main interface of CaPTk. ITK-SNAP is based on random forest classifiers, trained using a manual definition of tissues of interest, to produce an initial segmentation subsequently refined using level-sets. GLISTRboost,^15,16,18,19 available through the web-based image processing portal²⁰ of the Center for Biomedical Image Computing and Analytics (CBICA), leverages CaPTk as a user-friendly means of initialization. GLISTRboost performs multimodal brain glioma segmentation and atlas registration, which describes a semiautomatic hybrid generative-discriminative method. The generative part is based on an expectation–maximization framework to segment brain scans into tumor (i.e., edema, enhancing and nonenhancing tumor), as well as “healthy” tissue labels (i.e., white and gray matter, cerebrospinal fluid, vessels, and cerebellum), and incorporates the modeling of tumor growth and infiltration via reaction–diffusion–advection equations.^7,8 The discriminative part is based on a gradient boosting^21,22 multiclass classification scheme, to refine tumor labels based on information from multiple patients. Last, a Bayesian strategy²³ is employed to further refine and finalize the tumor segmentation labels, based on patient-specific intensity statistics from the multiple modalities available. GLISTRboost provides estimates of the segmentation labels, parameters of underlying tumor growth models, as well as the tumor anatomical location in a standardized anatomical system via deformable atlas registration.⁶ In addition to provide segmentation labels, such application-specialized segmentation approaches allow evaluation of cancer spatial patterns in standardized coordinate systems, which is receiving increasing attention as predictors of clinical outcome, as well as biomarkers of molecular characteristics of the underlying tumor.²⁴ Future work will result in tighter integration between the CaPTk desktop client and IPP web server.

2.2.2.

Imaging features

Analytic functions in CaPTk are based on an extensive QIP panel of features (Fig. 1), which are integrated into imaging signatures, using analytics and machine learning to produce a variety of diagnostic, prognostic, and predictive biomarkers. Although features can be common across many diagnostic and predictive tasks, the way in which they are integrated into a specific imaging signature depends entirely on the task of interest, such as predicting response to treatment¹³ or estimating underlying mutations.^10,12,25 Examples of QIP features implemented in CaPTk include: (i) multiparametric imaging signals of different coregistered protocols/modalities, such as native T1- and T2-weighted images, T1 with gadolinium (T1-Gd), T2-fluid attenuated inversion recovery (T2-FLAIR), diffusion tensor imaging (DTI), dynamic susceptibility contrast (DSC), or dynamic contrast-enhanced (DCE) MRI; (ii) textural features [e.g., co-occurrence, run-length, size-zone matrices, local binary patterns (LBPs), fractal dimensions, wavelets], capturing characteristics of the local microarchitecture of tissue. Such features have been used extensively in mammographic image analysis and breast cancer risk assessment,⁹ as well as in predictive modeling of glioblastoma;^3,4,26 (iii) histograms, reflecting various imaging signal distributions within different delineated tumor subregions. The shapes of these histograms express anatomical and functional changes caused by the tumor that result in signal changes and have demonstrated a connection to clinical endpoints, such as survival, risk factors, and underlying cancer molecular characteristics;⁴ (iv) temporal perfusion dynamics captured via principal component analysis (PCA), which have been related to recurrence and infiltration, as well as molecular tumor characteristics;^{3,10,11,25,26} (v) DTI-derived features, including fractional anisotropy (FA), radial diffusivity (RAD), axial diffusivity (AX), apparent diffusion coefficient (ADC), water-free diffusion, fiber tract connectivity and other properties; (vi) DSC-MRI derived features: peak height (PH), percent signal recovery (PSR), and relative cerebral blood volume (rCBV); (vii) spatial patterns of cancer distribution:²⁴ although previously relatively unappreciated, such spatial patterns that capture the spatial distribution and pattern of a tumor’s entirety (i.e., to which parts of the brain does the tumor extend), obtained via deformable registration to a standardized atlas space, receive increasing attention due to their relationship to prognosis and genotype.²⁷ Connectomic signatures, later discussed in Secs. 3.3 and 5.3.2, will also be incorporated in the immediate future.

An essential step in the aforementioned segmentation and feature extraction processes is image normalization and harmonization. In particular, image characteristics vary, often considerably, across different scanners, acquisition protocols, clinical centers, as well as patients. This renders the extraction of reproducible imaging features and signatures challenging. Appropriate histogram normalization is therefore particularly important. Toward this end, CaPTk provides the WhiteStripe approach²⁸ (Appendix A, Table 7), which normalizes conventional MRI by detecting a latent subdistribution of normal tissue and linearly scaling the histogram of images. In addition, histogram matching techniques available through the insight toolkit (ITK)²⁹ are also provided.

2.2.3.

Image registration

Coregistration of different imaging sequences, i.e., alignment of different anatomical regions, is very important in order to analyze in tandem voxel-by-voxel features, coming from the intensity signals of aligned imaging sequences. For example, building QIP features for analysis of brain tumors often requires the coregistration of conventional MRI scans, such as T1-weighted and T2-FLAIR but also various diffusion- and perfusion-based images. In addition to rigid registration (or fusion) of multichannel images, deformable registration is also important in cancer imaging and is used in two contexts: (i) the evaluation of temporal changes between longitudinal scans, as the tumor, its surrounding anatomy, and patient positioning may have changed between consecutive scans, and hence, appropriate registration can augment therapy response evaluation and prediction of survival and (ii) formation of population-based atlases of the spatial distribution and pattern of cancer, in order to evaluate the relationship between such distributions and molecular characteristics or clinical outcomes. CaPTk has access to a variety of image similarity metrics offered through ITK. Moreover, specialized deformable registration methods optimized for specific types of problems are available on CBICA’s IPP.²⁰ For example, the deformable registration described in Ref. 30, which is based on DRAMMS,^31,32 enables follow-up scans to be superimposed onto baseline scans in order to evaluate tumor volume changes over time as a measure of response to neoadjuvant chemotherapy for breast cancer. DRAMMS is a separate, standalone package, which is compatible with CaPTk and has been tested along with the breast CaPTk modules as a means for finding breast MRI changes.

5.2.1.

Predicting patient survival in glioblastoma patients

Glioblastoma is a disease with grim prognosis, of median survival of around 14 months after applying the standard of care, which comprises tumor resection and peritumoral radiation therapy along with chemotherapy. However, there is a fairly broad range of survival, from a few months to more than 2 years. Having baseline predictors of patient survival is important for patient management. It is also important for selecting patients with relatively homogeneous expected survival into a clinical trial, thereby likely increasing the trial’s ability to detect treatment effects, especially its ability to prolong survival. Toward this end, a preoperative multiparametric MRI of de novo glioblastoma patients was summarized in QIP signatures, capturing various characteristics of ET, NET, and ED, to estimate the likelihood of survival (Appendix A, Table 2).⁴ The initial features included (i) normalized volume of ET, NET, ED, and their combinations; (ii) distance of tumor (ET + NET) and ED to ventricles; (iii) mean and standard deviation of intensities of T1, T2, T1-Gd, T2-FLAIR, rCBV, PH, PSR, FA, RAD, AX, and ADC in ET, NET, and ED; (iv) frequency of intensities of T1, T2, T1-Gd, T2-FLAIR, rCBV, PH, PSR, FA, RAD, AX, and ADC in each distribution bin of ET, NET, and ED; (v) location of the tumor in the brain; and (vi) age. All features were integrated via a support vector machine (SVM) configuration to build two predictive classification models: a 6- and an 18-month SVM model to differentiate between patients surviving less/more than 6 months (short-survivors) and 18 months (long-survivors), respectively. The group of subjects having survival between 6 and 18 months was considered mid-survivor group. Forward feature selection was applied on the training set only to select important features. The SVM scores from each model were combined to calculate a composite survival prediction index (SPI), where higher SPI refers to relatively longer survival. We developed our predictive models on a discovery cohort ($n=110$) and tested these prospectively on a replication cohort ($n=57$) of glioblastoma patients. The two-class balanced accuracy of the 6-month and 18-month models was 75.06% ($\mathrm{AUC}=0.79$) and 77.85% ($\mathrm{AUC}=0.77$) in the replication cohort, respectively. Overall three-class classification accuracy into long, medium, and short survival groups was $\sim 70.18\%$ in the replication cohort. Kaplan–Meier survival curves [$p\text{-}\text{value}<0.001$, log-rank (Mantel–Cox)]⁵³ and hazard ratios were also computed for survival analysis at a 95% CI (Fig. 5).

Fig. 5

Kaplan–Meier survival curves for the replication cohort. Actual survival on $x$-axis is compared among each of the three survival groups based on predictions generated by the SPI. med, medium SPI HR, Hazard ratio.

5.2.2.

Predicting patient survival in breast neoadjuvant chemotherapy patients

The DCE-MRI images were analyzed for a subset of 106 women with complete imaging data available, which were recruited as part of the ACRIN 6657/I-SPY-1 trial.⁵⁴ A baseline model was created with age, race, hormone receptor status (ER/PR/Her2), and functional tumor volume (FTV) after neoadjuvant chemotherapy. Utilizing the DRAMMS deformable registration method^30,20 adapted to breast MRI, features regarding spatial/temporal changes between longitudinal scans (before and during the first patient visit after the initiation of neoadjuvant chemotherapy in this study), including voxel-wise volume ratio (Jacobian), as well as parametric response maps (PRM) for kinetic features, were estimated. Kinetic features included signal enhancement ratio, peak enhancement (PE) and wash-in/wash-out slope (WIS/WOS). To quantify heterogeneity for each feature, discrete wavelet transformation was used and then PCA was applied to reduce the number of estimated wavelet coefficient into the top two principal components, expressing 60% to 80% of the total variance (Appendix A, Table 6). The Cox proportional hazards model was utilized to perform a time-to-event analysis and predict recurrence-free survival (33 events) by estimating the c-statistic.⁵⁵ The baseline model using the standard clinical covariates and FTV was compared with the model, where registration-derived and PRM features were added. The c-statistic was 0.70 ($p\text{-}\text{value}<0.001$) for the baseline model, whereas the augmented model with the PRM features and the Jacobian information improved the c-statistics by 0.73 ($p\text{-}\text{value}<0.005$) and 0.74 ($p\text{-}\text{value}<0.001$), respectively. A model including both Jacobian and PRM features had the highest c-statistic of 0.77 ($p\text{-}\text{value}<0.001$; Fig. 6).

Fig. 6

Survival curves as function of (a) FTV, after first visit during neoadjuvant chemotherapy and (b) Jacobian heterogeneity when FTV is greater than the mean value (between images before and first visits).

5.2.3.

Predicting treatment response and survival of early-stage nonsmall cell lung cancer

To identify radiomic biomarkers for predicting treatment response and survival of early-stage nonsmall cell lung cancer (NSCLC) in patients, who received stereotactic body radiation therapy (SBRT), we carried out a radiomic analysis using CaPTk to distinguish patients with different treatment response and investigated the association between subclusters of tumor phenotypes and clinical outcomes. This study was performed based on a longitudinal fludeoxyglucose-positron emission tomography (FDG-PET)/computed tomography (CT) dataset of 80 patients, who were treated with SBRT for stage 1 NSCLC with over 2-years median follow-up. All patients in this dataset had a solid component of their NSCLC tumor, and some also had an additional ground glass component. Although all these patients were treated uniformly ($12.5\mathrm{Gy}\times 4\text{fractions}$, or $10\mathrm{Gy}\times 5\text{fractions}$), they had different primary tumor outcomes. From each patient’s standardized uptake values of FDG-PET scan collected before the treatment, we extracted 343 radiomic features, including intensity statistics, gray level co-occurrence matrix, and gray-level run-length matrices, as well as LBP within the tumor region. Then, the patients were grouped into two clusters with distinctive radiomic features using an unsupervised clustering analysis method³⁹ (Appendix A, Table 4). Kaplan–Meier survival analysis with respect to death and nodal failure at group level was performed for each cluster of patients (Fig. 7). Significant differences were observed for survival ($p\text{-}\text{value}=0.0004$, Log-rank test) and nodal failure ($p\text{-}\text{value}=0.001$).

Fig. 7

Survival analysis of two clusters of the early-stage NSCLC patients with respect to (a) death and (b) nodal failure.

5.3.1.

Predictive modeling of peritumoral infiltration and recurrence

Current practice in treating glioblastoma includes resection guided by imaging-based tumor margins (typically defined via enhancement in T1-Gd images), followed by uniform radiation of peritumoral brain tissue. It is well known, however, that glioblastomas infiltrate their surrounding brain tissue, especially in the peritumoral edematous tissue defined by high T2-FLAIR signal (i.e., ED). Ability to predict the regions in this peritumoral tissue that are heavily infiltrated and most likely to present early recurrence can considerably change clinical practice, by guiding aggressive supratotal resection, i.e., targeted resection of peritumoral tissue, as well as targeted elevated radiation dose in nonresected peritumoral tissue that is more likely to present early recurrence. Based on functions provided by CaPTk, it has been recently shown that the combination of multiparametric MRI and machine learning can lead to predictive models of tumor infiltration that highlight peritumoral tissue that is $\sim 10$ times more likely to recur³ (Appendix A, Table 1). In particular, two regions were selected within the edema region to train the model. The near region was defined as the area immediately adjacent to the tumor, whereas the distal edge of edema was designated as the far region. The signal intensity of T1, T1-Gd, T2, T2-FLAIR, AX, FA, RAD, ADC, and rCBV, and first five principal components derived from the DSC-MRI image were combined via Gaussian kernel function of SVM. The model was retrospectively cross-validated on a cohort of 31 patients [$\text{odds ratio}=11.17$ (99% CI: 10.71 to 11.64)] and was subsequently evaluated on a prospective cohort of 34 patients [$\text{odds ratio}=9.29$ (99% CI: 8.95 to 9.65)]. The rest of this section presents a case report from a clinical case that we recently processed and presented at the weekly Brain Tumor Conference of the Hospital of the University of Pennsylvania, using these methods and hence highlighting their potential for routine clinical use.

Case report: SN was 25 when she presented with headaches and dizziness. An outside brain MRI revealed a nonenhancing, expansile, right posteromedial occipital lesion, extending into gray matter. Low-grade glial neoplasm was top on the differential diagnosis. Three months later, a subsequent MRI with brain tumor protocol, including perfusion (DSC) and spectroscopy, was consistent with that diagnosis. Another 3 months later, a follow-up brain MRI demonstrated an increase in size of the lesion, with new intralesional enhancement, triggering a change in radiologic diagnosis to high-grade glial neoplasm. A repeat MRI confirmed that diagnosis, demonstrating increased rCBV within the enhancing component. The thought was that there had been malignant transformation of the original tumor, which was identified 6 months earlier. The lesion was resected 10 days later, and the postoperative MRI showed some residual abnormal T2 signal but no residual enhancement. The pathology diagnosis was glioblastoma. A follow-up MRI, 2.5 months after initial resection demonstrated new enhancement, is consistent with tumor progression/recurrence. Using multiparametric analysis with CaPTk, we determined the areas of high likelihood for recurrence on the preoperative scan and superimposed the corresponding probability maps [Fig. 8(a)]. The actual recurrence [Fig. 8(b)] was centered precisely on the predicted high likelihood area.

Fig. 8

(a) Sagittal postgadolinum T1-weighted images with recurrence probability maps on preop scan, calculated via CaPTk. (b) Actual recurrence scan, about 3 months later.

5.3.2.

Fiber tracking

Two of the major challenges faced by fiber tracking in the realm of neurosurgical planning are that the reconstruction of tracts is affected by mass effect, when the tracts are displaced and distorted, and edema and infiltration, when the tracts are broken as a result of the change of diffusion parameters due to the pathology. These underline the need for methods that can track through edema and reconstruct even partial and displaced tracts. We are developing tractography algorithms that are invariant to edema, based on multicompartment modeling of the diffusion data.^41,42 The improvement in tracking with this modeling, over-traditional tracking, can be seen in Fig. 2(a). We have developed automated atlas-based tract reconstruction called Confetti (connectivity-based fiber extraction and identification;⁴³ Appendix A, Table 8). The applicability, the reliability, and the repeatability of the Confetti were validated in a dataset of healthy individuals acquired repeatedly.⁴³ Compared to the clustering of fibers for each scan independently, our framework provided better test–retest reproducibility results, with decreased (25%) mean intraindividual distance (i.e., disagreement of clusters between different timepoints of the same individual), while preserving interindividual differences. Additionally, the Confetti was also tested in tumor patients⁴⁰ on six major fiber bundles: cingulum bundle; fornix, uncinate fasciculus (UF), arcuate fasciculus, inferior fronto-occipital fasciculus, and inferior longitudinal fasciculus (ILF). The agreement between clustering and experts as quantified by Cohen’s kappa ranged between 0.6 and 0.76. Except two tracts, ILF and UF, the agreement between clustering and experts was higher than agreement between experts themselves, highlighting the reliability of the paradigm. When the tumor demonstrated significant mass effect or shift, the automated approach was useful to provide an initialization to guide the expert with identification of the specific tract of interest.

We have developed a measure of injury called disruption index of the structural connectome (DISC), based on how the information transfer in the brain is affected by changes in different regions.⁵⁶ A representation of the vulnerability of the brain to injury can be seen in Fig. 2(b). We have tested the DISC on a dataset of traumatic brain injury (TBI) patients with moderate to severe TBI examined at 3 months postinjury.⁵⁶ DISC was significantly correlated with post-traumatic amnesia (Pearson $r=0.52$, $p$-value: 0.0007), verbal learning (Pearson $r=-0.42$, $p$-value: 0.0075), executive function (Pearson $r=-0.41$, $p$-value: 0.0083), and processing speed (Pearson $r=-0.58$, $p$-value: 0.0001), demonstrating that assessing structural connectivity alterations may be useful in development of patient-oriented diagnostic and prognostic tools.

These diffusion MRI-based tools are being developed as a separate suite in CaPTk because of special visualization needs for such data. A preliminary view of the proposed surgical plan can be seen in Fig. 2(c) in which the automatically extracted tracts are presented with regard to the tumor. The tracts, vulnerability maps, and recurrence maps will be incorporated in this plan in the future. Radiation plans can be created by using the resection cavity instead of the tumor.

5.4.1.

Imaging signatures of the EGFRvIII mutation, as well as transcriptomic subtype, in glioblastoma

CaPTk offers an imaging signature highly distinctive of the EGFRvIII mutation in glioblastoma. This imaging signature leverages the heterogeneity of DSC-MRI signals throughout the peritumoral region, which is depicted by abnormal/bright T2-FLAIR signal¹⁰ (Appendix A, Table 3). In particular, this signature is based on the observation that the gradient of perfusion between tissue immediate to the active tumor and tissue distant from the tumor but within this peritumoral edematous region (bright FLAIR) is significantly higher in tumors that do not harbor the mutation (i.e., EGFRvIII- glioblastoma). This observation is reflecting the more local invasion and neovascularization of the EGFRvIII- tumors, and vice versa, the potentially deeper invasion of EGFRvIII+ tumors. Principal component features are extracted from the DSC-MRI images using CaPTk and the Bhattacharya distance is calculated to form the peritumoral heterogeneity index (PHI). The method was evaluated in preoperative perfusion scans of independent discovery ($n=64$) and validation ($n=78$) cohorts. Analysis in cohorts demonstrated high accuracy (89.92%), specificity (92.35%), and sensitivity (83.77%), with significantly distinctive ability ($p\text{-}\text{value}=4.0033\times {10}^{-10}$, $\mathrm{AUC}=0.8869$). Figures 9(a)–9(c) show the clear separation between PHI values of EGFRvIII(+) and EGFRvIII(−) tumors.

Related work in glioblastoma was applied by obtaining a more extensive QIP integrating T1, T1-Gd, T2, T2-FLAIR, DTI extracted measurements, principal components of DSC-MRI, along with parameters of the tumor growth model, spatial location, and other features. The multiparametric obtained signature aimed to detect molecular subtypes of glioblastoma, as described in Ref. 57. Figure 9(d) shows ROC curves obtained in this four-way classification experiment (baseline “chance” level is $\sim 25\%$). The AUC for classical, mesenchymal, neural, and proneural subtypes, respectively, was 0.75, 0.89, 0.92, and 0.87.

5.4.2.

Breast DCE-MRI phenotypes correlate to gene-expression based breast tumor profiling

Preoperative breast DCE-MRI images from women with ER+ breast cancer were retrospectively analyzed.^25,58 All women had their primary tumor tested with Oncotype DX,⁵⁹ an assay that measures RNA expression of 21 genes from formalin-fixed paraffin-embedded tissue and provides a score recurrence risk 10 years after treatment ($\text{low}\le 17\%$, $18\%<\text{medium}<30\%$, $\text{high}\ge 31\%$). Validated morphologic, kinetic, and spatial heterogeneity features were extracted from each primary tumor: tumor area and perimeter were used to measure tumor size and ellipticity and convexity were computed to capture shape and structure.^25,58 Voxel-wise kinetic tumor features of PE, time-to-peak, WIS/WOS were also estimated, from which statistics were calculated to capture spatial kinetic tumor heterogeneity.²⁵ Multivariate linear regression was used to test associations between radiomic features and the gene-expression based recurrence score. To identify intrinsic imaging phenotypes, unsupervised clustering was applied on the extracted feature vectors.²⁵ There was significant correlation ($r=0.78$, $p\text{-}\text{value}<0.01$) between MRI features and the recurrence score.⁶⁰ Four imaging phenotypes were detected (Fig. 10), with two including only low and medium recurrence risk tumors. Tumors with a gene expression profile at high risk of recurrence showed a predominantly rapid contrast uptake, suggesting high levels of perfusion and vessel permeability. When phenotypes were used in a model to predict recurrence risk, the AUC reached up to 0.82 ($p\text{-}\text{value}<0.01$).