Lung Cancer Pathological Image Analysis Using a Hidden Potts Model.

Nowadays, many biological data are acquired via images. In this article, we study the pathological images scanned from 205 patients with lung cancer with the goal to find out the relationship between the survival time and the spatial distribution of different types of cells, including lymphocyte, stroma, and tumor cells. Toward this goal, we model the spatial distribution of different types of cells using a modified Potts model for which the parameters represent interactions between different types of cells and estimate the parameters of the Potts model using the double Metropolis-Hastings algorithm. The double Metropolis-Hastings algorithm allows us to simulate samples approximately from a distribution with an intractable normalizing constant. Our numerical results indicate that the spatial interaction between the lymphocyte and tumor cells is significantly associated with the patient's survival time, and it can be used together with the cell count information to predict the survival of the patients.

Introduction

Lung cancer is the most common human cancer and the deadliest in the United States and globally. Non–small-cell lung cancer (NSCLC) is the most common cause of lung cancer death, accounting for up to 85% of deaths from lung cancer. Within NSCLC, adenocarcinoma and squamous cell carcinoma are the 2 major subtypes, with distinct prognoses and therapeutic remedies.1,2 Current guidelines for treating lung cancer are largely based on clinical and pathological staging systems. However, the outcome varies widely.1,2 Identifying the biomarkers that are responsible for patient risk prediction could help with treatment options, as well as understanding features of the tumor. Pathological examination of tumor tissue slides is routine in lung cancer diagnosis. The pathological images are widely available from routine clinical practices. Recent studies3–5 have shown that the growth patterns of lung tumors are associated with patient survival outcomes. The pathological image features derived from the computer-aided pathological analysis have been used to predict the survival of patients with breast cancer6,7 and complement cancer genomic profiling.7,8 These studies demonstrate the feasibility of using digital pathological image analysis for objective and unbiased clinical prognosis. However, there still lacks such an analysis for lung cancer due to the complexity and heterogeneity of the disease.

Modeling spatial correlations in images is fundamental for pathological image analysis. In statistics, Markov random field models, such as the Ising model and Potts model, have been widely used to extract spatial correlation information for imaging data,9–11 The major difficulty with the Ising and Potts models is their intractable normalizing constant, which makes their parameters hard to be estimated.

In this article, we model the spatial distribution of different types of cells, namely, lymphocyte, stroma, and tumor cells, using a modified Potts model. The parameters of the Potts model are often called interaction parameters, which characterize the clustering behavior of the same types of spins (ie, cells in the context of this article). We estimate the parameters of the Potts model using the double Metropolis-Hastings (DMH) algorithm12 under a Bayesian framework. The DMH algorithm is very efficient for sampling from distributions with intractable normalizing constants, especially for the distributions defined on a large-scale lattice. We found that the interaction between lymphocytes and tumor cells is significantly associated with the patient’s survival time, and furthermore, it can be used together with the cell count information to improve the prediction of the patient’s survival time.

The remainder of this article is organized as follows. Section “Potts model” describes the modified Potts model and gives the details on how the DMH algorithm can be used to estimate its parameters. Section “Lung cancer imaging data” proposes a hidden, modified Potts model and presents our findings for lung cancer pathological image data. Section “Discussion” concludes the article with a brief discussion.

Potts Model

A modified Potts model

The q-state Potts model13 consists of a 2-dimensional lattice of spins, where each spin takes values from a set of q different elements. The energy function of the model is given by where = {x} denotes the collection of all spins of the model, denotes the neighboring pairs of the spins, J,′′ denotes the interaction parameter between the spin x and the spin x′′, δ(x, x′′) is an indicator function which is equal to 1 if x = x′′ and 0 otherwise, and the sum takes over all neighboring pairs of spins. If q = 2, the Potts model is reduced to the Ising model.14 The Potts model is also related to, and generalized by, several other models, such as the XY model,15 the Heisenberg model16 and the N-vector model.17 Generalizations of the Potts model have been used to model grain growth in metals and coarsening in foams.18,19 A further generalization of the model, known as the cellular Potts model,20 has been used to simulate static and kinetic phenomena in foam and biological morphogenesis.

In the standard form of Potts model (1), J,′′ represents the strength of interaction between the same type of neighboring spins as the function δ(x, x′′) takes a non-zero value if and only if x = x′′, where the value of x indicates the type of the spin. However, for the lung cancer pathological imaging data, we would like to study the interactions between different types of cells, which are coded as 1 for lymphocyte cells, 2 for stroma cells, and 3 for tumor cells. For this reason, we consider a modified Potts model with the energy function given by where k and l take values from the set {1,2,3}, θ12 represents the interaction between lymphocyte and stroma cells, θ13 represents the interaction between lymphocyte and tumor cells, and θ23 represents the interaction between stroma and tumor cells. This modified Potts model can also be viewed as a simplified cellular Potts model without the volume and surface constraints.

Corresponding to the energy function (2), the probability mass function of the modified Potts model is given by where = { θ12, θ13, θ23}, and Z() is the normalizing constant. As an exact evaluation of Z() needs to sum over the entire space of , which consists of 3 different elements with N denoting the total number of spins, Z() is intractable even for a small size model, say N = 100. How to estimate the parameters for such a model has been studied in the recent literature.9,12,21–23

DMH algorithm for the modified Potts model

Suppose that we are interested in estimating for the modified Potts model under a Bayesian framework. Let π() denote the prior density function of . Then, the posterior density function is given by

It is easy to see that the Metropolis-Hastings (MH) algorithm cannot be directly applied to simulate from this posterior as the acceptance probability would involve an unknown normalizing constant ratio Z()/Z(′), where ′ denotes the proposed value.

To address this issue, some auxiliary variable Markov chain Monte Carlo (MCMC) algorithms have been proposed, which aim to have the normalizing constant ratio Z()/Z(′) canceled in simulations by augmenting appropriate auxiliary variables to the target distribution and/or the proposal distribution. Along this direction, Møller et al21 proposed an algorithm which augments both the target and proposal distributions, and Murray et al22 proposed the so-called exchange algorithm which arguments only the proposal distribution. Although the underlying idea is very attractive, these 2 algorithms require the auxiliary variables to be drawn using a perfect sampler.24 As perfect sampling can be very expensive or impossible for many models with intractable normalizing constants, the applications of these algorithms are highly hindered. To overcome this difficulty, Liang12 proposed the DMH algorithm, and Liang et al23 proposed an adaptive exchange algorithm. The adaptive exchange algorithm generates auxiliary variables via an importance sampling procedure from a Markov chain running in parallel, and the DMH algorithm generates auxiliary variables through a short run of the MH algorithm initialized with the original observation. As noted in Liang,12 initializing the auxiliary MH chain with the original observation improves convergence of the algorithm. Other than the auxiliary variable MCMC algorithm, some approximation-based algorithms have been proposed in the literature, such as maximum pseudo-likelihood estimation,25 Monte Carlo maximum likelihood estimation,26 and adaptive kernel smoothing,27,28 which approximate the likelihood function, the normalizing constant Z(), or the normalizing constant ratio Z()/Z(′). A recent comparative review29 concludes that, compared with other algorithms, the DMH algorithm is very efficient for complex models with intractable normalizing constants, although the estimates are only asymptotically correct. The DMH algorithm is adopted in this article for estimating the parameters of the modified Potts model. The DMH algorithm can be described as follows.

Suppose that we are interested in simulating a sample from f(·|′) using the MH algorithm. If starting with the current state x, the transition probability, , is given by where K (·→ ·) is the MH transition kernel. Provided that the Markov chain has reached equilibrium states, then, by the detailed balance condition, we have

Let q(′|, ) denote the proposal distribution for drawing a new parameter vector ′. The DMH algorithm iterates between the following steps:

Draw a new sample ′ from the proposal density function q(′|, ).

Given , simulate an auxiliary variable starting from the observation and calculate the MH ratio:

Set +1 = ′ with probability min{1, r(, ′, |)}; set +1 = otherwise.

Suppose that a sequence of samples 1,…, has been collected from a run of DMH. An approximate Bayesian estimator of can then be obtained by averaging over the samples:

As implied by equation (6), the DMH sampler is almost exact for those parameters around the true value of . Hence, the estimator can be rather accurate even when m is not very large.

Simulation examples

We tested the performance of the DMH algorithm on the modified Potts model using simulated examples. We considered a variety of values of as given in Table 1. For each setting of , we simulated 30 data sets independently using the Gibbs sampler on a 50 × 50 lattice. To simulate each data set, the Gibbs sampler was run for 105 iterations with random starting configurations.

Table 1

Parameter estimates for the simulated data by the double Metropolis-Hastings algorithm, where SE(·) denotes the standard error of the corresponding estimate.

To conduct a Bayesian analysis for the simulated data, we let be subject to a uniform prior distribution, ie, π() ∝ 1 for ∈[0,1]3. For each of the simulated data sets, DMH was run for 6000 iterations, where the rst 1000 iterations were discarded for the burn-in process and the samples generated in the remaining iterations were used for inference. At each iteration, the auxiliary sample was simulated using the Gibbs sampler with a single sweep for all elements. Each run costs about 18 seconds central processing unit (CPU) time on a Dell OptiPlex 9020 computer. The estimates of are summarized in Table 1, where each estimate is obtained by averaging more than 30 independent data sets. Table 1 indicates that the DMH algorithm works well in parameter estimation for the modified Potts model.

Lung Cancer Imaging Data

The data set and accessibility

The pathological images and survival status from 205 patients with NSCLC in the National Lung Screening Trial (NLST) study30 were collected. The characteristics of the participants are summarized in Table 2. The Kaplan-Meier curve for survival of the whole set of patients is shown in Figure 1. For each patient, 1, 2, or 3 tissue slides were first taken, where the number of slides depends on the size of tumor. For patients with a large size of tumor, more slides are needed to have a more comprehensive characterization of the tumor, and vice versa. Then, each tissue slide was examined by a lung cancer pathologist, the regions of interest (ROIs) within the tumor region(s) were determined, and 5 ROIs were randomly selected from each patient for further analysis. In total, we had 1585 ROI images. In the ROIs, the nucleus of each cell and the cell boundary were determined using a watershed method.31 For each cell, a 160 pixels by 160 pixels image patch was extracted around the center, and the cell type was predicted using a convolutional neural network32 developed from another study. The prediction was evaluated by the lung cancer pathologist, and the accuracy was more than 94%. The cell locations (ie, the coordinates of the cell on the ROI image) and the predicted cell types were used as inputs of the proposed method in this article. We are making the code publicly available. Once the paper is accepted, the code will be linked to the published version of the article.

Table 2

NSCLC patient characteristics (N = 205).

Figure 1

Kaplan-Meier plot for 205 patients with non–small-cell lung cancer.

DMH for a hidden Potts model

As the cell locations are irregular, it is difficult to model the pathological image by a Potts model. In particular, it is difficult to identify the neighboring cells for each cell. For this reason, we model the image by a hidden Potts model by introducing an auxiliary lattice to the image, which is illustrated by Figure 2. We note that the idea of modeling spatial data via an auxiliary lattice has been explored in the literature.33,34

Figure 2

Illustration of the auxiliary lattice, where the circles represent cells. In real data sets, the location of a cell is given by a point so that a cell cannot belong to more than 1 square.

Consider a pathological image with n cells located at s1,…,s. Let y denote the type of the cell located at s, and it takes value 1 for lymphocyte cells, 2 for stroma cells, and 3 for tumor cells. Let = {y} denote the collection of observed cell types in the image. Let W = {(i, j) : i = 1,…M, j = 1,…,N} denote the auxiliary lattice, which partitions the image into (M + 1)(N + 1) squares. Denote the squares by C1,C2,…,C( + 1)( + 1). Let C denote the square that s belongs to. For convenience, we let each C be compact, ie, including all the boundary points of the square, while assuming that there are no cells belonging to 2 squares. If a cell is exactly on the boundary of some squares, then we randomly assign the cell to one of them. Let X,(i, j) ϵ W, denote the hidden cell types at the auxiliary lattice. Conditional on X, we model the distribution of {Y} by where γ is a projection parameter with a prespecified value. The larger γ is, the more similar the original and imputed images are.

Let f ( | ), as specified in equation (3), denote the likelihood function of the hidden Potts model. Let π() denote the prior density function of . Assume that all are independent conditional on , then we have

Therefore, the full conditional posterior of X is given by where qϵ {1,2,3}, C denotes the union of the squares that the spin (i, j) belongs to, and ∂ denotes the neighboring spins of the spin (i, j). In this article, a free boundary condition is assumed for the model, under which the boundary points have fewer neighboring spins than the interior spins. After the normalization for equation (9), we have

Therefore, given , the projection parameter, and the model parameter, the hidden Potts model can be imputed using the Gibbs sampler through iteratively drawing X s from distribution (10).

Similarly, we can get the full conditional posterior of :

As γ is a prespecified constant, the posterior can be simulated by iterating between the following 2 steps:

DMH algorithm for hidden Potts models

Impute the hidden Potts model by simulating from equation (10) for all (i, j) ∈W.

Simulate from distribution (11) using the DMH algorithm.

This algorithm consists of a few tunable parameters, including γ, M, and N. As mentioned previously, determine the similarity of the observed and imputed images. To make the 2 images more similar, we set γ to a large value. In all examples of this article, we set γ = 10. The parameters M and N determine the size of the auxiliary lattice. Following the suggestion of Park and Liang,33 we choose M and N such that n ≈ MN. To be precise, we set M and N in the following way for each image: Let d1 and d2 denote the ranges of the cell locations at x-axis and y-axis, respectively. Then, we set the side length of each square to , and set M = [d1 / l ]+ 1 and N = [d2 / l ]+ 1, where [z] denotes the integer part of z.

Numerical results

The algorithm described above was applied to the 1585 ROI images. For each image, the algorithm was run for 6000 iterations, where the first 1000 iterations were discarded for the burn-in process and the samples generated in the remaining iterations were used for inference, and it cost about 14 minutes of CPU time on a Dell OptiPlex 9020 computer. The CPU time may vary slightly according to the values of M and N. Figure 3 shows the observed (left panels) and imputed (right panels) images for 3 ROIs, where each of the imputed images is obtained by averaging over the samples generated in a single run of the DMH algorithm. As it takes a large value, the imputed image is almost the same at each iteration after the simulation has reached equilibrium. As shown in Figure 3, the imputed images are very similar to the observed ones.

Figure 3

Comparison of the observed (left panels) and imputed (right panels) images for 3 regions of interest which are all from different patients: lymphocyte cells are in blue, stroma cells are in red, and tumor cells are in light gray.

To assess the association between the spatial distributions of cells in pathological images and patients’ survival status, we fitted a Cox proportional hazards model on the survival time with respect to the estimates of the interaction parameters of the hidden Potts model. Here, the survival time is defined as the time from diagnosis (of lung cancer) to death from all causes; right-censored cases exist. A patient with censored survival time means the patient is still alive at the last follow-up or lost to follow-up. In the Cox regression model, the hazard function has the form where h(t) is the baseline hazard function, and Z1,…,Z are covariates having a multiplicative effect on the hazard function. For our model, we have k = 3 and the covariates , , and , which were produced by the proposed method with the projection parameter γ = 10. The parameters β1,…,β are estimated using the R package “survival.”35 As the data set contains multiple observations per patient, the generalized estimating equation method was used to compute a robust variance for each parameter estimate.36

Table 3 summarizes the estimates of the parameters of the Cox regression model. The overall P value for the significance of the model is .03374. The parameter estimates indicate that a higher value of the interaction between lymphocytes and tumor cells is significantly associated with a higher risk of death (at a significance level of .05); the hazard coefficient shows a negative correlation between the survival time and the value of the interaction between lymphocyte and tumor cells. In other words, widespread tumor cells indicate severity of the disease. Table 3 also shows that the interaction between stroma and tumor cells is also weekly associated with the risk of death (at a significance level of .1). However, there is no evidence suggesting any association between the interaction of lymphocytes and stroma cells and the risk of death. In summary, we may conclude that the proposed hidden Potts model is able to extract some useful information about the status of the disease.

Table 3

Survival analysis for lung cancer pathological images with the cell spatial interaction information (γ = 10).

To assess the sensitivity of the results to the projection parameter γ, different values of γ were tried. The results are reported in Tables 4 and 5, which indicate that our results are not sensitive to the value of γ.

Table 4

Survival analysis for lung cancer pathological images with the cell spatial interaction information (γ = 8).

Table 5

Survival analysis for lung cancer pathological images with the cell spatial interaction information (γ = 12).

We also assessed the proportional hazards assumption for the Cox regression model on this data set.37 The P value for the whole model is .219, and the P values of the respective parameters are all greater than .10, suggesting that the proportional hazards assumption is valid and the regression coefficients = {θ12, θ13, θ23} remain constant over time. Figure 4 shows the plot of scaled Schoenfeld residuals versus time.

Figure 4

Scaled Schoenfeld residuals versus time for ={θ12,θ13,θ23}.

Finally, we conducted a survival analysis based on the cell count information only. In particular, we considered the cell count ratios, lymphocyte/stroma and tumor/stroma. The results, which are shown in Table 6, indicate that the cell count information is also very useful in predicting patient survival; the hazard coefficient implies a negative correlation between the survival time and the ratio tumor/stroma. It is very interesting to point out that the spatial interaction information learned by the proposed method for different types of cells is complementary to the cell count information. As indicated in Table 7, the prediction for the survival can be further improved by using both of them. With both information, the overall P value (in likelihood ratio test) of the Cox regression model has been improved to .00539 from .02367 (with the cell count information only). In addition, compared with Tables 3 and 6, the significance levels of θ13, θ23, and tumor/stroma ratio are also improved.

Table 6

Survival analysis for lung cancer pathological images with cell count information only.

Table 7

Survival analysis for lung cancer pathological images with both the cell count information and the cell spatial interaction information.

Discussion

In this article, we have proposed to model pathological images using a hidden Potts model and applied the DMH algorithm to estimate the model parameters. The introduction of auxiliary lattice makes the proposed method very general, which can be used for any type of imaging data with or without regular observations. The auxiliary lattice also helps reduce the complexity of imaging data and defines a concise and explicit neighborhood for each spin of the hidden Potts model. Other auxiliary variable MCMC algorithms, eg, the adaptive exchange algorithm,23 can potentially be applied to this problem. However, it would be more time-consuming given the hidden structure of the proposed Potts model.

For the lung cancer pathological imaging data, our study shows that the survival time of NSCLC patients might be significantly associated with the strength of interactions between lymphocyte and tumor cells. The spatial interaction parameter together with the cell count information can potentially be used as a biomarker for prognosis and personalized treatments of patients with NSCLC. It would be of great interest to extend the proposed method to other pathological imaging data.