Multi-kernel graph embedding for detection, Gleason grading of prostate cancer via MRI/MRS

2.1. DSS strategies for CaP detection via T2w MRI

Madabhushi et al. (2005) presented a supervised DSS system for detection of CaP from 4 Tesla (T) ex vivo prostate T2w MRI where 33 3D texture features (statistical, gradient, and Gabor) were quantitatively extracted at each voxel (T2w MRI spatial resolution) These extracted features were then used to train a number of supervised classifiers (Adaboost, Bayes, and Decision Trees) which were employed to assign a probability of CaP presence at each image voxel. Viswanath et al. (2012) developed a classifier strategy to obtain distinct quantitative imaging signatures for central and peripheral CaP on 25 patients who underwent T2w endorectal in vivo MRI.

2.2. DSS strategies for CaP detection via MRS

Kelm et al. (2007) presented a comparative study of classification techniques for prostate MRS data based on pattern recognition methods such as Principal Component Analysis (PCA) and Independent Component Analysis against quantification based feature extraction methods using Support Vector Machines, Random Forests, and Gaussian process classifiers. Tiwari et al. (2009) presented a DSS for CaP detection using 1.5 Tesla in vivo prostate MRS where each prostate spectrum was classified, on a per voxel basis, as either belonging to cancerous or non-cancerous classes using a hierarchical, clustering scheme in conjunction with non-linear dimensionality reduction (NLDR) methods. NLDR schemes were employed to obtain a low dimensional representation of high dimensional MR spectra, followed by hierarchical k-means clustering to identify CaP locations in the prostate.

2.3. DSS strategies for CaP detection via MP-MRI channels

In (Chan et al., 2003), a statistical classifier which integrated texture features from 1.5 T in vivo MP MRI was presented to generate a statistical probability map representing likelihoods of cancer for different regions within the prostate. Liu et al. (2009) examined in vivo MP MRI maps (T2w, DCE, DWI) within a fuzzy Markov Random Field (MRF) framework. The maps were generated via curve fitting of data from each of the protocols with the region of interest (ROI) limited to the peripheral zone of the prostate, while the evaluation of the results was done against manually delineated CaP regions on MRI (with corresponding whole-mount histology and ex vivo MRI data used for reference). Ampeliotis et al. (2008) reported a statistically significant improvement in detection accuracy in the combined use of image intensity features from both DCE and T2w MRI data for the classification of CaP, compared to the use of individual modalities. A comparison of different supervised and unsupervised methods for CaP segmentation using MP MRI was presented in (Ozer et al., 2010). However, most of these methods involve integrating imaging-imaging protocols, and it is not apparent that these schemes can be extended to accommodate imaging and non-imaging parameters (e.g. MRS), especially when the different channels to be combined reside in different dimensional spaces.

Recently, Tiwari et al. (2012) presented a DSS method combining information from the T2w MRI and MRS channels for CaP diagnosis via a combination of wavelets and PCA. However, PCA which is a linear DR scheme is typically unable to handle situations where the overlapping classes are embedded on non-linear manifolds – a situation that is more amenable to multi-kernel learning (MKL) (Lanckriet et al., 2004) and graph embedding (GE) (Shi and Malik, 2000) schemes. Kernels (Scholkopf and Smola, 2001) are positive definite functions which transform the input data to a dot product similarity space such that K(F(c_i), F(c_j)) = 〈Φ ((F(c_i)),Φ(F(c_j)))〉, where Φ is the implicit pairwise embedding between input feature vectors F(c_i) and F(c_j) associated with points c_i, c_j, and 〈.〉 denotes the dot product operation. MKL (Lanckriet et al., 2004) involves computing similarity matrices for kernels derived from the individual modalities being combined, so that fused classifiers (within the fused kernel space) can be built. Lanckriet et al. (2004) transformed data from amino acid sequences, protein complex data, gene expression data, and protein interactions into a common kernel space. The kernels were linearly combined and used to train a SVM classifier for classifying functions of yeast proteins. However, when a large amount of information is present in each input source, most COD methods, including MKL, suffer from the curse of dimensionality (Duda et al., 2000).

To avoid the curse of dimensionality, several DR schemes have been proposed (Roweis and Saul, 2000, Tenenbaum et al., 2000, Shi and Malik, 2000). One such DR strategy is GE (Shi and Malik, 2000), which accounts for the non-linearities in the data by constructing a similarity graph G = (V, W), where V corresponds to the set of vertices that connect pairwise points and W is a n × n weight matrix, reflecting the feature similarity between the n data points. The similarity graph is then reduced to lower dimensions by solving a simple Eigenvalue decomposition problem. However GE, like most other DR schemes, is unsupervised and often leads to overlapping clusters resulting in poor class discriminability. Recently several semi-supervised DR (SSDR) schemes have been proposed where some class label information is employed in the construction of a pairwise similarity matrix. Sugiyama et al. (2008) applied semi supervised-learning (SSL) to Fisher’s discriminant analysis in order to find projections that maximize class separation. Sun and Zhang (2009) implemented a semi-supervised version of PCA by exploiting between-class and within-class scatter matrices. Semi-supervised GE (Zhao, 2006) involves constructing a weight matrix (W) by leveraging the known class labels such that higher weights are given to within-class points and lower weights to points from different classes. The proximity of labeled and unlabeled data is then used to construct the low dimensional manifold.

4.1. Notation

Let ℱ = [F(c₁), F(c₂),…,F(c_n)] ∈ ℝ^D be a data matrix of n objects, c_i, i ∈ {1,…, n}, with dimensionality D. The corresponding class labels for these objects are given as ω_i ∈ [0,1]. Let G = {ℱ, W} be an undirected weighted graph with vertex set ℱ and similarity matrix W ∈ ℝ^n×n. W = [w_ij] assigns edge weight similarities in a pairwise fashion between objects c_i and c_j, i,j ∈ {1,…, n}. The diagonal matrix 𝒟 and Laplacian matrix L of a graph G is defined as: 𝒟_ii = Σ_jw_ij, where L = 𝒟 – W. A kernel gram matrix defining similarities between n points is given as K_m for protocol m, m ∈ {1,…,M}, where M is the total number of protocols (or data channels).

4.2. Graph embedding framework

The aim of GE (Shi and Malik, 2000) is to reduce the data matrix ℱ ∈ ℝ^D into a low-dimensional space y ∈ ℝ^d (D ≫ d), such that object adjacencies are preserved from ℝ^D to ℝ^d. GE attempts to find the optimal low dimensional vector representations among the vertices of G that best characterize the similarity relationship between the vertex pairs in G. The low dimensional representation y = [y₁, y₂,…,y_n] can be obtained by solving an Eigenvalue decomposition problem as,

with the constraint y^T𝒟y = 1. W = [w_ij] is a similarity matrix which assigns edge weights to characterize pairwise similarities between points c_i and c_j, i,j ∈ {1,…,n} such that $w_{ij}=\frac{-\parallel F(c_i)-F(c_j)\parallel }{\sigma }$, where σ is a scaling parameter.

4.3. Kernel Graph Embedding (KGE) framework

KGE is a technique to extend linear projections of data to a nonlinear dot product space using the kernel trick (Muller et al., 2001), which maps data from the original input space to an alternative higher dimensional space as K(F(c_i),F(c_j)) = 〈(Φ(F(c_i)), Φ(F(c_j)))〉, where Φ is the implicit pairwise embedding between F(c_i) and F(c_j). A kernel gram matrix K_m for each protocol m may be obtained as K_m = [K(F(c_i),F(c_j))], ∀i,j ∈ {1,…, n}, where K_m may be expressed as

According to the Representer Theorem (Schlkopf et al., 2001), to calculate the kernel representation K(F(c_i), F(c_j)) of input data, it is assumed that the optimal embedding y lies in the input space such that,

where, α = [α₁, α₂,…, α_n]^T is the low dimensional matrix representation for KGE and K(ℱ) = [K(ℱ(.), F(c₁)), K(ℱ(.),F(c₂)),…, K(ℱ(.), F(c_n))]_m, which simplifies to,

where K_m is the kernel matrix K_ij = K(F(c_i), F(c_j)). Using Eq. (4), the objective function for KGE can be reduced to,

where K is a valid positive semi-definite kernel and α is the d dimensional Eigenvector of the objective kernel function in Eq. (5). As shown in (Schlkopf et al., 2001), optimization function of KGE can again be solved by an Eigenvalue decomposition problem as,

with the constraint α^TK𝒟Kα = 1.

4.4. Semi-Supervised Multi-Kernel Graph Embedding (SeSMiK-GE)

4.5. SeSMiK-GE optimization

In Eq. (10), two variables, ᾱ and β_m (within K̂) need to be optimized simultaneously. The low dimensional representation ᾱ is optimized using a kernel ridge regression function, while the weights β_m are optimized using a hierarchical brute force algorithm. Each individual optimization step is explained below.

4.6. SeSMiK-GE algorithm

Algorithm. SeSMiK-GE

K_m and W̃_m (via ω_l, where ω_l are the exposed labels while computing W̃_m) are obtained for each of the data channels ℱ_m, m ∈ {1,…,M} in Steps 1 and 2, respectively. Combined weight matrix Ŵ is obtained by averaging the modified weight matrices (obtained in Step 2) across the M data channels (Step 4). Similarly, in Step 5, the contributions of each K_m are individually weighted using the optimal weights β̂_m across the M data channels. Ŵ and K̂ are then substituted in the generalized KGE framework in Step 6 which is solved in Step 7 to obtain the d-dimensional fused representation of the multi-modal data. A supervised classifier can then be trained on the fused low dimensional representation ᾱ for subsequent object classification. Note that when v new samples are introduced in the original data matrix ℱ of size n, the kernel matrix K_m, and the weight matrix W̃_m (with original ω_l exposed) are recomputed to obtain matrices of size n + v × n + v, and the classifiers are retrained using the recomputed low dimensional representation ᾱ for samples n + v.

5.1. Data description and preliminaries

A total of 29 1.5 Tesla (T) endorectal T2w MRI, MRS studies were obtained prior to radical prostatectomy. The 29 1.5 T studies comprised a total of 960 CaP and 1365 benign metavoxels. Of the 29 studies, 12 studies were found to have high grade CaP with 188 low and 310 high grade CaP metavoxels. All of these studies were biopsy proven prostate cancer patient studies that were clinically referred for a CaP MR staging exam for improved therapeutic selection. MRI was performed by using a 1.5-T whole-body MRI unit (Signa; GE Medical Systems, Milwaukee, Wisconsin). The patients were imaged while in the supine position by using a body coil for signal excitation and a pelvic phased-array coil (GE Medical Systems) combined with a balloon-covered expandable endorectal coil (Medrad, Pittsburgh, PA) for signal reception. Data sets were acquired as 16 × 8 × 8 phase-encoded spectral arrays (1024 voxels) by using a nominal spectral resolution of 0.24−0.34 cm³, 1000/130, and a 17-min acquisition time. Three-dimensional, MR spectroscopic imaging data were processed and aligned with the corresponding T2w imaging data using a combination of in-house software and Interactive Display Language (Research Systems, Boulder, Colorado) software tools (Tran et al., 2000). The raw spectral data were apodized with a 1-Hz Gaussian function and Fourier transformed in the time domain and in three spatial domains. Choline, creatine, and citrate peak parameters (i.e., peak area, peak height, peak location, and line width) were estimated by using an iterative procedure that was used to first identify statistically significant peaks (those with a signal-to-noise ratio higher than 5) in the magnitude spectrum. The frequency shift that best aligns the spectral peaks with the expected locations of choline, creatine, citrate, and residual water is then estimated. Subsequently, the spectra are phased by using the phase of the residual water and the metabolite resonances. Baseline values were corrected by using a local non-linear fit to the non-peak regions of the spectra. Subsequent feature extraction and classification steps were performed using algorithms developed within the MATLAB (The MathWorks, Inc.) programming environment.

The 3D prostate T2w MRI scene is represented by 𝒞̂ = (Ĉ,f̂), where Ĉ is a 3D grid of voxels ĉ ∈ Ĉ and f̂(ĉ) is a function that assigns an intensity value to every ĉ ∈ Ĉ. We also define a MR spectral scene 𝒞 = (C, F) where C is a 3D grid of MRS metavoxels, c ∈ C, and F is a spectral vector associated with each c ∈ C. Note that on account of differences in resolution of MRI/MRS, a single spectral metavoxel is several times larger compared to the size of a corresponding T2w MRI voxel. We define a CaP classifier output as h(c). The corresponding PPMM classifier is defined as h̃(c), and a high grade CaP classifier output is defined as ĥ(c). Similarly, notation for a classifier trained in conjunction with different feature vectors is identical to the corresponding notation for the feature vectors and involves replacing the F with h (e.g. a CaP classifier that leverages T2w MRI features in F^T2 (c) is denoted as h^T2(c), corresponding MRF CaP classifier as h̃^T2(c) and the corresponding high grade CaP classifier as ĥ^T2(c)). The common notations used in this paper are listed in Table 1, and the major acronyms are listed in Table 2.

5.2. Manual ground truth annotation of voxels

For all the studies considered in this work, ex vivo whole mount histological sections obtained from radical prostatectomy specimens were available (Fig. 3). An expert pathologist manually delineated CaP and grade extent on whole mount histological sections. The CaP and grade extent on the MRI were further delineated and graded by an expert prostate imaging expert (JK) (more than 25 year experience in the field of MP MRI for CaP detection and grading) by visually registering corresponding histological (Fig. 3a) and radiological sections (Fig. 3b); correspondence between sections having been determined manually by visually determining anatomical fiducials (urethra, verumontanum, large benign prostatic hyperplasia (BPH) nodules) on the histology and the imaging. Based on the CaP extent on visually registered histology, hypo intense appearance on T2w MRI, and reduced choline concentrations on MRS, MR spectral voxels were annotated as CaP, benign, or excluded (as described in Section 5.2.1). Similarly, correspondences across pathological scores from visually registered histology, hypo intense appearance on T2w MRI, and reduced choline concentrations on MRS, were used to further label the CaP MR spectral voxels as low- or high-grade (as described in Section 5.2.2). The class labels thus obtained for each of the individual MR spectral voxels were used as the surrogate ground truth for CaP detection and grading on the MRI/MRS, and were employed for subsequent training and evaluation of the SeSMiK-GE classifier.

5.3. Feature extraction from MRI and MRS

5.4. Hierarchical classification of high grade CaP via MP SeSMiK-GE signature

5.5. Implementation details and classifier training

5.6. Comparative data integration strategies

In the following sub-sections, we evaluate and compare SeSMiK with other feature extraction schemes used in the context of automated CaP detection for (i) individual MRS (Tiwari et al., 2009), T2w MRI (Madabhushi et al., 2005) modalities, (ii) a data integration (COD) scheme for combining MRS with T2w MRI features, and (iii) a decision integration strategy where the individual uni-modal classifier outputs are fused to obtain a combined classification output (Rohlfing et al., 2005). A summary of all the comparative strategies used for CaP detection and grading in this work is given in Table 4.

5.7. Experimental evaluation

6.1. Experiment 1: CaP detection via SeSMiK-GE

ROC curves were generated by thresholding CaP probabilities at different thresholds for every cross-validation run (29 CV runs for LOO and 25 runs accumulated over threefold CV runs). AUC and accuracy values were reported at the operating point of each ROC curve for both LOO and threefold CV. Fig. 6a and c show AUC results (φ^AUC), while Fig. 6b and d show accuracy results (φ^Acc) obtained across 29 iterations of LOO CV, and 25 iterations of threefold CV respectively for various feature extraction and classifier strategies (h̃^T2,h̃^MRS,h̃^Int,h̃^IntD,h̃^T2MRS) at the operating point using box-and-whiskers plots. The mean $({\phi }_{\mu }^{AUC})$ and standard deviation (ζ^AUC) of AUC values, and classifier accuracy $({\phi }_{\mu }^{Acc})$ across 29 iterations for LOO, and 25 iterations of threefold CV were also recorded and are shown in Table 5.

The results across the two CV strategies reflect that although the LOO CV strategy yielded higher φ^AUC and φ^Acc values, variance across φ^AUC and φ^Acc was significantly reduced with a threefold CV strategy. It is worth noting however that φ^AUC and φ^Acc across both CV strategies were found to be very similar, and superior to using individual modalities and other comparative strategies (h̃^T2,h̃^MRS,h̃^Int, and h̃^IntD).

Both Table 5 and Fig. 6 suggest that MP COD based data integration strategies (h̃^Int, and h̃^T2MRS) yield higher CaP detection φ^AUC and φ^Acc as compared to uni-modal h̃^T2 and h̃^MRS classifiers (Chan et al., 2003, Vos et al., 2010, Liu et al., 2009) and a COI based data integration classifier.

As suggested in (Demsar, 2006), to overcome the issue of multiple comparisons and biases with CV and parametric paired t-test, we incorporated a non-parametric Wilcoxon signed ranked test while comparing different classifier strategies across the two CV schemes. To further make the statistical significance test more stringent, the p-value was appropriately adjusted at p = 0.00125 via a Bonferroni test to account for Type 1 errors due to multiple testing. Table 6 shows the p values of statistical significance obtained using a non-parametric Wilcoxon signed test (Demsar, 2006) for comparing ${\phi }_{\mu }^{AUC}$ and ${\phi }_{\mu }^{Acc}$ for h̃^T2MRS against h̃^T2, h̃^MRS, h̃^Int, and h̃^IntD across LOO and threefold-CV strategies at d = 15. It is worth noting that the results across different d ∈ {10,15,20} were found to be consistent, where h̃^T2MRS consistently significantly outperformed (p < 0.001) the other classifiers across all reduced dimensions (results not shown). We believe that the high sensitivity and specificity of SeSMiK-GE compared to the other data integration strategies (COI and COD) is due to (a) combining and weighting individual kernel contributions within the MKL formation, and (b) employing partial label information for SSL which improves class differentiability.

Finally, it should be noted that employing PPMM classifier h̃^T2MRS $({\mathit{\text{varphi}}}_{\mu }^{Acc}=0.72\pm 0.16)$ significantly improved CaP detection accuracy (~20% improvement for d = 15) over not imposing any spatial constraints h^T2MRS$({\phi }_{\mu }^{Acc}=0.89\pm 0.07)$. The improvement was consistent across h̃^Int, h̃^IntD, h̃^MRS, and h̃^T2 over both LOO and threefold-CV (Table 5), and resonates with similar findings in (Liu et al., 2009).

6.2. Experiment 2: CaP grading via SeSMiK-GE

Fig. 7 shows qualitative results for the hierarchical classification strategy for identifying high grade CaP. Fig. 7a, d and g show CaP and high grade CaP ground truth as annotated by an expert. Yellow outline in Fig. 7a, d and g defines the CaP extent while the black outline denotes the high grade CaP extent for three different T2w MRI sections from three different patients. The corresponding probability heat maps for CaP classification are shown in Fig. 7b, e and h, where the spatial locations shown in red were identified by the respective classifiers as having a higher probability of CaP presence. Locations shown in blue were those identified as having a higher probability of being benign by the SeSMiK-GE classifier. Within the high probability CaP regions (red) in Fig. 7b, e and h, a probability of high grade CaP is further assigned, the corresponding probability heat maps for which are shown in Fig. 7c, f and i. Here, the spatial locations shown in red (Fig. 7c, f and i) were identified as having a higher probability of high grade CaP, while locations shown in blue were identified as having a higher probability of being low grade using ĥ^T2MRS. Note the high detection accuracy in accurately identifying CaP and high grade CaP using our hierarchical classification strategy across all three T2w MRI sections.

Fig. 8a shows φ^AUC while Fig. 8b shows φ^Acc values across various feature extraction and classifier strategies (ĥ^T2, ĥ^MRS, ĥ^Int, ĥ^IntD, ĥ^T2MRS) for high grade CaP using box-and-whiskers plots obtained via a LOO-CV strategy. Corresponding results for φ^AUC and φ^Acc obtained via a threefold CV strategy are shown in Fig. 8c and d respectively. The mean $({\phi }_{\mu }^{AUC})$ and standard deviation (ζ^AUC) of AUC values at the operating point of the ROC curve for LOO and threefold CV are shown in Table 7. Table 8 shows the p-values obtained by performing a paired non-parametric Wilcoxon test for comparing φ^AUC obtained from ĥ^T2MRS with other classifier strategies (ĥ^T2, ĥ^MRS, ĥ^Int, ĥ^IntD) using LOO and threefold-cross-validation strategies at d = 15.

Fig. 9 shows the average ROC curves obtained using h^T2, h^MRS, h^Int, h^T2MRS and h̃^T2MRS for CaP versus benign classification (Fig. 9a) and using ĥ^T2, ĥ^MRS, ĥ^Int, and ĥ^T2MRS for low versus high grade classification (Fig. 9b) for d = 15. Again note the improvement in AUC using the MRF based SeSMiK-GE CaP classifier (h̃^T2MRS) compared to SeSMiK-GE CaP classifier with no spatial constraints (h^T2MRS) (Fig. 9a).

We attribute the high φ^AUC and φ^Acc obtained via ĥ^T2MRS to the hierarchical classification of high grade CaP which systematically hones-in on the CaP region of interest by eliminating other confounders that might otherwise affect classification (Tiwari et al., 2009). The hierarchical scheme is especially relevant in cases where the morphologic differences between the two classes are subtle. These differences may not be appreciable in the presence of other object classes which first need to be eliminated (benign confounders) thereby allowing for the subtle differences to be accentuated. Our results are also consistent with studies that have demonstrated a higher accuracy of high grade CaP detection using MP MRI/MRS compared to individual protocols (Langer et al., 2010; Shukla-Dave et al., 2007; Shukla-Dave et al., 2009).

Fig. 10b shows the classification result obtained for a single 2D section from a 3D volume for a T2w MRI/MRS sections (Fig. 10a), the classifier results having been plotted back onto the corresponding MRS grid Fig. 10b. Fig. 10b shows the corresponding classification result obtained via SeSMiK-GE on a MRS grid, which demonstrates the high sensitivity and specificity obtained via SeSMiK-GE for detection of high grade CaP. Note that the spectra identified as low and high grade CaP (in red) by the classifier appear to be qualitatively different in terms of the relative concentrations of choline, creatine, and citrate.