Quantitative imaging of cancer in the postgenomic era: Radio(geno)mics, deep learning, and habitats.

Although cancer often is referred to as "a disease of the genes," it is indisputable that the (epi)genetic properties of individual cancer cells are highly variable, even within the same tumor. Hence, preexisting resistant clones will emerge and proliferate after therapeutic selection that targets sensitive clones. Herein, the authors propose that quantitative image analytics, known as "radiomics," can be used to quantify and characterize this heterogeneity. Virtually every patient with cancer is imaged radiologically. Radiomics is predicated on the beliefs that these images reflect underlying pathophysiologies, and that they can be converted into mineable data for improved diagnosis, prognosis, prediction, and therapy monitoring. In the last decade, the radiomics of cancer has grown from a few laboratories to a worldwide enterprise. During this growth, radiomics has established a convention, wherein a large set of annotated image features (1-2000 features) are extracted from segmented regions of interest and used to build classifier models to separate individual patients into their appropriate class (eg, indolent vs aggressive disease). An extension of this conventional radiomics is the application of "deep learning," wherein convolutional neural networks can be used to detect the most informative regions and features without human intervention. A further extension of radiomics involves automatically segmenting informative subregions ("habitats") within tumors, which can be linked to underlying tumor pathophysiology. The goal of the radiomics enterprise is to provide informed decision support for the practice of precision oncology.

Introduction

Although cancer is commonly described as a “disease of the genes,” it also is indisputable that the genetic and epigenetic properties of individual cancer cells are highly variable, even in different regions of the same tumor, or between metastatic sites within the same patient. By generating heterogeneous subpopulations, this genetic diversity can lead to the failure of targeted therapies or chemotherapies, even with validated targets and drugs, as resistant clones commonly emerge and proliferate. Even with the most promising of immune therapies, only approximately 10% to 20% of patients exhibit durable responses. Hence, approaches to characterize and quantify the extent of intratumoral heterogeneity in individual patients might be useful for guiding therapies that adapt during the course of treatment. Herein, we propose that a form of image‐based quantitative analysis known as “radiomics” can be used to quantify heterogeneity, predict outcome, and longitudinally monitor responses.

The practice of radiology always has been focused on the interpretation of images, which, until recently, was facilitated entirely by human observation and recollection. Although there are exceptions (eg, in nuclear medicine, in which localized metabolic activity can be quantified as a specific uptake value), the characterizations of object location, shape, sharpness, and intensity are accomplished by highly trained human observers. These interpretations are critical for disease and patient management, including diagnosis, prognosis, staging, and assessment of treatment response. However, to our knowledge, with few exceptions, this is largely a subjective operation with limited specificity and high interreader variability.1

In 2007, Segal et al showed that careful characterization of the appearance of liver lesions in contrast‐enhanced computed tomography (CT) scans could be used to reach a conclusion regarding their molecular properties.2 Similarly, in 2008, Brown et al demonstrated that quantitative analysis of gray scale image texture on magnetic resonance imaging (MRI) scans of patients with oligodendroglioma could be used to predict their genetic signatures, specifically codeletion of chromosomes 1p and 19q.3 These early examples indicated that not only could the appearance of noninvasive imaging be quantified, but that these data could illuminate the fundamental molecular properties of cancers. Early examples used “semantic features” (ie, categorizations of human observations using a controlled vocabulary),4, 5, 6 and these have expanded to the use of “computational features,” which are direct mathematical summarizations of image regions. Together, semantic and computational feature classes are the basis of the new field of “radiomics,”7, 8 defined as the “high‐throughput extraction of quantitative features that results in the conversion of images into mineable data,”9, 10 and feature prominently in what today is called quantitative imaging.11, 12, 13 The following section describes the radiomics workflow, including its strengths and its challenges, and highlights the integration of radiomics with clinical data and the molecular characterization of tumors, also known as radiogenomics,14, 15, 16, 17, 18, 19, 20, 21 for the building of predictive models. This includes both “conventional radiomics” (ie, machine computation of human‐engineered image features) and artificial intelligence, or “deep learning,” (DL) to automatically discover the most informative features for image phenotype description from the image and clinical data. In the second section, we discuss an emerging field, called “habitat” imaging, which combines radiomics with multiparametric images to automatically classify the pathophysiological nature of subregions within tumors and surrounding tissues. All 3 approaches–conventional radiomics, DL, and habitat imaging–have shown high promise in predicting, making prognoses, and monitoring responses of patients with cancer to therapies.

Conventional Radiomics

Figure 1 illustrates the radiomics workflows for a volumetric data set. Component parts are: 1) the identification of the location of the subject tumor; 2) annotation of the tumor with semantic features; 3) tumor segmentation (ie, identification of the entire volume to be analyzed); and 4) image feature computation via human‐engineered image features (“conventional radiomics”) and/or via DL approaches (“DL radiomics”). In some cases, “delta” features22, 23, 24, 25, 26 may be computed by comparing individual feature values derived from different imaging phases or images obtained at different times during the course of therapy, and/or feature vectors may be created that combine computed features from multiple imaging methods.27, 28, 29, 30

Figure 1

Radiomics workflows. The conventional workflow computes human‐engineered features, whereas the deep learning workflow explores features during training using labeling from clinical and/or demographic data. Either workflow can be used, or both can be used and the resulting features combined. The deep learning workflow can use, but may not require, tumor identification or segmentation (hence “alternate pathway”). Both workflows also can use semantic features for the final classification stage, which generates an output such as benign/malignant, responder/nonresponder, probability of 5‐year survival, etc.

Tumor identification

Each and every tumor to be processed first must be identified, either semiautomatically or manually by a radiologist or automatically using computer‐aided detection approaches. When multiple tumors are present in a single imaging study, human effort generally is required to identify those that are clinically relevant, as in the case of index lesions scored with RECIST (Response Evaluation Criteria in Solid Tumors). When delta radiomics features are required, matching of tumors across observations also will be required.

Annotation with semantic features

Semantic features are descriptive observations of image content. For example, semantic features of a lung tumor might include “left lower lobe,” “pleural attachment,” “spiculated,” “ground‐glass opacity,” etc. However, the extraction of semantic content from unstructured radiology reports is not appropriate due to inconsistent and/or ambiguous vocabulary across observations and observers. Furthermore, structured reports may not support the types of detailed observations required for making fine distinctions among tumor characteristics that may prove useful in classification or assessing response. Advances in structured reporting include the description of semantic features using a controlled vocabulary, such as RAD‐LEX.31 However, although much effort has been expended to develop more structured reporting, it is still not common in the broad radiology community. Currently, there are at least 2 systems that facilitate semantic annotation of radiological images using controlled vocabularies: electronic Physician Annotation Device (ePAD) (Fig. 2)32, 33 and the Annotated Image Markup (AIM) data service plug‐in34 to ClearCanvas.35 Some semantic features, such as location, other morbidities, etc, are meant to be complementary to computational features, whereas others, such as “spherical,” “heterogeneous,” etc, are correlated with computational features.36, 37 One advantage of semantic annotations is that they are immediately translatable (ie, they can be ascribed in clinical environments without specialized algorithms [eg, segmentation] or workstations). As such, they have shown interesting results in several radiogenomic studies.20, 38, 39, 40, 41, 42, 43

Figure 2

Example of semantic annotation of a lung tumor that is part solid and part ground glass using ePAD. After tumor segmentation (green), either manually within ePAD or created via other means, an observer selects annotations using a custom template (eg, one that was built using the AIM Template Builder33) or one of several available templates. The example in this figure shows a subset of the annotation topics that are required to complete this annotation.

Tumor segmentation

Radiomics features can be extracted from arbitrary regions within the image volume. A given region may contain an entire tumor, a subset of the tumor (eg, a habitat), and/or a peritumoral region thought to be involved with or affected by the tumor. In all cases, these regions must be unambiguously identified (segmented) and input to the radiomics feature computation algorithms. This segmentation step is the single most problematic aspect of conventional radiomics workflows, because the features computed from tumor volumes may be extremely sensitive to the specification of the volume to be analyzed. Each combination of tumor type and image modality presents its own challenges (including volume averaging of tissues within each voxel, tumor contrast with surrounding/adjacent structures, image contrast‐to‐noise characteristics, and variations of image quality across vendor implementations and time). In addition, many algorithms require operator inputs, such as bounding boxes and/or seed points, and the segmentation outlines and radiomics features computed from them may be sensitive to these inputs. For example, a recent study demonstrated wide variations in lung nodule segmentation outlines across 3 different algorithms, and even across individual algorithms initialized with different user inputs.44 Another study using the same data set revealed that many radiomics features, particularly those that quantify shape and margin sharpness, are quite sensitive to segmentation.45 Thus, the state of the art today is such that each segmentation must be reviewed and possibly edited by a human observer in order for the radiomics features computed from it to be trusted. One potential mitigation is to ignore segmentation altogether and to compute only features that do not require complete edge‐to‐edge coverage of the tumor (ie, histogram and texture features), which may be less sensitive to the exact tumor definition, and to ignore shape and margin sharpness features, which require accurate and consistent edge delineations.46, 47

Image feature computation

Conventional or human‐engineered computational image features can be divided into 4 classes, namely those that describe shape48, 49, 50; margin sharpness51, 52; histogram features (eg, mean, variance, kurtosis, maximum, and minimum); and texture features,53, 54, 55, 56, 57, 58, 59, 60 which describe the spatial variation of gray values within the tumor. Within each class, there are hundreds to thousands of individual features: for example, some texture features quantify the spatial variation of gray values across multiple scales and orientations, and shape features can similarly quantify edge irregularity at multiple scales. It is important to recognize that many image features are intercorrelated and, as a result, not all features may add independent predictive power to radiomics models. A standard approach in many studies is to generate an autocorrelation matrix and combine correlated features into a single descriptor.61 Several groups have made available computer code and processing pipelines for the calculation of image features from volumetric image data and segmentations (or at least volumes of interest), for example the Imaging Biomarker Explorer,62 the Quantitative Image Feature Engine,63 and pyRadiomics.64

Challenges to conventional radiomics

The ultimate usefulness of conventional radiomics workflows is to generate image features that then can be integrated with other medical data such as diagnosis, survival, response to therapy, mutations or genomic profiles, and demographics for the purposes of building predictive models for ≥1 of the clinical variables. Space does not allow room for details, but we note that a primary concern in these analyses is overfitting (ie, building statistically significant models based on hundreds to thousands of image features requires data from many multiples of that number of individuals). As noted, many of the features are correlated and therefore feature reduction techniques based on redundancy, relevance, and/or sparse regression should be used wherever possible.65, 66 Although many investigators use cross‐validation methods to avoid training and testing on the same data, it still is important to test each predictive model on a completely independent cohort to assess generalizability.

An additional challenge is the need to standardize the methods for the calculation of radiomics features so that identically intended features computed from the same data by different algorithms have the same name and values. Indeed, the same study referenced above that compared features computed from multiple segmentations also indicated that implementations from separate institutions of purportedly the same feature sometimes produced different values.45 Much effort currently is underway to standardize feature naming and computation conventions, predominately led by the Image Biomarker Standardization Initiative67 and the Quantitative Imaging Network.11, 12, 13

In addition to the challenges raised by segmentation, a final challenge is the sensitivity of radiomics features to image acquisition and reconstruction (ie, the heterogeneity of image acquisitions). Each clinical study uses their own combination of acquisition parameters, such as slice thickness, reconstruction kernel, MR pulse sequences, etc. In addition, many acquisition parameters are optimized for the particular patient being studied (eg, kilovoltage, mass, and field of view). Although radiologist interpretations are somewhat immune to these differences, computational radiomics features are by design sensitive to these choices. In addition, random image noise also may affect many radiomic computations.61, 68 Thus, a predictive model built for lung nodule characterization may be validated in one cohort and fail in another because of differences in the data acquisition and reconstruction. Although there have been some efforts to compensate for these differences (eg, using phantoms and image interpolation),68, 69, 70, 71, 72, 73 and herculean efforts on protocol standardization by the Quantitative Imaging Biomarkers Initiative (QIBA)74, 75, 76, 77, 78 currently are underway, much work remains before radiomics analyses can derive usefulness from large databases collected by multiple institutions. Another possibility to mitigate these issues is through a branch of artificial intelligence known as DL.

Deep Learning

In conventional radiomics workflows, accurate lesion segmentation and user‐defined features play key roles. As described above, there are challenges to semiautomated segmentation, such as interoperator variability and time consumption.79, 80, 81, 82, 83, 84, 85, 86 Although there are many hundreds of descriptive features designed by prior knowledge, the current feature sets still may not be optimal for a given task.87

DL can mitigate these limitations because it does not require accurate segmentation and it creates its own features through multiple layers of learning.87, 88 Only recently applied to radiomics, DL has proven to be valuable in both differential diagnosis89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103 and prognosis.104, 105, 106, 107 DL involves abstraction by building networks with >2 processing layers.108 For medical images, convolutional neural networks (CNNs) are the most common models because they can accept 2‐dimensional (2D) or 3‐dimensional (3D) images as input. The first CNN was proposed by LeCun et al in 1998,109 but its success was limited until the advent of graphic processing units and the development of learning algorithms.93, 110, 111, 112

Figure 3 shows examples of CNNs, which can be used to classify data types as diverse as purchasing preferences and satellite images. With medical imaging, the input layer during training includes images or subregions of labeled images, which then are convolved in sublayers along with their known classifiers (for example, benign or malignant) to identify those image features that are most related to the classification. Because images contain many common features that are relevant to classification (eg, edges, shapes, colors), the core layers of CNNs relating to these can be transferred to other classification tasks. Indeed, although only recently applied in radiomics, DL has proven valuable in both differential diagnosis89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103 and prognosis prediction.104, 105, 106, 107

Figure 3

Schematic representation of convolutional neural network (CNN) architecture and the concept of “transfer learning.” Because images from multiple sources have common salient features (borders, shapes, etc), the core of a CNN trained for 1 task (Top row: eg, cat from pot‐au‐feu from camping tent) can be “transferred” (ie, used without modification) for a second task (Bottom row: benign from malignant tumors in computed tomography images). Because only the distal layers (fully connected [FC] I and FCII) remain to be trained, much less training data are required for the second task. Max indicates maximum.

Two major factors influencing CNN applications in diagnostic imaging are computational power and the availability of training data. Although the solution to computing power is simply time and/or money, the solution to the limited availability of well‐annotated training data sets is not so simple. One solution can be to limit the size of the input data. For example, small 2D patches can be used to significantly reduce computational costs and required training cases. Even using only 2D images, CNNs have been shown to improve accuracy compared with conventional radiomic approaches. Examples include tumor grading with MRI,100 classification of pulmonary nodules with 2D CT,113 and predicting response to neoadjuvant therapy with positron emission tomography (PET) images.114

To better use volumetric information from the 3D medical images, information from 2D patches of the 3 orthogonal views (axial, sagittal, and coronal) can be used, which is a method referred to as 2.5D.115 These 3 views can be inputted into existing CNNs trained on 3 channels (often the red‐green‐blue channels of natural images).116 If enough computational power and training cases are available, full 3D CNNs can be used to provide better results.105, 117, 118 In the past 2017 Data Science Bowl focusing on lung cancer detection,119 full 3D CNNs obtained the best performance. In addition to the different input and dimensions of the networks, many different architectures also have been attempted by incorporating unique task‐specific properties such as multiview120 and multiscale to achieve higher accuracy.98, 99, 101, 102, 121

An approach for solving the issue of limited data is to use “transfer learning.”122, 123 This approach uses an unrelated training set to preconfigure the CNN to reduce the size of the task‐specific training set. There currently are 2 popular transfer learning strategies, one of which uses an “off‐the‐shelf” set of CNN features and the second of which performs unsupervised pretraining on unrelated images/classification task (eg, dogs vs cats in images from the Internet) and fine‐tuning on medical target images. In either strategy, a pretrained network from a larger image data set could be used as a feature extractor. The first n layers of a pretrained CNN model are used as the first n layers of a target network, and the remaining layers of the target network then are retrained toward the target task (Fig. 3). The errors from the new task can be back‐propagated into the pretrained features to fine‐tune them to the new task.124 Thus, the outputs extracted from the initial layers are regarded as meaningful features that can be used to train a separate classifier, which can include image data as well as clinical, demographic, and genomic data.101 This has been used, for example, using the gastrointestinal stromal tumors standard graphics library125 or bag of visual words126 descriptors to train a CNN for chest pathology,91, 92 followed by combining generic CNN features with handcrafted features to significantly improve results in the detection of lung nodules.127 Using a similar approach, a multiview CNN model based on a publicly available pretrained CNN‐F model128 was used for estimating a patient’s risk of developing breast cancer.94

There also are solutions to create computerized “pseudodata” to improve the performance of DL models. This involves generating additional examples by perturbing an image by multiple transformations while keeping the class unchanged. A generic and accepted current practice is to perform geometric augmentation including image reflection, cropping, translation, and rotation. Noting that the final aim of these studies is to achieve a classification (eg, “normal” vs “cancer,” “benign” vs “malignant”), when 1 class is overrepresented, only a small percentage of these cases are informative. Because the pathological classes generally are underrepresented in the cancer imaging population, using multiple transforms of these images for training the CNN can improve efficiency and accuracy.103

A common problem with CNNs is “overfitting,” which occurs when the training set is too small compared with the number of layers or a trained model fails to generalize. This can be mitigated with regularization techniques, which add penalties for complexity or extreme parameter values.129, 130 More recently, CNN performances have been enhanced by propagating errors between layers to focus on the most robust and relevant features.131 Several other network architectures, such as Boltzmann machine,132 deep belief network,133 or stacked autoencoders,134 also have been used successfully in the detection, diagnosis, and prognosis of cancer. Although these networks can be applied to solve multiple classification problems, it is critical that the network is trained for the specific task at hand, and to test the classifier in a completely independent data set, preferably from a different institution.

Finally, and relevant to the previous and next sections, it is well known that multimodality images often provide additional information compared with those from a single modality; consequently, recent work using deep CNNs with multimodality images have attracted much attention.132, 135, 136, 137 For example, for the task of pulmonary nodule detection, Teramoto et al136 combined CT and PET images of pulmonary nodules into a CNN and the outputs were used as inputs to 2 support vector machines, resulting in 50% fewer false‐positive detections compared with a previous study.138 In another multimodal study, a 2D CNN was trained for low‐grade glioma tumor segmentation using T2, fluid‐attenuated inversion recovery (FLAIR), and T1 contrast‐enhanced images, and image features were obtained from its final convolutional layers. These “learned” features demonstrated improved accuracy compared with conventional radiomics for predicting IDH1 mutation status in patients with brain cancers.137

Habitat Imaging

One of the most well‐studied hallmarks of cancer is dysregulated angiogenesis, in which the tumor cells locally promote new blood vessels to provide nutrients and oxygen to maintain their fitness. However, because cancer cells are independently promoting local angiogenesis without regard for the rest of the tumor, a chaotic and leaky vasculature results.139, 140, 141 This leads to stochastic and cyclical changes in blood flow142 that, in turn, have profound effects on the local environmental properties, and the resulting evolutionary selection pressure applied to cancer cells.143 Spatial variation in environmental stressors drives a natural selection for cells that are most fit within their specific patterns of oxygen and nutrient availability, and the buildup of acidic waste products.143, 144, 145, 146 Thus, these subregions, or habitats, are expected to contain cells with similar genotypes and phenotypes, including those that confer therapy resistance.147, 148, 149, 150 In fact, intratumoral spatial heterogeneity has been associated with resistance to therapy,151 tumor growth,152 and poor prognosis,153 and therefore has direct clinical consequences. Herein, we describe approaches to characterize the presence of these habitats using multiparametric MRI (mpMRI) and emergent analyses of PET.

In this context, mpMRI can go beyond its common use to characterize tumor size and anatomy. Although response to therapy is conventionally measured by tumor size (RECIST),154, 155 functional MRI allows the assessment of tissue metabolism and physiology.156, 157 Thus, different MRI sequences can be used to measure the spatial differences in cell density, tissue organization, perfusion, and metabolism, which potentially provide earlier indications of therapeutic efficacy than tumor shrinkage as observed in size measurements.158 The combined application of images produced by multiple MR pulse sequences (mpMRI) can simultaneously assess multiple functional processes at the cellular and molecular levels.159 Anatomic (T2‐weighted) MRI commonly is used in combination with diffusion‐weighted MRI (DWI‐MRI) and dynamic contrast‐enhanced MRI (DCE‐MRI), which are sensitive to cellularity and vascular perfusion, respectively.160 Normally, these are evaluated independently and/or are combined into an average parameter for the entire region of interest, discarding important spatial information.144, 161 With habitat imaging, these orthogonal data sets are combined into a single parametric map. However, multiparametric analysis on a per‐pixel basis to generate these maps remains challenging, because it requires spatial registration of all imaging sequences and must take into consideration resolution differences and deformations that can occur between scans.160, 162, 163

Based on the assumption that MRI scans are analogous to satellite images, the image processing technology from the National Aeronautics and Space Administration (NASA) first was adapted to segment MRI data into tissue populations as separate “scene elements.” To the best of our knowledge, Vannier et al were the first to apply this methodology in a series of articles164, 165, 166 wherein mpMRI data served as input to pattern recognition procedures implemented as feature extraction followed by classification. This enabled the production of “theme maps” through supervised and unsupervised pattern recognition procedures. By plotting the gray scale frequency distribution of each image, the data tended to form clusters, which were examined against the original MRI scans and classified with tissue characteristics. The supervised classification technique started by defining a region of interest over the known structure, which served as a signature with which to train for the recognition of other structures that presented with the same classification. The unsupervised method analyzed mpMRI to classify tissue components without any prior information regarding the identification of the tissues. These analyses were able to differentiate brain parenchyma from hematoma, skin from cortical bone, fat from air, and cerebrospinal fluid from its surroundings. More recently, studies have evaluated intratumoral heterogeneity using specific MR images (eg, parameters derived from DWI‐MRI or DCE‐MRI).

DWI‐MRI measures the random movements of water molecules, which are sensitive to cellular density in tissues.159 It can be used alone or can be analyzed by calculating the apparent diffusion coefficient (ADC) using multiple diffusion weightings, which can be used as a noninvasive biomarker for diagnosis, tumor grading, and prognosis.167, 168 In addition, because changes in ADC can be observed earlier than changes in the tumor size, DWI‐MRI is a useful technique for the assessment of therapy response.159, 169 Low ADC values indicate restricted diffusion, which commonly is observed in solid tumors.170, 171 Conversely, areas of necrosis and edema are expected to demonstrate high ADC values.159, 160 In preclinical studies, these data can be spatially coregistered with histology (eg, demonstrating an increase in ADC in the one‐half of a prostate xenograft that received radiotherapy, which was associated with a decreased nuclear count and increase in extracellular space and nuclear size).169 DWI‐MRI scans of MDA‐MB‐231 and MCF‐7 xenografts were used to provide quantitative texture information regarding intratumoral heterogeneity by correlating parameters derived from ADC maps with histological features.172 In this study, histogram‐based measurements of ADC values, such as standard deviation, skewness, and kurtosis, were found to be positively correlated with Ki‐67 staining. Histologically defined microvessel density was found to be inversely correlated with ADC values. Thus, the less heterogeneous texture parameters derived from DWI‐MRI and observed in the MDA‐MB‐231 xenografts corresponded to higher vascularity (CD34) and diffuse multifocal necrosis identified in histology. The incorporation of texture features derived from the combination of ADC and T2‐weighted MRI increased the accuracy of the classification of Gleason patterns in patients with prostate cancer when compared with classifiers using only the ADC mean instead of the texture features.173

DCE‐MRI reflects vascular flow and permeability by tracking the passage of a bolus injection of intravenous gadolinium‐based contrast through the tissue of interest using T1‐weighted images before and after the contrast administration.174 DCE‐MRI can be analyzed quantitatively by fitting the data to ad hoc or pharmacokinetic models that measure the time course of the presence of gadolinium within the tissues.171, 175 Quantitative parameters commonly reported by pharmacokinetic model–based analysis include the extracellular volume fraction (ve) and the volume transfer constant (ktrans) between blood plasma and extravascular space, which is sensitive to blood flow rate and vessel permeability.160, 174 These parameters have been related to tumor angiogenesis and hypoxia by histological coregistration.158, 176, 177, 178 For example, in patients with head and neck tumors, several DCE‐MRI parameters, including ve and ktrans, were found to be correlated with positive staining for pimonidazole, a hypoxia tissue marker.177 These parameters can be clustered to identify regions with different blood flow dynamics.145, 179, 180 However, relevant to habitat imaging is that the lack of perfusion would be expected to correlate to hypoxia only in areas with relatively high cell density.

The above discussion suggests that tumor habitats can be identified using different mpMRI combinations.151, 181, 182 Preclinical studies are valuable because these habitats can be identified by coregistration with histology. An example is shown in Figure 4, which demonstrates that data from T2, T2* (T2star), DWI‐MRI, and DCE‐MRI images can be combined to identify 4 different habitats that reflect different underlying phenotypes. It is worth mentioning here that the image contrast on T2‐weighted MRI is determined mainly by tissue‐specific T2 relaxation times. Such images provide anatomic information and typically are used qualitatively and in mpMRI protocols. Although quantitative information also can be obtained from T2‐weighted images, T2 maps provide complete quantification because signals are captured during multiple echo times and then are fit to an exponential decay curve on a pixel‐by‐pixel basis, yielding quantitative T2 information. The habitat map derived from mpMRI in Figure 4 shows good visual agreement with the corresponding histological slice. In this example, a 3D‐printed tumor mold was used to facilitate the coregistration of the histological slices with MRI.

Figure 4

Multiparametric magnetic resonance imaging (mpMRI) data used to identify tumor subregions (“habitats”) in a 4T1 breast cancer allograft in a mouse model. (A) Parameters derived from T2 map, T2* map, diffusion‐weighted (DWI)‐MRI, and dynamic contrast‐enhanced (DCE)‐MRI were clustered to create a habitat map, demonstrating 4 intratumor subpopulations with different patterns of cellularity. Blue indicates necrotic; yellow, nonviable hypoxic; pink, viable hypoxic; green, viable nonhypoxic. (B) A coregistered habitat map derived from histology demonstrating good spatial correspondence with the habitat map derived from MRI. Histological slices were obtained by cutting the tumors using a 3‐dimensional printed mold, which ensured the coregistration with MRI slices. Hematoxylin and eosin (H&E) and immunohistochemistry images were automatically segmented and superimposed to create habitat maps derived from histology, showing necrosis (blue), viable normoxic (green), pimonidazole in a viable region (pink), and pimonidazole in a nonviable region (yellow).

In other examples, k‐means clustering of the ADC from DWI‐MRI, T2, and M0 (proton density) images from colorectal tumor xenografts were used to follow the temporal evolution of the tumor populations in response to therapy.181 In particular, 2 necrotic clusters were identified, both demonstrating high ADC values but differing by short or long T2 values. Necrotic regions with a long T2 value and high ADC value were confirmed as acellular, whereas the short T2 value in the other necrotic regions was ascribed to the presence of deoxyhemoglobin due to either hemorrhage or leaky vasculature. Low ADC values were found in areas of viable tumor cells and subcutaneous adipose tissue; however, viable tumor tissue exhibited lower M0 values compared with adipose. In a similar analysis, viable and necrotic areas were identified in primary tumors and metastases of a 4T1 syngeneic mouse model of breast cancer.183 Based on the k‐means clustering of ADC and T2 data, 4 clusters (2 viable and 2 necrotic areas) were identified and confirmed by histological hematoxylin and eosin staining.182 Necrotic regions were characterized by high ADC values, but could be clustered into 2 different habitats based on their quantitative T2 values. One viable region (V1) presented with low ADC and qT2 values, whereas the other (V2) demonstrated intermediate ADC and lower qT2 values. Immunostaining of hypoxia‐inducible factor 1 differentiated the 2 viable regions as well oxygenated or hypoxic subpopulations. Expression of hypoxia‐inducible factor 1 was concentrated in the perinecrotic regions, classified previously as V2, which are expected to be hypoxic. The evaluation of these different viable tumor populations was informative for treatment response. Viable subregions classified as V1, which presumably are well oxygenated, demonstrated regression after irradiation, whereas the V2 subregions, characterized as hypoxic, increased, suggesting a selective survival of a subpopulation of cells resistant to radiotherapy.151 In another study, DCE‐MRI parameters were evaluated over the viable and necrotic subpopulations previously classified by the multispectral analysis of DWI‐MRI, T2, and proton density. After an antiangiogenic therapy (antivascular endothelial growth factor antibody G6‐31), a reduction in the Ktrans values was observed in the viable tumor regions, whereas Ktrans did not change in the areas classified as necrosis, and only provided noise to the DCE estimates.184

Clinically, mpMRI approaches are being used to identify tumor habitats, which can be used to predict outcomes for patients with cancer. In 36 patients with soft tissue sarcoma who were treated with neoadjuvant doxorubicin, the T1 postcontrast and T2 images each were divided into “high” and “low” using Otsu segmentation, and then combined to generate 4 distinct habitats: high‐high, high‐low, low‐high, and low‐low (Fig. 5).185, 186 The goal of the Otsu thresholding method is to choose the threshold that minimizes the intraclass variance. The volumes of each of these regions were entered along with multiple clinical features in a meta‐analysis to predict progression‐free survival (PFS) and overall survival (OS). In this, only 1 univariate predictor, the volume of the (high T1 postcontrast/low T2) habitat, was able to accurately predict OS and PFS. In cases in which this habitat occupied >18% of the tumor volume, it predicted a short PFS (P = .05) and lower OS (P = .036). In addition, these habitats could predict the existence of occult metastases with >80% accuracy. Another important aspect of these studies is that soft tissue sarcomas are a histologically heterogeneous group of diseases, but this classifier was applicable to all, regardless of histological subtype.

Figure 5

Multiparametric magnetic resonance imaging (mpMRI) habitats in patients with soft tissue sarcoma. Otsu thresholding was used to separate postcontrast T1 images into low (L) and high (H) intensities. Within each of these subsets, Otsu thresholding again was used to separate quantitative T2 images into L and H, generating 4 habitats: L‐H (indicated in green), L‐L (indicated in blue), H‐L (indicated in yellow), and H‐H (indicated in red). The occurrence of the yellow habitat was the only significant predictor of outcome in univariate analyses (in cases in which this habitat occupied >18% of the tumor volume, progression‐free survival and overall survival were significantly worse).

In patients with glioblastoma, habitat imaging can be used to segment normal from abnormal tissue. As shown in Figure 6, ADC, precontrast T1 and postcontrast T1, T2‐weighted, and FLAIR images can be combined and clustered to classify tissue classes.186, 187 It has been proposed that such classifiers can be used instead of segmentation, possibly mitigating the problems discussed previously. To predict outcomes, habitat images have been generated based on postcontrast T1 and FLAIR sequences.188 First, tumors were separated into regions of low or high contrast enhancement, based on the postcontrast T1‐weighted data. Then, FLAIR values were added, thus classifying the images into 5 habitats according to different combinations of blood flow and cellularity. In this study, patients initially were divided into short‐term or long‐term survival groups. Long‐term survivors had a predominant tumor habitat with high enhancement and intermediate cell density, whereas the tumors of short‐term survivors had habitats with low enhancement and high FLAIR signals (indicating necrosis). It is interesting to note that habitats with low blood flow also exhibited cellularity comparable to that of habitats noted in regions with high blood flow, which could indicate the presence of viable cells that had adapted to survive in hypoxic and acidic conditions. More recently, these same authors presented a computational framework with which to quantitatively extract image features from the tumor habitats that demonstrated high accuracy to potentially predict the survival of patients with glioblastoma.189 In related studies from another group, habitats in patients with glioblastoma multiforme were defined using postcontrast T1 and FLAIR data, and the spatial diversity features obtained from the habitats were associated with survival190 as well as with the epidermal growth factor receptor expression status in patients with glioblastoma.191

Figure 6

Five coregistered magnetic resonance imaging (MRI) sequences, namely apparent diffusion coefficient, T1 precontrast, T1 postcontrast, T2 weighted, and fluid‐attenuated inversion recovery (FLAIR), were combined to generate this map containing 7 habitats. The entire brain was clustered to identify normal and abnormal tissue classes. Note that this provides an automatic delineation of the tumor (arrow) along with its accompanying subregions.

An emerging area of interest considers the use of PET data in combination with CT or MRI to identify habitats with distinct combinations of metabolic activity (PET) with cellularity or perfusion (MRI). In patients with lung, colorectal, or breast carcinomas, subregions based on metabolic rate of glucose, which can be calculated from the pharmacokinetic analysis of dynamic PET images and is used for measuring local glucose consumption, and CT revealed a relationship between high or low metabolic regions with blood volume, cellular uptake, and washout as well as phosphorylation rate constants.192 In a novel application of PET‐MRI, 13 patients with primary or metastatic cancers were examined with DCE‐MRI and PET using [18F]‐galacto‐arginine‐glycine‐aspartic acid (18F‐galacto‐RGD) (which binds to αvβ3, an integrin upregulated during angiogenesis) and [F]fludeoxyglucose (18F‐FDG) (glycolysis).193 The total tumor area was divided into 4 habitats with high and low RGD and FDG uptake, respectively. Perfusion (with DCE‐MRI) was noted to be highest in tumor areas with simultaneously high RGD and FDG uptake, and was restricted in areas with low uptake of both RGD and FDG. In a recent PET‐MRI study of preclinical and clinical breast cancers, Gaussian mixture models combining FDG‐PET and DCE‐MRI were found to be able to distinguish between cystic hyperplastic, solid acinar, and solid nodular malignancies. These subregions were validated with coregistered hematoxylin and eosin staining and ex vivo autoradiography in a preclinical setting.194

Conclusions

In this era of personalized medicine in oncology, we have a responsibility to collect as much meaningful information from different modalities as possible, which can help to make better informed decisions. Although mutational analysis and gene expression data are indisputably important to this process, there are additional relevant data that can be used. Quantitative imaging is able to contribute significantly to decision support for 3 major reasons: 1) virtually every patient with cancer is imaged with CT, MRI, and/or PET; 2) these images are obtained from the entire tumor, along with metastases, and thus can be used to describe and classify heterogeneity; and 3) these images can be obtained routinely longitudinally to monitor responses and to guide specific therapies. As the technologies of radiomics mature, including conventional radiomics, DL, and their combination, we expect that the classifications will become more granular and will be able to identify subtumor regions (“habitats”) with specific compositions of cancer cell phenotypes and genotypes, along with immune and stromal involvement, and thereby provide better options for diagnosis, staging, treatment, and monitoring of response.

Funding Support

Supported by National Institutes of Health grants R01CA187532 (to Robert J. Gillies), U01CA143062 (to Robert J. Gillies), U01CA187947 (to Sandy Napel), R01CA160251 (to Sandy Napel), U24CA194354 (to Hugo J.W.L. Aerts), and U01CA190234 (to Hugo J.W.L. Aerts).

Conflict of Interest Disclosures

Sandy Napel has acted as a paid member of the medical advisory boards for Fovia Inc, Radlogics Inc, and EchoPixel Inc and has acted as a paid consultant for Carestream Inc for work performed outside of the current study. In addition, Dr. Napel has a patent issued (US 9,721,340). Robert J. Gillies has received a grant from the National Cancer Institute for work performed as part of the current study and is an investor and board member of HealthMyne and has patents issued for Systems, Methods and Devices for Analyzing Quantitative Information Obtained from Radiological Images (PCT 9,721,340; issued August 1, 2017) and Systems and Methods for Diagnosing Tumors in a Subject By Performing a Quantitative Analysis of Texture‐Based Features of a Tumor Object in a Radiological Image (PCT 9,940,709; issued April 10, 2018), and a patent pending for Radiologically Identified Tumor Habitats (15/124,718).