Integrated Multiparametric Radiomics and Informatics System for Characterizing Breast Tumor Characteristics with the OncotypeDX Gene Assay

Simple Summary Artificial Intelligence methods using machine learning and radiomics is an emerging area of research for radiological and oncological applications for patient management. Recent evidence from breast cancer suggests that different breast cancer subtypes may respond differently to adjuvant therapies. The use of a 21-gene array assay called OncotypeDX can predict potential recurrence of cancer in patients with estrogen positive breast cancer after treatment, however, there are potential cost disadvantages that hamper its widespread use. Multiparametric magnetic resonance imaging can simultaneously identify key functional parameters and provide unique imaging phenotypes of breast cancer, which is used in radiomic analysis. Radiomics provide quantitative information of different tissue types. We have developed a new machine learning radiomic informatics tool that integrates clinical and imaging variables, single, and multiparametric radiomics to compare with the OncotypeDX test to stratify patients into three risk groups: low, medium, and high risk of breast cancer recurrence. Abstract Optimal use of multiparametric magnetic resonance imaging (mpMRI) can identify key MRI parameters and provide unique tissue signatures defining phenotypes of breast cancer. We have developed and implemented a new machine-learning informatic system, termed Informatics Radiomics Integration System (IRIS) that integrates clinical variables, derived from imaging and electronic medical health records (EHR) with multiparametric radiomics (mpRad) for identifying potential risk of local or systemic recurrence in breast cancer patients. We tested the model in patients (n = 80) who had Estrogen Receptor positive disease and underwent OncotypeDX gene testing, radiomic analysis, and breast mpMRI. The IRIS method was trained using the mpMRI, clinical, pathologic, and radiomic descriptors for prediction of the OncotypeDX risk score. The trained mpRad IRIS model had a 95% and specificity was 83% with an Area Under the Curve (AUC) of 0.89 for classifying low risk patients from the intermediate and high-risk groups. The lesion size was larger for the high-risk group (2.9 ± 1.7 mm) and lower for both low risk (1.9 ± 1.3 mm) and intermediate risk (1.7 ± 1.4 mm) groups. The lesion apparent diffusion coefficient (ADC) map values for high- and intermediate-risk groups were significantly (p < 0.05) lower than the low-risk group (1.14 vs. 1.49 × 10−3 mm2/s). These initial studies provide deeper insight into the clinical, pathological, quantitative imaging, and radiomic features, and provide the foundation to relate these features to the assessment of treatment response for improved personalized medicine.


Introduction
Integrating clinical health information with radiological imaging and other biomarkers could be beneficial for different types of cancer. This integration of seemingly disparate data may improve our understanding of the complex nature of cancer and potentially provide predictive markers with clinical benefit in certain cancer phenotypes. For example, in breast cancer, there is active research on how to predict the potential of local recurrence after conservative treatment. One such method to predict local and systemic recurrence of breast cancer is the OncotypeDX assay [1][2][3]. OncotypeDX is based on the mRNA expression by RT-PCR in estrogen receptor (ER) positive disease without the human growth factor receptor 2 (HER2-neu) overexpression and is tested on tissue obtained at biopsy or core diagnostic or surgical samples [4,5]. OncotypeDX has been validated in prospective studies as a prognostic tool predictive of excellent outcomes in patients with ER-positive disease treated with endocrine therapy [1]. It has since shown to be a predictive tool to identify patients with breast cancer, most likely to benefit from the addition of adjuvant chemotherapy to endocrine therapy [1][2][3][4][6][7][8][9][10] Multiparametric (mp) radiological imaging can accurately detect and characterize breast lesions using advanced quantitative parameters [11][12][13][14]. Dynamic contrast-enhanced (DCE)-magnetic resonance imaging (MRI), which is a marker of the vascularity and permeability of breast lesions, can characterize malignant lesions by the rapid uptake of a contrast agent, followed by fast washout over time, and benign lesions by slow uptake and persistent or plateau washout [11,[15][16][17][18][19]. Moreover, by interrogating the movement of water within the intra-and inter-cellular environments of normal and lesion tissue using diffusion-weighted imaging (DWI), with the apparent diffusion coefficient (ADC) of water map, can provide characterization of breast and other lesions [12,[20][21][22][23][24]. Recent research had demonstrated that textural and shape analysis of multiparametric radiomics (mpRad) could potentially add a unique insight into the underlying tissue pathology [25]. Radiomics is an area of research, which deals with quantitative and qualitative textural analysis of the underlying tissue pathology [26][27][28][29][30][31][32][33]. Radiomic analysis of breast mpMRI has been investigated for breast cancer diagnosis and treatment response assessment in the research setting across multiple studies [25][26][27][28]31,32,34]. However, to date, no one has investigated using mpRad coupled with breast mpMRI and informatic data and testing with OncotypeDX for breast lesion characterization.
The challenge is to accurately combine mpMRI with radiomic, clinical, and pathologic features to stratify patients and identify the potential for cancer recurrence, similar to the OncotypeDX. To answer this challenge, we have developed a new machine-learning informatic method termed Integrated Radiomics Informatic System (IRIS), which can be applied to multiparametric MRI and radiomics, clinical, and pathologic descriptors, as well as a gene array analysis [32,[35][36][37][38]. The purpose of this study is to test the IRIS algorithm by combining data from imaging, radiomics, and electronic medical health records (EHR) to stratify patients into three risk groups: low, medium, and high risk, and then compare these groups with the OncotypeDX 21-gene assay scores.

Clinical Subjects
All studies were performed in accordance with the institutional guidelines for clinical research under a protocol approved by our Institutional Review Board (IRB number: NA_00001113 and NA_00022703), and all Health Insurance Portability and Accountability Act (HIPAA) agreements were followed for this retrospective study and informed consent was waived. All procedures performed in studies involving human participants were in accordance with the ethical standards of the institutional and/or national research committee and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. Patients were selected from the Johns Hopkins Integrated Breast Cancer Research Database, developed by one on the authors (C.U.), and who had underwent MRI as part of the clinical health record review. Of the patients with biopsy proven breast cancer who presented to our facility for bilateral breast MRI, 123 patients were identified to have both the OncotypeDX and an advanced MRI exam, which included DCE and DWI. Our inclusion criteria were: (1) breast imaging on 3T MRI scanner, (2) dynamic contrast enhanced and DWI MRI sequences, and (3) pathology proven diagnosis of Estrogen receptor positive (ER+) breast cancer and (4) OncotypeDX having been performed on lesion tissue samples. There were 80 patients with 83 lesions that satisfied the inclusion criteria.

Histological Phenotyping
All breast cancers were categorized by histological phenotyping based upon immunohistochemistry (IHC). Estrogen and progesterone receptors (ER and PR), HER2-neu by Fluorescence in situ hybridization (FISH), and Ki-67 proliferation index (%). The Elston tumor grades for each lesion were distributed as Grade 1 (9%), Grade 2 (78%), and Grade 3 (13%). Histopathological data was obtained from the breast pathology database (C.B.U.). All patient demographics matched the current clinical criteria for the OncotypeDX test (ER+).

Multi-Subspace Embedding and Clustering
We have developed a machine learning model that integrates different types of clinical and imaging parameters, which allows for the construction of a clinical decision support model [36,38,39]. The inherent high dimensionality of the clinical and imaging parameters and any complex correlations within the data presents significant challenges for integration and visualization of the data. These challenges were solved using the nonlinear dimensionality reduction (NLDR) method [39]. Nonlinear dimensionality reduction algorithms transform and embed a D dimensional space into a lower d dimensional manifold representation of D's intrinsic dimensionality, where d < D. The goal of the multi-subspace embedding and clustering method is to transform the patient space, represented as X = x 1 , x 2 , . . . , x n p ∈ R D where, x i represents the ith patient, n p represents the number of patients, and D represents the number of clinical and imaging parameters, into an Integrated Radiomics Informatics System (IRIS) visualized by a heatmap as shown in Figure 1. The steps of the IRIS system are outlined below [32,36].
In the first step, n p d-dimensional subspaces, where d ∈ {1, 2, 3}, are extracted from each D-dimensional patient vector, x i , by selecting all combinations of one, two, or three IRIS parameters from D. Here, n p is given by (D, 3). Each subspace can be represented as in Equation (1): The second step involves transformation of each d-dimensional subspace, S i ∈ R d ∀ i ∈ 1, 2, . . . , n p into a one dimensional embedding, Y i = y 1 , y 2 , . . . , y n p ∈ R 1 ∀i ∈ 1, 2, . . . , n p using a nonlinear dimensionality reduction algorithm [40]. In the third step, each one-dimensional embedding, Y i is Cancers 2020, 12, 2772 4 of 23 evaluated against a ground truth (e.g., OncotypeDX scale) based on the correlation coefficient, R i between Y i and the ground truth. In this paper, we employed the correlation coefficient as the evaluation metric. The aim of this step is to identify the optimal set of one-dimensional embeddings, U, such that, R i ≥ 0.5, as shown in the following Equation (2): In the fourth step, a hierarchical clustering algorithm is used to cluster the set, U to produce two clustering configurations. The first clustering configuration is between different informatics and radiomics parameters. The first clustering configuration provides a visualization of relationships between the clinical and imaging parameter embedding, each parameter's importance, and identify redundant embeddings. The second clustering configuration is between each patient, visualizing the relationship between different patients. This allows IRIS to identify patients clustered into different risk groups and classify any unknown patient into a relevant risk group. In the first step, np d-dimensional subspaces, where ∈ {1,2,3}, are extracted from each Ddimensional patient vector, , by selecting all combinations of one, two, or three IRIS parameters from D. Here, np is given by ∁ ( , 3). Each subspace can be represented as in Equation (1): The second step involves transformation of each d-dimensional subspace, ∈ ∀ ∈ {1,2, … , } into a one dimensional embedding, = { 1 , 2 , … , } ∈ 1 ∀ ∈ {1,2, … , } using a nonlinear dimensionality reduction algorithm [40]. In the third step, each one-dimensional embedding, Yi is evaluated against a ground truth (e.g., OncotypeDX scale) based on the correlation coefficient, Ri between Yi and the ground truth. In this paper, we employed the correlation coefficient as the evaluation metric. The aim of this step is to identify the optimal set of one-dimensional embeddings, U, such that, Ri ≥ 0.5, as shown in the following Equation (2):

Feature Importance
In IRIS, the importance of each clinical or imaging parameter, i is calculated as the percentage of embeddings in U that include parameter i. The contribution of each parameter to the IRIS model allows for assessment of which parameters to keep and which to discard.

Complex Network Analysis of Informatics Parameters
Using IRIS, the high dimensional relationship between different clinical and imaging parameters was explored by modeling and analyzing a complex informatics network. Before modeling the complex network, the raw values corresponding to each clinical and imaging parameter were transformed into Cancers 2020, 12, 2772 5 of 23 risk prediction score normalized in the range 0-1, such that zero corresponds to low risk and one corresponds to high risk. The normalized risk prediction score for every parameter was calculated using the following steps: First, a correlation coefficient, r i between the values, y i spanned by the parameter, i across all the patients and corresponding OncotypeDX scores are calculated. b.
Second, the range of clinical and imaging parameter values across all the patients are normalized from zero to one according to the following formula (3): Here, z i represents the resultant risk prediction score for each parameter, i.

Construction of Network Model of Informatics Parameters
The complex informatics network, G was constructed from the multidimensional data, Zi = {z 1 , z 2 , . . . , z K } ∈ R n p , where, z i represents the risk prediction score of the clinical and imaging parameter, i; K is the number of clinical and imaging parameters; and n p is the number of patients [38]. The network G is represented as G = (V,E) with V = {v 1 ,v 2 , . . . ,v K } being the set of K vertices representing the clinical and imaging parameters and E being the adjacency matrix indicating the interactions between the clinical and imaging parameters in the form of edge weights. The edge weight between any two vertices v i and v j was computed by Equation (4) as follows The edge weights represent the distances between any two parameters. We defined a neighborhood parameter, k, to define the number of nearest neighbors each parameter could be connected to. The connectivity in the resulting complex network is dependent on the value of k chosen. If the value of k is chosen to be too large, the complex network may produce short circuit or spurious connections while a low value of k would produce a disconnected network [40]. The value of k selected as three by empirical analysis.

Statistics and Topological Characteristics of the Complex Network
The complex network was analyzed using graph summary metrics and centrality metrics [41,42] The average path length and diameter are the basic statistical metrics computed for any complex network. A path is defined as the set of edges connecting any two nodes and the sum of weights of these edges represent the path length. Average path length, as the name suggests, is the average of the path lengths across all pairs of nodes or clinical and imaging parameters. Diameter is the maximum value among all the path lengths. Graph centrality metrics identify the most important clinical and imaging parameters in the complex network and called the hub nodes [43,44]. These hub nodes influence the network properties. Furthermore, the probability of any incoming node connecting to these hub nodes is significantly higher than connecting to other nodes [42]. The hub nodes may correspond to the key clinical and imaging parameters that are predictors of breast cancer recurrence risk. In total, the following metrics were extracted from the complex network: degree distribution, average path length, diameter, clustering coefficient, and different centrality measures such as degree centrality, harmonic centrality, and betweenness centrality (see Appendix A). We define a tissue signature (TS) that represents the composite feature representation of a tissue type based each of the different imaging sequences [25]. Mathematically, for N different imaging parameters with TS at a voxel position, S p is defined as a vector of gray level intensity values at that voxel position, p across all the (N) images in the data sequence for different tissue types and is given by the Equation (5), where, I p is the intensity at voxel position, p on each image, and T is the transpose. Then we define the tissue signature probability matrix (TSPM) as an N dimensional matrix with each cell representing a tissue signature configuration. The TSPM characterizes the spatial distribution of tissue signatures within a Region of Interest (ROI). From the TSPM, we derive three radiomic features, the TSPM entropy, uniformity, and mutual information (MI). The tissue signature first order statistics (TSFOS) features, which characterize the distribution of voxel intensities across all the imaging parameters are calculated, similar to traditional first order radiomics. The second order mpRad features are calculated from the tissue signature co-occurrence matrix features (TSCM). The TSCM characterizes the spatial relationship between tissue signatures within an image or ROI.
The input parameters were determined from empirical analysis for the statistical texture analysis methods of FOS, GLCM, GLRLM, and NGTDM were set as follows: a.
The distance d for GLCM was set to one voxel. d.
Both GLCM and GLRLM were evaluated in all the four directions-0 • , 45 • , 90 • , and 135 • . Rotational invariance was achieved by extracting the radiomic features from GLCM and GLRLM averaged across all the directions.
We also computed radiomic feature maps (RFMs) [32] corresponding to the FOS and GLCM features. The mean of the RFM feature map metrics were used in the prediction model. Finally, the fractal dimension for the tumor intensity profile and tumor boundary were evaluated in addition to the convexity of the tumor boundary. All of the radiomic features were extracted using our in-house software developed using MATLAB (Version 19a).

Patient Classification
We implemented patient classification using the hybrid IsoSVM feature transformation and classification algorithm [32] based on the Isomap [40] and the Support Vector Machine (SVM) algorithms [50]. We evaluated the following four models: a.
Low vs. intermediate risk group. b.
Low vs. high-risk group. c.
Intermediate vs. high-risk group. The imbalance in the number of patients in different risk groups was overcome by setting different misclassification penalties for different risk groups while training the SVM [51]. The optimal values for the Isomap neighborhood parameter and the misclassification penalty were estimated using leave one out cross validation.

Multiparametric Breast Imaging
Patients were scanned on a 3T MRI system (3T Achieva, Philips Medical Systems, Best, The Netherlands) using a bilateral, dedicated four-channel, phased array breast coil (Invivo, Orlando, FL, USA) with the patient in the prone position.

Pharmacokinetic Contrast Enhancement Metrics
Pharmacokinetic kinetic DCE MRI provides quantitative metrics of the volume transfer constant (K trans (min −1 )), which characterize uptake of the contrast agent, the leakage within the extracellular extravascular space (Ve (%)), and the transfer rate constant (kep (min −1 )). Post-processing of the DCE exam was performed by a combined Brix and Tofts model [15,54,55] using DynaCAD (Invivo, FL, USA) software from the identified breast lesions, and detailed in these manuscripts [11,[56][57][58].

ADC Mapping
Regions of Interest (ROI) were drawn on normal appearing glandular tissue and breast lesions defined by DCE MRI. Means and standard deviations were calculated for both tissue types. Ratios of lesion ADC to glandular tissue ADC (L/GT) were calculated from equation (6) below using the lesion and glandular tissue [23].

Statistical Analysis
We computed summary statistics (mean and standard deviation of the mean) for the quantitative imaging parameters from the mpMRI. An unpaired two-sided t-test was performed between each pair of risk groups imaging parameters to determine statistical significance. Sensitivity and specificity, and Area Under the Curve (AUC), were calculated to determine the classification of the patients in the different groups. Statistical significance was set at p ≤ 0.05.

Radiological Findings
Representative multiparametric breast imaging for each risk group as defined by OncotypeDX are illustrated in Figure 2. The high-risk group had the largest tumor size (2.9 ± 1.7 mm). Followed by the low-risk group tumor size (1.9 ± 1.3 mm) and the intermediate risk group (1.7 ± 1.4 mm). These differences were not significant (p > 0.5). For advanced MRI parameters, there were clear differences in each parameter and Oncotype risk groups. The PK-DCE parameter K trans values for the intermediate-and high-risk groups were higher (0.46 and 0.49 (1/min) (p = 0.26)) compared to the low-risk group (0.30 (1/min) (p = 0.02). Similar results were noted for the other PK-DCE parameters (extracellular extravascular space (EVF) and k ep ). The maximum contrast enhancement from DCE was largest for the low-risk group (503 ± 33 s), compared to the intermediate-risk (461 ± 24 s), and high-risk groups (468 ± 31 s). Similarly, the ADC map values from the high-and intermediate-risk patients in the lesion tissue were significantly lower (p < 0.05) than those for the low-risk patients (1.14 vs. 1.49 × 10 −3 mm 2 /s). However, the ADC map values in glandular tissue remained constant across all groups (2.14-2.17 × 10 −3 mm 2 /s). The bar graphs are shown in Figure 3.

Radiological Findings
Representative multiparametric breast imaging for each risk group as defined by OncotypeDX are illustrated in Figure 2. The high-risk group had the largest tumor size (2.9 ± 1.7 mm). Followed by the low-risk group tumor size (1.9 ± 1.3 mm) and the intermediate risk group (1.7 ± 1.4 mm). These differences were not significant (p > 0.5). For advanced MRI parameters, there were clear differences in each parameter and Oncotype risk groups. The PK-DCE parameter K trans values for the intermediate-and high-risk groups were higher (0.46 and 0.49 (1/min) (p = 0.26)) compared to the low-risk group (0.30 (1/min) (p = 0.02). Similar results were noted for the other PK-DCE parameters (extracellular extravascular space (EVF) and kep). The maximum contrast enhancement from DCE was largest for the low-risk group (503 ± 33 s), compared to the intermediate-risk (461 ± 24 s), and high-risk groups (468 ± 31 s). Similarly, the ADC map values from the high-and intermediate-risk patients in the lesion tissue were significantly lower (p < 0.05) than those for the low-risk patients (1.14 vs. 1.49 × 10 −3 mm 2 /s). However, the ADC map values in glandular tissue remained constant across all groups (2.14-2.17 × 10 −3 mm 2 /s). The bar graphs are shown in Figure 3.   . Bar graphs of quantitative multiparametric MRI parameters for each risk group from the OncotypeDX (OnDx). There are significant differences between each group of patients in the apparent diffusion coefficient (ADC (×10 −3 mm 2 /s)) of water and the Pharmacokinetic Dynamic Contrast Enhanced (PK-DCE) metrics. The PK-DCE metrics are the volume transfer constant (K trans (min −1 )) and the fractional volume of the extracellular extravascular space (EVF (Ve)). * p < 0.05.

Single and Multiparametric Radiomics
The single and multiparametric radiomic features demonstrated significant differences between the low, intermediate, and high-risk groups and shown in Figure 4. The single and multiparametric radiomic features, first order entropy feature map and the GLCM entropy feature map, which are measures of heterogeneity, were significantly higher for the low risk group as compared to the intermediate and high-risk groups. In contrast, the first order uniformity and GLCM energy, both of which are measures of homogeneity demonstrated an opposite trend. These radiomic results are summarized in Table 2.  . Bar graphs of quantitative multiparametric MRI parameters for each risk group from the OncotypeDX (OnDx). There are significant differences between each group of patients in the apparent diffusion coefficient (ADC (×10 −3 mm 2 /s)) of water and the Pharmacokinetic Dynamic Contrast Enhanced (PK-DCE) metrics. The PK-DCE metrics are the volume transfer constant (K trans (min −1 )) and the fractional volume of the extracellular extravascular space (EVF (V e )). * p < 0.05.

Single and Multiparametric Radiomics
The single and multiparametric radiomic features demonstrated significant differences between the low, intermediate, and high-risk groups and shown in Figure 4. The single and multiparametric radiomic features, first order entropy feature map and the GLCM entropy feature map, which are measures of heterogeneity, were significantly higher for the low risk group as compared to the intermediate and high-risk groups. In contrast, the first order uniformity and GLCM energy, both of which are measures of homogeneity demonstrated an opposite trend. These radiomic results are summarized in Table 2.
Cancers 2020, 11, x 10 of 23 Figure 3. Bar graphs of quantitative multiparametric MRI parameters for each risk group from the OncotypeDX (OnDx). There are significant differences between each group of patients in the apparent diffusion coefficient (ADC (×10 −3 mm 2 /s)) of water and the Pharmacokinetic Dynamic Contrast Enhanced (PK-DCE) metrics. The PK-DCE metrics are the volume transfer constant (K trans (min −1 )) and the fractional volume of the extracellular extravascular space (EVF (Ve)). * p < 0.05.

Single and Multiparametric Radiomics
The single and multiparametric radiomic features demonstrated significant differences between the low, intermediate, and high-risk groups and shown in Figure 4. The single and multiparametric radiomic features, first order entropy feature map and the GLCM entropy feature map, which are measures of heterogeneity, were significantly higher for the low risk group as compared to the intermediate and high-risk groups. In contrast, the first order uniformity and GLCM energy, both of which are measures of homogeneity demonstrated an opposite trend. These radiomic results are summarized in Table 2.

Integrated Radiomics Informatic System (IRIS) Model
The IRIS heatmap risk profile of each patient is shown in Figure 5. The top surrogates for mpRad and single radiomics, imaging, and histological parameters determined from the clinical and imaging model based on the IRIS heatmap are summarized in Tables 3 and 4.

Integrated Radiomics Informatic System (IRIS) Model
The IRIS heatmap risk profile of each patient is shown in Figure 5. The top surrogates for mpRad and single radiomics, imaging, and histological parameters determined from the clinical and imaging model based on the IRIS heatmap are summarized in Tables 3 and 4.      For the patient classification, the quantitative imaging parameters, ADC map values and the PK-DCE parameters were combined with histopathological parameter of Ki-67, and the radiomic parameters of FOS/GLCM entropy, FOS uniformity, and GLCM energy features. For the single parameter-based IRIS model, the IsoSVM model classified the low risk group from the intermediate and high-risk groups with a sensitivity and specificity of 95% and 88%, respectively, with an AUC of 0.93 (CI = 0.88-0.99). The optimal neighborhood parameter, dimensionality, and imbalance ratio were found at 60, 10, and 2:1, respectively. The ROC curves from different inter-group IsoSVM classifiers are shown in Figure 6. The lowest and non-diagnostic AUC (0.64 (CI = 0.47-0.81)) resulted from the comparison between the intermediate and high-risk groups. Figure 7 demonstrates the ROC curves from different inter-group IsoSVM classifiers using mpRad based IRIS model. There was no significant difference between the ROC curves from the two radiomic feature sets. The p values for comparison between the ROC curves for single and multiparametric radiomics were p = 0.33 (low vs. intermediate    The topological graph theoretic metrics for integrated centrality and other different centrality measures for each informatics parameter are summarized in Tables 4 and 5. The average path length between each parameter was 2.1 and diameter of the complex informatics network was 5.1. The average clustering coefficient was 0.53, much higher than the clustering coefficient of Erdos-Renyi random graph (CCER = 0.0228). Figures 8 and 9 illustrate the complex interaction network and the inter-parametric relationships between all the variables for single and multiparametric radiomics models. The hub IRIS parameters for single and multiparametric radiomics IRIS models have been summarized in Tables 5 and 6.

Discussion
We have introduced and demonstrated an advanced NLDR integrated clinical and imaging model (IRIS) to analyze the relationships and interactions between mpMRI parameters, radiomics, clinical heath records, and histological variables and compared these results with the OncotypeDX assay for risk assessment of breast cancer recurrence. The data integration by the IRIS model using radiomic, radiological, and clinical variables were able to group patients into the three different categories of low, intermediate, and high risk of breast cancer recurrence. Importantly, using IRIS, we defined several mpMRI and mpRad variables that were predictive of potential local and systemic tumor recurrence compared with categorization as defined by the OncotypeDX risk score. This integration of mpMRI, clinical data, and radiomics compared favorably to OncotypeDX and may lead to an accurate non-invasive assessment for risk of local and systemic recurrence. For example, we found that the most important radiological and histological parameters were the ADC map values, PK-DCE metrics, and Ki-67. These quantitative biological metrics reflect the cellularity, vascularity, and proliferation status of the tumor. Therefore. By using an integrated machine learning model of imaging, radiomics, clinical heath records, and histopathology data has the potential to describe unseen features of cancer and may provide data for precision personalized care as shown from the IRIS visualization heatmap. This is the one of the first studies to use an integrated graph theoretic and machine learning model based on quantitative mpMRI, clinical variables, and both single and multiparametric radiomics in breast cancer compared with gene array data [35,37].
The categorization of the different risk groups from our model was strikingly consistent based on the combined imaging, radiomics, and pathological variables. Indeed, the ADC map values were lower in the high and intermediate risk group consistent with the increased cell proliferation metric, Ki-67 from histological analysis. The PK-DCE parameters and lesion size demonstrated similar characteristics in these groups. Moreover, the radiomics parameters revealed an interesting underlying pattern across each of the three risk groups.
This report demonstrates that using an advanced unsupervised machine learning method in breast cancer and integration of several variables can accurately separate those cancers into different risk stratifications consistent with the OncotypeDX results. Finally, by using the complex interaction mapping, one can visualize the connections formed between each variable, which may form the basis for even further predictive modeling using the heatmap visualization tool or newer surveillance methods of recurrence in clinical use. For example, in some low risk patients, there were potential markers of increased risk, which may be missed by just looking at "single points" of data and these cases could be followed more closely over time. Similarly, potentially high-risk markers were seen in the intermediate risk cases, which could be useful for a surveillance program in these patients.
The clinical and radiological parameters utilized in this study were derived from our clinical experience, since current treatment decision algorithms are based on standard clinicopathologic prognostic, and predictive factors, large datasets using clinical measures such as tumor size, node status, grade, ER, and HER2-neu [9,10,[59][60][61][62][63]. Similarly, imaging features such as, breast density, lesion morphology, size, enhancement patterns as well as quantitative metrics (ADC map values and PK-DCE) are routinely used in practice and are familiar to the radiologist and readily available. Finally, radiomic feature analysis methods are increasing bringing a new potential quantitative biological biomarker to different cancers. The OncotypeDX assay has been shown to be a predictive tool to identify patients most likely to benefit from the addition of adjuvant chemotherapy to endocrine therapy with validation in prospective studies [1][2][3]. Thus, the ability to combine these quantitative measures would be an important step in ensuring that "the right patient receives the right treatment".
We developed an integrated informatics-radiomics decision support system (IRIS) based on multi-subspace embedding and clustering method for the purpose of diagnosis or prognosis [36]. Furthermore, the complete multi-subspace embedding, and clustering method is unsupervised and does not require any training data. The IRIS heatmap provides a visualization of relationship between different cancers along with quantifiable embedding metrics. Using the IRIS heatmap, we would be able to identify a patient or a group of patients with the most similar informatics embedding metrics and use these metrics for a new patient with an unknown risk of recurrence. Understanding the complex relationships between different embeddings can provide an insight on how these metrics are related at biological level predicting recurrence of breast cancer. Interestingly, the lesion size was not an accurate feature for categorizing cancers into the risk groups in this study. The lesion size was largest for the highest risk group, but was smallest for the intermediate risk group, suggesting that lesion size alone is not an accurate predictor. Differentiating and characterizing benign from invasive breast cancer is an important issue that was the focus of many different studies [12,19,[63][64][65][66].
The complex interaction network provided insight into how the different IRIS parameters relate to each other. For example, radiomic features of entropy and uniformity have higher integrated centrality in the subgraph indicating these features provided complementary information about the underlying network. In contrast, glandular ADC and Ki-67 have lower integrated centrality and more associated with other distinct features in the network. The highest integrated centrality measures were the mpRad and vascular features suggested the mpMRI provides more information about breast lesions.
Prior studies developing risk scores comparable to the OncotypeDX, typically used a single or just few MRI parameters [10,34,67,68]. Some differences with those studies and the present study are several; we employed machine learning and graph theory methods to differentiate these entities. The IRIS model incorporated pathophysiological and imaging characteristics of different breast tissue types, enabling a more predictive model of the tumor environment. Finally, we used the standard of care BI-RADS information from the radiologist report in clinical records, which is an informatic risk assessment and provides diagnostic characteristics of the lesion and surrounding tissue. The addition of BI-RADS increased the dimensionality to our model by the combination of the breast density, mass shape, margins, enhancement patterns, and other factors, in conjunction with the radiomics features, a more complete picture of the lesion and surrounding texture characteristics. One limitation to the study is our small sample size of patients, but this report provides encouraging preliminary data for further testing on a larger cohort in a prospective study. For example, the non-diagnostic AUC between the intermediate and high-risk group in this study may be attributed to the large class imbalance of one to seven between the groups. Recent reports have shown non-diagnostic AUCs (0.50-0.77) results [69] in a larger group of similar patients studied using the OncotypeDX and single parameter radiomics consistent with our results. However, the authors in those studies used only a single MRI parameter (DCE) and developed some different features based on image analysis, but they did not utilize quantitative PK-DCE, ADC mapping, or other common MRI parameters, which are more representative of the tumor microenvironment.

Conclusions
The incorporation of multiparametric radiomics into the IRIS model provided more quantitative metrics for better characterization and complete picture of breast lesions. Ongoing work is underway to see if IRIS could be used as a potential predictive or selection tool for patients for adjuvant treatment and assist in decision making with pathology for those tough cases in the "low to intermediate" risk group. IRIS could be used to identify high-risk features in the low to intermediate groups for potential active surveillance, in terms of risk defined by the integration of several variable over time.
In conclusion, these initial studies provide insight into the molecular underpinning of the surrogate imaging and clinical features and provide the foundation to relate these changes to histologic and molecular pathology parameters. The integration of these clinical and imaging parameters may help refine available prognostic and predictive markers, and improve clinical decision-making.  The clustering coefficient defines the connectedness of the neighborhood of an informatics parameter. Clustering coefficient ranges from zero to one with zero representing completely disconnected neighborhood and one representing completely connected neighborhood. Mathematically, the clustering coefficient, i is defined by Equation (A1),

Patents
Here e i is the number of connected edges neighborhood of i, k i are the number of nodes in the neighborhood of i making k i (k i −1) 2 the maximum possible number of edges in the neighborhood.

Degree Distribution
The degree distribution is the most fundamental metric calculated for any complex network. The degree distribution, P(k) is determined using the following equation Here, deg(i) represents the degree of the parameter, i, which is defined as the total number of parameters that are directly connected to it; N is the number of parameters.
The degree distribution enables us to identify whether the network is a scale-free network or not. For scale free networks, the degree distribution follows the power law P(k) ∝ k −γ , with the value of γ ranging between 2 and 3. Scale free networks are characterized by the presence of few highly connected hub nodes that influence the network properties and may correspond to the key parameters that are predictors of breast cancer recurrence risk [42].

Appendix A.2. Centrality Measures
The Centrality measures determine the importance of each informatic parameter in the complex network. The most widely used measures of centrality are betweenness, harmonic, and degree centrality [43,44,70]. We use three centrality measures, each centrality measure highlighting a different importance property of the network. Appendix A.2.1. Betweenness Centrality Betweenness centrality quantifies the amount of information that flows through each parameter or node [43]. It is defined as the number of shortest paths that pass through a parameter, given by Equation (A3): Here, N st is the total number of shortest paths between the parameters, s and t and N st (i) is the total number of shortest paths between s and t that pass through i.

Appendix A.2.2. Harmonic Centrality
Harmonic centrality of a parameter is defined as the sum of inverse of all geodesic distances (path lengths) from that parameter to all the other parameters [44]. Harmonic centrality measures the nodes that are influential and how information spreads on the local network. Mathematically, harmonic centrality of a parameter, i is given by Equation (A4): Here G(s,i) is the geodesic distance or path length between the parameters s and i.

Appendix A.2.3. Degree Centrality
Degree centrality or degree is defined as the total number of parameters linked to a node, i.e., the number of neighbors, where each parameter is connected in the complex network.

Appendix A.2.4. Integrated Centrality
Each centrality measure signifies the importance of each parameter based on a pre-define characteristic. Where, the significance of each parameter across all centrality indices can be defined using integrated centrality [17] as shown in Equation (A5):