Keywords

These keywords were added by machine and not by the authors. This process is experimental and the keywords may be updated as the learning algorithm improves.

1 Introduction

Currently screening of prostate cancers is usually performed through routine prostate-specific antigen (PSA) blood tests and/or a rectal examination. Based on positive PSA indication, a biopsy of randomly sampled areas of the prostate can then be considered to diagnose the cancer and assess its aggressiveness. Biopsy may miss sampling cancerous tissues, resulting in missed or delayed diagnosis, and miss areas with aggressive cancers, thus under-staging the cancer and leading to under-treatment.

Studies have shown that the tissue stiffness described by the tissue properties may indicate abnormal pathological process. Ex-vivo, measurement-based methods, such as [1, 11] using magnetic resonance imaging (MRI) and/or ultrasound, were proposed for study of prostate cancer tissue. However, previous works in material property reconstruction often have limitations with respect to their genericity, applicability, efficiency and accuracy [22]. More recent techniques, such as inverse finite-element methods [6, 13, 17, 21, 22], stochastic finite-element methods [18], and image-based ultrasound [20] have been developed for in-vivo soft tissue analysis.

In this paper, we study the possible use of tissue (i.e. prostate) elasticity to help evaluate the prognosis of prostate cancer patients given at least two set of CT images. The clinical T-stage of a prostate cancer is a measure of how much the tumor has grown and spread; while a Gleason score based on the biopsy of cancer cells indicates aggressiveness of the cancer. They are commonly used for cancer staging and grading. We present an improved method that uses geometric and physical constraints to deduce the relative tissue elasticity parameters. Although elasticity reconstruction, or elastography, can be used to estimate tissue elasticity, it is less suited for in-vivo measurements or deeply seated organs like prostate. We describe a non-invasive method to estimate tissue elasticity values based on pairs of CT images, using a finite-element based biomechanical model derived from an initial set of images, local displacements, and an optimization-based framework.

Given the recovered tissue properties reconstructed from analysis of medical images and patient’s ages, we develop a multiclass classification system for classifying clinical T-stage and Gleason scores for prostate cancer patients. We demonstrate the feasibility of a statistically-based multiclass classifier that classifies a supplementary assessment on cancer T-stages and cancer grades using the computed elasticity values from medical images, as an additional clinical aids for the physicians and patients to make more informed decision (e.g. more strategic biopsy locations, less/more aggressive treatment, etc.). Concurrently, extracted image features [8–10] using dynamic contrast enhanced (DCE) MRI have also been suggested for prostate cancer detection. These methods are complementary to ours and can be used in conjunction with ours as a multimodal classification method to further improve the overall classification accuracy.

2 Method

Our iterative simulation-optimization-identification framework consists of two alternating phases: the forward simulation to estimate the tissue deformation and inverse process that refines the tissue elasticity parameters to minimize the error in a given objective function. The input to our framework are two sets of 3D images. After iterations of the forward and inverse processes, we obtain the best set of elasticity parameters. Below we provide a brief overview of the key steps in this framework and we refer the interested readers to the supplementary document at http://gamma.cs.unc.edu/CancerClass/ for the detailed mathematical formulations and algorithmic process to extract the tissue elasticity parameters from medical images.

2.1 Forward Simulation: BioTissue Modeling

In our system, we apply Finite Element Method (FEM) and adopt Mooney Rivlin material for bio-tissue modeling [3]. After discretization using FEM, we arrive at a linear system,

$$\begin{aligned} \mathbf {K}\mathbf {u}=\mathbf {f} \end{aligned}$$
(1)

with \(\mathbf {K}\) as the stiffness matrix, \(\mathbf {u}\) as the displacement field and \(\mathbf {f}\) as the external forces. The stiffness matrix \(\mathbf {K}\) is not always symmetric possitive definite due to complicated boundary condition. The boundary condition we applied is the traction forces (shown in Fig. 7(a) of the supplementary document) computed based on the displacement of the surrounding tissue (overlapping surfaces shown in Fig. 7(b) of the supplementary document). We choose to use the Generalized Minimal Residual (GMRES) [16] solver to solve the linear system instead of the Generalized Conjugate Gradient (GCG) [14], as GMRES can better cope with non-symmetric, positive-definite linear system.

The computation of the siffness matrix \(\mathbf {K}\) in Eq. 1 depends on the energy function \(\varvec{\varPsi }\) of the Mooney Rivlin material model [15, 19].

$$\begin{aligned} \varvec{\varPsi }=\frac{1}{2}\mu _1((\mathbf {I}_{1}^{2}-\mathbf {I}_{2})/\mathbf {I}_{3}^{\frac{2}{3}}-6)+\mu _2(\mathbf {I}_1/\mathbf {I}_3^{\frac{1}{3}}-3)+v_1(\mathbf {I}_3^{\frac{1}{2}}-1)^2, \end{aligned}$$
(2)

where \(\mu _1\), \(\mu _2\) and \(v_1\) are the material parameters. In this paper, we recover parameters \(\mu _1\) and \(\mu _2\). Since prostate soft tissue (without tumors) tend to be homogenous, we use the average \(\bar{\mu }\) of \(\mu _1\) and \(\mu _2\) as our recovered elasticity parameter. To model incompressibility, we set \(v_1\) to be a very large value (\(1+e7\) was used in our implementation). \(v_1\) is linearly related to the bulk modulus. The larger the bulk modulus, the more incompressible the object.

Relative Elasticity Value: In addition, we divide the recovered absolute elasticity parameter \(\bar{\mu }\) by the that of the surrounding tissue to compute the relative elasticity parameter \(\hat{\mu }\). This individualized relativity value helps to remove the variation in mechanical properties of tissues between patients, normalizing the per-patient fluctuation in absolute elasticity values due to varying degrees of hydration and other temporary factors. We refer readers to our supplementary document for details regarding non-linear material models.

2.2 Inverse Process: Optimization for Parameter Identification

To estimate the patient-specific relative elasticity, our framework minimizes the error due to approximated parameters in an objective function. Our objective function as defined in Eq. 3 consists of the two components. The first part is the difference between the two surfaces – one reconstructed from the reference (initial) set of images, deformed using FEM simulation with the estimated parameters toward the target surface, and one target surface reconstructed from the second set of images. This difference is measured by the Hausdorff distance [4]. In addition we add a Tikhonov regularization [5, 7] term, which improves the conditioning of a possibly ill-posed problem.

With regularization, our objective function is given as:

$$\begin{aligned} \varvec{\mu }=\underset{{\varvec{\mu }}}{{\text {argmin}}}\sum \Vert \mathbf {d}(\mathbf {S}_l,\mathbf {S}_t)\Vert ^2+\lambda \varvec{\varGamma }\mathbf {S}_l, \end{aligned}$$
(3)

with \(\mathbf {d}(\mathbf {S}_l,\mathbf {S}_t)\) as the distance between deformed surface and the reference surface, \(\lambda \) as the regularization weight, and \(\varvec{\varGamma }\) as the second-order differentiatial operator.

The second-order differential operator \(\varvec{\varGamma }\) on a continuous surface (2-manifolds) \(\mathbf {S}\) is the curvatures of a point on the surface. The curvature is defined through the tangent plane passing that point. We denote the normal vector of the tangent plane as \(\mathbf {n}\) and the unit direction in the tangent plane as \(\mathbf {e}_\theta \). The curvature related to the unit direction \(\mathbf {e}_\theta \) is \(\kappa (\theta )\). The mean curvature \(\kappa _{mean}\) for a continuous surface is defined as the average curvature of all the directions, \(\kappa _{mean}=\frac{1}{2\pi }\int ^{2\pi }_0\kappa (\theta )\mathrm {d}\theta \). In our implementation, we use triangle mesh to approximate a continuous surface. We use the 1-ring neighbor as the region for computing the mean curvature normal on our discrete surface \(\mathbf {S}_l\). We treat each triangle of the mesh as a local surface with two conformal space parameters u and v. With these two parameters u and v the second-order differential operator \(\varvec{\varGamma }\) on vertex \(\mathbf {x}\) is, \(\varDelta _{u,v}\mathbf {x}=\mathbf {x}_{uu}+\mathbf {x}_{vv}\).

2.3 Classification Methods

For classification of cancer prognostic scores, we develop a learning method to classify patient cancer T-Stage and Gleason score based on the relative elasticity parameters recovered from CT images. Both the prostate cancer T-stage and the Gleason score are generally considered as ordinal responses. We study the effectiveness of ordianl logistic regression [2] and multinomial logistic regression [12] in the context of prostate cancer staging and grading. For both cases we use RBF kernel to project our feature to higher dimentional space. We refer readers to supplementary document for method details and the comparison with the Random Forests method.

Fig. 1.
figure 1

Real Patient CT Image and Reconstructed Organ Surfaces. (a) shows one slice of the parient CT images with the bladder, prostate and rectum segmented. (b) shows the reconstructed organ surfaces.

Full size image

3 Patient Data Study

3.1 Preprocessing and Patient Dataset

Given the CT images (shown in Fig. 1a) of the patient, the prostate, bladder and rectum are first segmented in the images. Then the 3D surfaces (shown in Fig. 1b) of these organs are reconstructed using VTK and these surfaces would be the input to our elasticity parameter reconstruction algorithm. Our patient dataset contains 113 (29 as the reference and 84 as target) sets of CT images from 29 patients, each patient having 2 to 15 sets of CT images. Every patient in the dataset has prostate cancer with cancer T-stage ranging from T1 to T3, Gleason score ranging from 6 to 10, and age from 50 to 85. Gleanson scores are usually used to assess the aggressiveness of the cancer.

3.2 Cancer Grading/Staging Classification Based on Prostate Elasticity Parameters

We further study the feasibility of using recovered elasticity parameters as a cancer prognostic indicator using our classifier based on relative tissue elasticity values and ages. Two classification methods, ordinal logistic regression and multinomial logistic regression, were tested in our study. We test each method with two sets of features. The first set of features contains only the relative tissue elasticity values \(\hat{\mu }\). The resultant feature vector is one dimension. The second set of features contains both the relative tissue elasticity values and the age. The feature vector for this set of features is two dimensional. Our cancer staging has \(C=3\) classes, T1, T2 and T3. And the cancer grading has \(G=5\) classes, from 6 to 10. In our patient dataset, each patient has at least 2 sets of CT images. The elasticity parameter reconstruction algorithm needs 2 sets of CT images as input. We fix one set of CT images as the initial (reference) image and use the other M number of images \(\mathcal {T}, where\ |\mathcal {T}|=M\) as the target (deformed) images. By registering the initial image to the target images, we obtain one elasticity parameter \(\hat{\mu }_i,i=1\dots M\) for each image in \(\mathcal {T}\). We perform both per-patient and per-image cross validation.

Per-Image Cross Validation: We treat all the target images (\(N=84\)) of all the patients as data points of equal importance. The elasticity feature for each target image is the recovered elasticity parameter \(\hat{\mu }\). In this experiment, we train our classifier using the elasticity feature of the 83 images then cross validate with the one left out. Then, we add the patient’s age as another feature to the classifier and perform the validation. The results for cancer staging (T-Stage) classification are shown in Fig. 2a and that for cancer grading (Gleason score) classification are shown in Fig. 2b. The error metric is measured as the absolute difference between the classified cancer T-Stage and the actual cancer T-Stage. Zero error-distance means our classifier accurately classifies the cancer T-Stage.

The multinomial method outperforms the ordinal method for both cancer staging (T-Stage) and cancer aggression (Gleason score) classification. The main reason that we are observing this is due to the optimization weights or the unknown regression coefficients \(\varvec{\beta }\) (refer to supplementary document for the definition) dimension of the multinomial and ordinal logistic regression method. The dimension of the unknown regression coefficients of the multinomial logistic regression for cancer staging classification (with elasticity parameter and age as features) is 6 while that of ordinal logistic regression is 4. With the ‘age’ feature, we obtain up to \(\mathbf {91\,\%}\) accuracy for perdicting cancer T-Stage using multinomial logistic regression method and \(\mathbf {89\,\%}\) using ordinal logistic regression method. For Gleason score classification we achieve up to \(\mathbf {88\,\%}\) accuracy using multinomial logistic regression method and \(\mathbf {81\,\%}\) using ordinal logistic regression method.

Fig. 2.
figure 2

Error Distribution of Cancer Grading/Staging Classification for Per-Image Study. (a) shows error distribution of our cancer staging classification using the recovered prostate elasticity parameter and the patient’s age. For our patient dataset, the multinomial classifier (shown in royal blue and sky blue) outperforms the ordinal classifier (shown in crimson and coral). We achieve up to \(\mathbf {91\,\%}\) accuracy using multinomial logistic regression and \(\mathbf {89\,\%}\) using ordinal logistic regression for classifying cancer T-Stage based on recovered elasticity parameter and age. (b) shows the correlation between the recovered relative elasticity parameter and the Gleason score with/without the patient’s age. We achieve up to \(\mathbf {88\,\%}\) accuracy using multinomial logistic regression and \(\mathbf {81\,\%}\) using ordinal logistic regression for classifying Gleason score based on recovered elasticity parameter and age.

Full size image

Per-Patient Cross Validation: For patients with more than 2 sets of images, we apply Gaussian sampling to \(\hat{\mu }_i,i=1\dots M\) to compute the sampled elasticity parameter as the elasticity feature of the patient. We first train our classifier using the elasticity feature of the 28 patients then test the trained classifier on the remaining one patient not in the training set. We repeat this process for each of the 29 patients. Then we include the patient age as another feature in the classifier. The error distribution for cancer staging (T-Stage) classification results are shown in Fig. 3a and the error distribution of cancer grading (Gleason score) classification are shown in Fig. 3b. We observe that the multinomial method in general outperforms the ordinal method. More interestingly, the age feature helps to increase the classification accuracy by \(2\,\%\) for staging classification and \(7\,\%\) for Gleason scoring classification). With the age feature, our multinomial classifier achieves up to \(\mathbf {84\,\%}\) accuracy for classifying cancer T-Stage and up to \(\mathbf {77\,\%}\) accuracy for classifying Gleason scores. And our ordinal classifier achieves up to \(\mathbf {82\,\%}\) for cancer T-Stage classification and \(\mathbf {70\,\%}\) for Gleason score classification. The drop in accuracy for per-patient experiments compared with per-image ones is primary due to the decrease in data samples.

Among the \(16\,\%\) failure cases for cancer staging classification, \(15\,\%\) of our multinomial classification results with age feature is only 1 stage away from the ground truth. And for the failure cases for scoring classification, only \(10\,\%\) of the classified Gleason scores is 1 away from the ground truth and \(13\,\%\) of them are 2 away from the ground truth.

Fig. 3.
figure 3

Error Distribution of Cancer Aggression/Staging Classification for Per-Patient Study. (a) shows the accuracy and error distribution of our recovered prostate elasticity parameter and cancer T-Stage. For our patient dataset, the multinomial classifier (shown in royal blue and sky blue) outperforms the ordinal classifier (shown in crimson and coral). We achieve up to \(\mathbf {84\,\%}\) accuracy using multinomial logistic regression and \(\mathbf {82\,\%}\) using ordinal logistic regression for classifying cancer T-Stage based on our recovered elasticity parameter and patient age information. (b) shows the correlation between the recovered relative elasticity parameter and the Gleason score. We achieve up to \(\mathbf {77\,\%}\) accuracy using multinomial logistic regression and \(\mathbf {70\,\%}\) using ordinal logistic regression for classifying Gleason score based on our recovered elasticity parameter and patient age information.

Full size image

4 Conclusion and Future Work

In this paper, we present an improved, non-invasive tissue elasticity parameter reconstruction framework using CT images. We further studied the correlation of the recovered relative elasticity parameters with prostate cancer T-Stage and Gleason score for multiclass classification of cancer T-stages and grades. The classification accuracy on our patient dataset using multinormial logistic regression method is up to \(84\,\%\) accurate for cancer T-stages and up to \(77\,\%\) accurate for Gleason scores. This study further demonstrates the effectiveness of our algorithm for recovering (relative) tissue elasticity parameter in-vivo and its promising potential for correct classification in cancer screening and diagnosis.

Future Work: This study is performed on 113 sets of images from 29 prostate cancer patients all treated in the same hospital. More image data from more patients across multiple institutions can provide a much richer set of training data, thus further improving the classification results and testing/validating its classification power for cancer diagnosis. With more data, we could also apply our learned model for cancer stage/score prediction. And other features, such as the volume of the prostate can also be included in the larger study. Another possible direction is to perform the same study on normal subjects and increase the patient diversity from different locations. A large-scale study can enable more complete analysis and lead to more insights on the impact of variability due to demographics and hospital practice on the study results. Similar analysis and derivation could also be performed using other image modalities, such as MR and ultrasound, and shown to be applicable to other types of cancers.