A Bayesian Nonparametric Approach for Functional Data Classification with Application to Hepatic Tissue Characterization.

Computed tomography perfusion (CTp) is an emerging functional imaging technology that provides a quantitative assessment of the passage of fluid through blood vessels. Tissue perfusion plays a critical role in oncology due to the proliferation of networks of new blood vessels typical of cancer angiogenesis, which triggers modifications to the vasculature of the surrounding host tissue. In this article, we consider a Bayesian semiparametric model for the analysis of functional data. This method is applied to a study of four interdependent hepatic perfusion CT characteristics that were acquired under the administration of contrast using a sequence of repeated scans over a period of 590 seconds. More specifically, our modeling framework facilitates borrowing of information across patients and tissues. Additionally, the approach enables flexible estimation of temporal correlation structures exhibited by mappings of the correlated perfusion biomarkers and thus accounts for the heteroskedasticity typically observed in those measurements, by incorporating change-points in the covariance estimation. This method is applied to measurements obtained from regions of liver surrounding malignant and benign tissues, for each perfusion biomarker. We demonstrate how to cluster the liver regions on the basis of their CTp profiles, which can be used in a prediction context to classify regions of interest provided by future patients, and thereby assist in discriminating malignant from healthy tissue regions in diagnostic settings.

Introduction

Functional imaging techniques have been progressively used in oncology to improve tumor detection, characterization, staging, prognosis, and treatment of cancer patients.1 In particular, computed tomography perfusion (CTp) has recently emerged as a widely available, inexpensive and noninvasive, imaging technique to evaluate changes in tumor vascular physiology and tumor biology. CTp can potentially lead to biomarkers for tumor characterization, as well as for early assessment of therapeutic response in multiple cancers.2,3 Any modern CT scanner system can perform CTp: a small bolus of iodine-based intravenous contrast is injected rapidly over a target region, and repeated images are taken at regular frequent time points before, during, and after the passage of the contrast agent through the tumor vasculature. Hence, this dynamic acquisition allows to measure temporal changes in tissue density after injection with the contrast agent. Tissue perfusion can then be estimated in the target region through the use of different kinetic models.4

The importance of assessing tumor perfusion and vascular permeability is related to their association with tumor angiogenesis, ie, the formation of new blood vessels, in several tumors.5–8 Neovascularization is considered an important process in cancer progression and tumor growth. For example, high tumor angiogenesis activity has been often associated with distant metastases and adverse clinical outcomes.9–11

The angiogenic vasculature of malignant solid tumors is usually characterized by dense, dilated, and tortuous microvessels. In particular, blood volume (BV) and blood flow (BF) are higher in malignant tumors due to the proliferation of new blood vessels through the process of tumor angiogenesis.12,13 Therefore, it has been suggested that perfusion imaging can yield biomarkers of angiogenesis and tumor growth and thereby greatly enhance the clinical development of antiangiogenic therapies. This is because changes in tumor perfusion appear to occur soon after therapy initiation with antiangiogenic drugs.14–16

CTp parameters are commonly calculated after the scan session by using commercially available postprocessing software platforms (CT Perfusion 4, version 4.3.1, Advantage Windows 4.4; GE Healthcare). Typical CTp parameters include BF (mL/100 g/minute), BV (mL/100 g), mean transit time (MTT, seconds), and permeability surface area product (PS, mL/100 g/minute), which are obtained using standard deconvolution physiological modeling, based on the central volume principle (BF = BV/MTT) first described in the context of cerebral perfusion.17,18 PS is a measure of capillary permeability and it reflects the flux of solutes from blood plasma to the interstitial space.19 As an illustration, Figure 1 reports the time course of the PS parameter observed during the imaging session for two representative normal and tumor regions of interest (ROIs) in a single patient. The disruption and decline in PS values within the first 100 seconds of the scan in Figure 1B are usually the result of neoangiogenesis and microvascular attenuation.

Figure 1

Illustrative plot of the observed measurements on the PS perfusion characteristic in two ROIs for a representative patient, as a function of scan time, zoomed in the first 100 seconds of the scan: (A) normal tissue and (B) tumor tissue.

The trajectory of the perfusion parameters over the scan period can be affected by many factors, which may pertain to patient and tumor characteristics. For example, in squamous cell carcinoma, BF and PS values have been shown to be significantly higher in subjects with longer local tumor control than in subjects with local recurrence.20 This consideration highlights the heterogeneity of the perfusion parameter–time curves both between ROIs within a patient as well as between patients, in addition to the fact that the temporal resolution and duration of the CT acquisition can also impact the CT perfusion parameter values.21,22 Yet, heretofore, the clinical implications that have been inferred from perfusion CT studies have predominately relied on simple hypothesis testing approaches, which discard much of the acquired information. By way of contrast, more advanced statistical approaches that facilitate characterization of distinct perfusion signatures that are attributable to different types of tissues (eg, normal tissues vs different stages of tumor growth), as well as similar cluster patterns of responses within and between patients, could provide new insights for enhancing the use of this technology in diagnostic settings. Moreover, analytical tools that facilitate the study of CTp heterogeneity could be critical for the development and effective application of perfusion-based surrogate markers of therapeutic response.23

In this manuscript, we consider a functional data analytic approach24 and assume that the observations collected on each CTp parameter during the scan period are a realization, up to a measurement error, of an underlying perfusion process. The corresponding perfusion curves are estimated in a flexible way through the use of Gaussian processes (GP),25 which are characterized by possibly different parameters for the normal and tumor tissues. In order to take into account the heterogeneity of the functional CTp responses, we further propose to borrow strength in the estimation of the inferred time course patterns by using a Bayesian nonparametric mixture model, more specifically by employing a functional Dirichlet process prior.26–28 This modeling choice also allows us to achieve dimension reduction by clustering the observed dynamics of the CTp parameters into a smaller set of canonical response curves. On the contrary, estimation of the CTp parameters using any parametric higher order polynomial basis function would suffer from nonlocal influence, and thus would be unable to accommodate the inherently local temporal trends that are exhibited in the data.

Bayesian hierarchical models and Dirichlet Process (DP) mixtures have been successfully exploited to represent functional curves, eg, clustering spline coefficients in Bayesian multivariate adaptive regression splines models.29 In computer modeling30 and machine learning,25 Gaussian process realizations are often used as a basis to model random functions.31 Finite mixtures and Dirichlet process mixtures of Gaussian processes have also been proposed to model a sample of curves, possibly encoding complex spatiotemporal or covariate dependencies.32–37

With respect to those contributions, the method we propose for the analysis of functional CTp responses takes explicitly into account prior information about the perfusion imaging experiment and the vascular physiology of normal and tumor tissues. Because the time course mappings of CT perfusion parameters are expected to stabilize some time after injection of the contrast agent, in order to obtain a reliable quantification of the perfusion characteristics,21 our proposal accommodates change-points in the distribution of the functional responses over time, namely, in their temporal covariance structure. Moreover, since the perfusion time courses may vary across tissue types, we also allow for distinct change-points in the curves characterizing the normal and tumor tissues.

In addition, we demonstrate how to use the results from the posterior inference obtained from our Bayesian modeling framework to characterize ROIs from future patients as malignant tumor tissues to reflect diagnostic clinical settings. The classification is based on posterior predictive computations and the use of the Bayes factor, which weights the evidence in the data for the tested hypotheses.38 More specifically, we describe how the classification can be conducted by compounding the contribution of time courses from multiple perfusion parameters. Figure 2 summarizes the main features of our contribution.

Figure 2

General scheme of our modeling framework and inferential objectives.

The remainder of the manuscript is organized as follows. The “Methods” section describes the details of Bayesian nonparametric approach. Methods for posterior inference and classification are succinctly described in “Gaussian processes with varying autocorrelation” and “Posterior inference” sections. In the “Results” section, we describe the posterior and clustering results for analysis of tissue permeability in a sample of 16 patients with neuroendocrine liver metastases. We also provide an illustration of the classification performance of our algorithm on the basis of two studies: one aimed at comparing the predictive ability of our model with respect to commonly used classifiers in our small dataset, and a more comprehensive simulation study that further illustrates its general properties. We provide some concluding remarks in the “Discussion” section.

Methods

Bayesian mixture for modeling the heterogeneity of CTp characteristics

Without loss of generality in this section, we detail our modeling strategy for a single perfusion characteristic. More specifically, let ij = (Y(t1i), …, Y(T))T denote the vector of perfusion values collected on each ROI j = 1, …, n, and on each patient i = 1, …, n. Notice that we allow for scans obtained at different temporal frequencies across regions and patients. Furthermore, let z be a binary variable, such that z = 1 if the jth ROI corresponds to a tumor region, and z = 0 if the jth ROI is from normal tissue. In this section, we assume to know the tissue type (ie, tumor or normal) of each sampled region. However, this might not always be the case. In the “Classification performance” section, we provide a way to predict from the perfusion data if a tissue is normal or tumor based on the knowledge accumulated from previous data within the two conditions.

Following the typical functional data analysis framework,24 we envision that for either the normal or tumor tissue, the CTp responses can be formally described by the following model, where the ∈ij are independent realizations of a Gaussian white noise process with variance dependent on the tissue type, ie, , with d indicating the number of elements in the observation vector, Y, and I a d-dimensional identity matrix. Equivalently, given the θ, we can write

Arguably, the inferential interest often pertains to the modeling of the mean vector θ, which can be seen as the realization of a random function θ(t) on the observation points, ie, the data are supposed to be noisy realization of the underlying perfusion curves, y(t) = θ(t) + ∈(t), with , t ∈ ℝ+. A typical assumption specifies θ(t)’s as independent realizations of a Gaussian process. Instead, we propose to borrow strength in the estimation by introducing probabilistic dependence across the θ’s. More specifically, we assume that the θ’s are sampled from a prior (unknown) probability measure on ℝ+, which is also tissue-dependent, that is, . The vector of observed CTp parameter values can then be described by the following location mixture of Gaussian densities, where, for notational simplicity, we have used the symbol to also indicate the finite-dimensional distribution of the random function observed at time points t1, …, T.

Functional Dirichlet process priors

We follow a Bayesian nonparametric approach and further assign a prior probability to the mixing distributions , z = 0, 1. More specifically, we assume a functional Dirichlet process prior,27 by taking as the realization of a Dirichlet process.26 The Dirichlet process has been widely employed in Bayesian nonparametric (BNP) models due to its various properties. For example, the use of a Dirichlet process prior allows borrowing of information across observations in the estimation of the parameters of the model by automatically providing an unsupervised clustering of the data. More formally, the Dirichlet process formulation implies that is almost surely a discrete probability measure of the form, where δ denotes a probability measure degenerate on the atom . In our framework, this implies that concentrates probability masses p on a set of canonical perfusion curves in ℝ+, so that curves characterized by similar trajectories can be clustered together. Both the weights p’s and the curves are random functions, with probability laws specified as follows. For the sequence of weights (p1, p2, …), we assume a stick-breaking prior with parameter a, ie, p1 =V1, , Beta (1, α).39 The ’s are i.i.d. from a nonatomic probability measure G0 on Θ, independent of the p’s. In particular, G0 is commonly regarded as a parametric centering (base) distribution, since E(G) = G0, whereas the parameter α > 0 is a precision parameter, since it controls the variability of G around G0, with larger values of a resulting in realizations of G that are closer to G0. By assuming a random mixing distribution G with a DP prior, we do not restrict the model to a specific parametric form, and we increase its flexibility to capture different types of trajectories of the perfusion values. In the following, we denote our functional prior by G ∼ f DP(α, G0).

Gaussian processes with varying autocorrelation

We further take G0 as a stationary Gaussian process, that is, we assume that the distribution of the canonical curves on any finite set of time points in ℝ+ is Gaussian.25 Just as a multivariate Gaussian distribution is fully specified by its mean and covariance matrix, a Gaussian process is also specified by the mean and a covariance function. In symbols, we write , where θ0 may be either known (eg, a constant function) or itself assigned a prior. In the following, we assume θ0 ∼ N(0, η2), with η2 large so to encode vague information on the marginal expectation of the process. R(φ) represents the covariance function of the process, as a function of a set of parameters φ. The covariance function encodes the smoothness properties of the process (eg, mean square continuity and differentiability). A class of covariance functions that has proven to be attractive in various respects is the so-called Matérn covariance function.40 Let φ = (ψ, v τ), with ψ > 0, v > 0, and τ > 0, and let ∆ > 0 denote a time inter val of length ∆ from t. Then, the covariance function can be written as where Hν() is the modified Bessel function of order v.,41 ch. 9. The parameter ψ can be viewed as a decay parameter, since it governs the rate at which the covariance drops as a function of the distance, and v is a parameter that controls the degree of smoothness of the process. For , we get the exponential covariance function, R, (∆) = τ2 exp(−ψ∆). The resulting stochastic process is continuous but not differentiable at the origin; thus, it may be appropriate for modeling curves that can possibly vary abruptly in their gradients. Finally, τ represents a variance term. In Bayesian inference for Gaussian processes, it is often common to assume τ = 1, due to the poor identifiability of the model parameters, especially when the inferential interest is on estimating the correlation structure of the data.42

In order to obtain a reliable quantification of the perfusion characteristics, since the perfusion values are expected to show low autocorrelation soon after injection of the contrast agent, but subsequently they become more stable and highly correlated,21 we allow for changes in the correlation structure of the process over time. More specifically, we assume that the decay parameter in the correlation function may vary with time. Mathematically, we could describe these changes as different states s of the process, for k = 1, …, K (K finite). Each state would be characterized by a specific correlation function R(φ), and we would allow the process to move across the different correlation states by employing a general hidden Markov model framework.43,44 However, based on our knowledge of the perfusion characteristics, we only expect one change point during the scan period. Let indicate such (unknown) change point. Therefore, we consider two states s1 and s2 and assume that, for , the process cannot revert anymore to the preceding state. We further assume that the change point may be typically different for the normal and tumor tissues due to the different environment the contrast agent faces. Let , , where we set T = max for notational simplicity. Then, we assume that the base distribution of the functional Dirichlet prior (3) is characterized as where is a block diagonal matrix of the form with defined on T1 and defined on T2. Note that our formulation implicitly assumes conditional independence of the observations before and after the change point. For z = 0, 1, we can summarize our model as follows where , and is typically an inverse gamma prior, is a gamma distribution, and and are both (independent) uniform distributions.

Posterior inference

Lai and Xiang45 have recently considered a simple Bayesian model for multiple parameter changes in a multiparameter exponential family, developing explicit formulas for the estimators of the change-points. Our Dirichlet process formulation does not allow the use of the formulas of Lai and Xiang. Thus, we rely on a Markov chain Monte Carlo (MCMC) algorithm to obtain samples from the posterior distribution of the parameters of interest and conduct posterior inference. In this section, we provide a very brief description of the basic steps of the MCMC algorithm. Since perfusion parameters can be acquired at different time points across patients, we start by considering the union set of all observed time points. We then treat the values of the perfusion parameters at the times not observed as missing at random.46,47 Following the standard Bayesian approach,48 we impute the missing values within each MCMC iteration by drawing from the relevant posterior predictive distributions conditionally on the observed data and the currently sampled values of the parameters. Due to the normality and conjugacy of the model, the corresponding posterior predictive distribution is multivariate normal. More precisely, let denote the set of observed values for ROI j and patient i, and let indicate the missing values, which need to be imputed at each iteration. Let , with Hmis,mis, Hobs,obs, and Hmis,obs denoting the corresponding submatrices. Then, the distribution of conditional on is multivariate normal with mean and covariance matrix . The remainder of the model parameters is updated as described below, where, for notational simplicity, we omit the subscript z.

Sampling the canonical curves . We update the canonical curves using the algorithm in Ref. 49, which is based on the Pólya urn representation of the Dirichlet process.26,50 The canonical curves are then resampled within each cluster to improve the mixing of the chain, similar to the algorithm originally described in Ref. 51.

Updating the change point . We use a Metropolis–Hastings step to update t. Since the set of time points is discrete, we consider a multinomial proposal, such that the proposed change point is sampled with probability proportional to the values of the unnormalized likelihoods at each time t, calculated assuming .

Updating the covariance parameters . The update of depends on the particular form of the covariance function we specify. For the data analysis described in “Discussion” section, we use a Matérn with parameter v = 1/2. Then, the resulting exponential correlations are parameterized by the decay parameters and for the process in states s1 and s2, respectively. Conditional on the current sampled value of at the MCMC iteration, are updated with standard Metropolis–Hastings techniques.

The updates of the other model parameters are standard.52 For example, the full conditional of θ0 is multivariate normal and that of σ2 is an inverse gamma. The update of a follows the procedure described in Ref. 49. We refer to the Supplementary File for a more detailed description of the MCMC algorithm.

The posterior samples so obtained can be used to provide posterior estimates of the parameters of interest, as well as highest posterior credible intervals. For example, the MCMC sample average provides an ergodic estimate of the posterior expected value of a parameter. Let , b = 1, …, B, indicate the posterior draws of θ at time t after burn-in. Then . The accuracy of the estimate will increase as B, the length of the chain, increases. We found that in our applications, assuming B = 10, 000 was generally enough to produce good estimates. Similarly, we can acquire 95% posterior credible intervals of the parameters by estimating the posterior density of the MCMC draws and determining the corresponding 2.5% and 97.5% quantiles.

Classification via Bayes factors

We can use the information obtained from posterior inference to classify new observations, eg, from ROIs for which the tissue class is heretofore unknown. Let indicate the CT perfusion time courses, measured for a characteristic m = 1, …, M on a given ROI. A clinician may be interested to classify the corresponding tissue based on the values of and the knowledge gathered from previous datasets. This can be accomplished by taking advantage of the updating scheme typical of the Bayesian framework. More specifically, we can compute the Bayes factor [BayesF, 38] to discriminate between normal and tumor tissues. The Bayes factor is related to the posterior odds of one hypothesis relative to another (eg, tumor vs normal tissue), and it actually coincides with them if the prior probability of each hypothesis is the same. The log10 of the Bayes factor has been traditionally referred to as the weight of evidence in the data for the tested hypothesis.53 Then, if log10(BayesF) > 0, the data indicate some evidence against the hypothesis, with values log10(BayesF) > 0.5 and higher denoting substantial and increasing evidence. Let znew be a binary indicator, such that znew = 1 if the newly examined tissue is classified as tumor, and znew = 0 if normal. Then, we can compute the Bayes factor for each characteristic m = 1, …, M as follows, where the quantities and Pr(θ| data) denote, respectively, the posterior predictive density of and the posterior probability of the model parameters if we assume that the data are from tissue type z, based on the already available CT perfusion data.

For clinical purposes, the classification of targeted ROIs shall rely on the Bayes factor computed either on a single perfusion characteristic (eg, blood flow) or from a combination of the perfusion characteristics, eg, only those which have been previously found in the literature to be associated with angiogenesis and tumor growth. For M different characteristics, the decision will be based on the evaluation of the following compound log-Bayes factor, where w ≥ 0 represent weights that can be chosen to reflect the relative importance that a clinician might attribute to the different perfusion parameters in the decision process. If w = 0, then the mth perfusion characteristic will not contribute to the classification. Due to the lack of available prior information, the choice w1 = w2 = … = w = 1 effectuates equivalence among the perfusion parameters. The choice of the weights could also be informed by existing knowledge about the dependence between perfusion characteristics, eg, by choosing weights inversely proportional to the strength of association of each CTp parameter with the others.

Results

A case study from patients with neuroendocrine liver metastases

In this section, we apply our Bayesian nonparametric approach to a retrospective clinical study. Both our study and the retrospective clinical study were approved by the MD Anderson Institutional Review Board (IRB). The study included 16 patients (6 men and 10 women) with metastatic neuroendocrine tumors who received optional CT perfusion imaging between April 2007 and September 2009 as part of two IRB approved clinical trials. CT perfusion was performed for a target lesion in the liver, which was clinically or radiologically determined to be malignant. Images were obtained with a 64-row multidetector CT scanner (VCT; GE Healthcare) and with acquisition time of 12–590 seconds.21 We consider the time courses of four perfusion characteristics, as follows: BF, BV, MTT, and PS. The measurements were obtained on 27 separate ROIs in liver metastases and 25 separate ROIs in normal liver tissue. More specifically, for each of the eight axial slice locations of each dataset, a liver tumor ROI was drawn free hand around the periphery of the primary target lesion, using an electronic cursor and mouse, with reference to the source cine CT images and perfusion parametric maps, displaying the images at soft tissue windows (width = 350 HU, level = 40 HU). Wherever possible, a second tumor ROI was delineated, provided it was greater than 1.5 cm in diameter. ROIs were placed in the abdominal aorta and in the portal vein on the source images to provide these vascular inputs. Refer to Ref. 21 for more details. In this section, we illustrate the results of our analysis for the time courses of the PS characteristic, since PS values have been previously associated with the markers of angiogenesis.54 More specifically, we apply model Equation (5) to a normalizing log-transformation of the PS perfusion values, and then perform posterior inference as described in the “Posterior inference” section. The log scale is frequently used in CT perfusion analysis to adjust for conditionally asymmetric residual error at a given acquisition time and to alleviate the effect of heteroskedasticity as a function of time. In terms of prior choices, we assume that the mean of the base G0 , is the null vector 0 with diagonal variance–covariance matrix 4I. Given the range of values observed in CT perfusion data, this setting corresponds to a vague (noninformative) prior on the elements of . The prior on the range parameter of the Matérn covariance function, , is taken to be uniform over the values from 0 to 10. The prior specification is completed by assuming a gamma prior α ∼ Ga(l, 1) on the precision para meter of the functional Dirichlet process, a vague inverse gamma prior on the sampling variance, , and a uniform prior across the full range observed time domain (0, 600) for the change point, .

Figure 3 displays the posterior distributions of the number of clusters for the normal (left panel) and the tumor liver tissue (right panel). In both cases, the posterior is characterized by high variance and long tails. However, the modal number of clusters is estimated around 5 in the normal tissue and 6 in the tumor tissue type. Hence, we may conclude that there is evidence of heterogeneous clustering patterns among the CT perfusion time courses. The MCMC algorithm typically samples possibly different cluster configurations at each iteration. Several methods have been used in the literature to obtain a single-point estimate of clustering draws from the posterior distribution.55,56 In the following, we consider the maximum a posteriori (MAP) estimate. The MAP clustering is the clustering estimate that maximizes the posterior density.

Figure 3

Posterior of the number of clusters of PS time courses for the normal (left) and tumor (right) liver ROIs.

Figure 4 (top) displays a cubic spline interpolation of the original PS values, with different colors to highlight the different MAP clusters. In both normal and tumor samples, the clusters appear to represent different time course trajectories, either because of their behavior within the first 300 seconds of acquisition or because of the resultant level that they achieve at the end of the scanning period. Figure 5 shows the posterior distribution of the change point in the covariance structure of the CT perfusion values. It is well known that for many of the correlation functions typically used for Gaussian processes, the decay parameters are only weakly identifiable.57,58 Accordingly, the posterior distribution of the change-points is quite diffuse. Nevertheless, one can recognize a mode within the first 100 seconds of the acquisition time. More specifically, the posterior expectation is 89.6 seconds (67.8 SD) for the normal tissue and 93.6 seconds (68.4 SD) for the tumor tissue. These values reflect the heteroskedasticity patterns typically observed in CT perfusion values.

Figure 4

Posterior clustering: cubic spline interpolation of the observed log-PS values, color coded according to the MAP estimate for normal tissue type (left) and tumor tissue type (right).

Figure 5

Posterior distribution of the change point in the correlation structure of log PS values for the normal (left) and tumor (right) liver tissues.

Classification performance

In order to assess the classification performance of our Bayesian nonparametric modeling approach, we devise two different strategies. First, we consider the dataset analyzed in “A case study from patients with neuroendocrine liver metastases” section and we select one normal and one tumor curve for each of the 16 patients. We then use formula (6) to classify the selected curves either as normal or tumor, separately for each perfusion parameter. Ng et al.21 have recently discussed the minimal acquisition duration that is necessary to attain stability and good quantification of the CT perfusion parameter values. There is clear motivation to reduce the overall duration of CT acquisition to the shortest time possible to reduce radiation exposure without compromising quantification of the CT perfusion parameter values. A few studies have suggested that acquisition times of 30–60 seconds might be satisfactory for some of the CT perfusion parameters.59,60 However, Ng et al.21 suggest that for most parameters stabilization with moderate confidence is attained only between 220 and 360 seconds of acquisition. Those studies did not address the predictive ability of the measurements collected over the different acquisition times to discriminate between tissue types. For this reason, in the following, we will evaluate the predictive ability of observations collected over 30–100 seconds in addition to the full-time courses over the complete 590 seconds of acquisition duration.

For comparison, we also implement the commonly used classifiers, linear discriminant analysis (LDA), quadratic discriminant analysis (QDA), and support vector machine (SVM) algorithms. LDA assumes that the measurements from each class are normally distributed with the same variance–covariance matrix for both tumor and normal tissues, whereas QDA relaxes the homoskedasticity assumption. As these classifiers require the same number of time points for each observation, we first fit cubic spline interpolations to the data in order to be able to evaluate all the classifiers on the same acquisition intervals. This is particularly true for the PS and the BF parameter, which is consistent with previous literature. The discriminant analysis techniques (run after curve interpolation) do not appear to perform consistently well, especially when using the first 30–100 seconds of the acquisition scans. This is particularly relevant considering that higher acquisition duration would require high overall radiation exposures.

Table 1 shows the true negative and true discover rates for the prediction of normal and tumor curves, separately for the four CT perfusion parameters. The BNP method we propose generally seems to perform better than the alternative methods in most cases.

Table 1

Prediction true negative (normal) and true discover (tumor) rates using LDA, QDA, and the proposed BNP between 30–100 seconds and 0–590 seconds.

	LDA	LDA	QDA	QDA	SVM	SVM	BNP	BNP
Normal	(30–100 sec)	(0–590 sec)	(30–100 sec)	(0–590 sec)	(30–100 sec)	(0–590 sec)	(30–100 sec)	(0–590 sec)
logBF	75%	94%	81%	94%	75%	63%	100%	100%
logBV	69%	88%	63%	88%	75%	50%	94%	75%
logMTT	75%	88%	69%	94%	50%	56%	100%	100%
logPS	69%	88%	88%	94%	69%	69%	100%	100%
Tumor
logBF	50%	88%	56%	94%	63%	56%	75%	56%
logBV	81%	88%	81%	88%	75%	50%	69%	56%
logMTT	69%	88%	88%	82%	50%	56%	82%	63%
logPS	94%	100%	94%	94%	69%	56%	100%	82%

To further test our classification method, we also devise a simulation strategy where we consider a larger dataset, obtained by generating 1,000 distinct time courses for each of the four CT perfusion parameters, as follows. First, we obtain spline interpolations of the observed time courses in the liver metastases case study considered in “A case study from patients with neuroendocrine liver metastases” section. The spline interpolations are fit using the R package spline on each of the 52 curves and for each of the four characteristics.61 Then, we select one of the 52 original observations by sampling with repetition, separately for each characteristics. A new time course is generated as a draw from a multivariate normal distribution with the selected interpolation as the mean and a diagonal covariance matrix, with variance η2. This simulation scheme allows us to evaluate the effect of the noise η2 on the prediction and also to more accurately assess the asymptotic performance of our classifier under multiple realizations of the generating models.

Tables 2 and 3 report the results of the classification for varying degrees of the noise η2 (η2 = 0.1, 0.2, 0.5) using the first time frame of 30–100 seconds of acquisition duration and the full scanning period, respectively. The column corresponding to Combination refers to the compound log-Bayes factor in equation (7), where we assume all weights w1 = w2 = w3 = w4 = 1. Such a choice can be considered as default, unless prior information leads to alternative choices. As expected, the prediction accuracy decreases with higher values of noise. The best predictor appears to be represented by the PS values, whose performance also appears more robust to higher levels of noise. Interestingly enough, the accuracy of the prediction with PS values increases slightly when using the full acquisition duration from 0 to 590 seconds, whereas for all other perfusion characteristics, a relatively better performance is obtained by considering the values between 30 and 100 seconds. This is consistent with the findings in the study by Ng et al.21, where it is shown that the PS values require longer acquisition times than do BF, BV, and MTT to reach comparable levels of stabilization in CT perfusion values. Among the other parameters, the BF and the MTT also appear to be relatively robust to higher values of noise. BV appears to be the less informative of the parameters in terms of prediction performance. Finally, the compound Bayes factor (7) appears to be the best classifier only when considering the first 30–100 seconds of acquisition duration. If the full-time course is considered, the performance is notably reduced. This result appears to be mainly due to the negative accuracy of the BV, and therefore, one might consider assigning a weight w close to zero to this perfusion parameter.

Table 2

Classification results for the large simulation with 1,000 simulated time courses described in the “Classification performance” section, for varying degree of noise ŋ2 = 0.1, 0.2, 0.5 and considering the first 30–100 seconds of acquisition duration.

30–100 SECS	BF	BV	MTT	PS	COMBINATION
ŋ2 = 0.1
Normal	342/503 (68%)	332/503 (66%)	397/503 (79%)	412/503 (82%)	423/503 (84%)
	<0.001	<0.001	0.03	0.36
Tumor	357/497 (72%)	333/497 (67%)	378/497 (76%)	402/497 (81%)	442/497 (89%)
	<0.001	<0.001	<0.001	0.07
ŋ2 = 0.2
Normal	337/503 (67%)	322/503 (64%)	377/503 (75%)	357/503 (71%)	392/503 (78%)
	<0.001	<0.001	<0.001	0.01
Tumor	338/497 (68%)	323/497 (65%)	347/497 (70%)	400/497 (80%)	407/497 (82%)
	<0.001	<0.001	<0.001	0.57
ŋ2 = 0.5
Normal	287/503 (57%)	206/503 (41%)	292/503 (58%)	312/503 (62%)	347/503 (69%)
	<0.001	<0.001	<0.001	0.02
Tumor	293/497 (59%)	193/497 (39%)	293/497 (59%)	318/497 (64%)	359/497 72%)
	<0.001	<0.001	<0.001	0.005

Note: The P-values of the comparison of each individual CT characteristics and the combination (two-sample test for equality of proportions) are reported under each result.

Table 3

Classification results for the large simulation with 1,000 simulated time courses described in the “Classification performance” section, for varying degree of noise ŋ2 = 0.1, 0.2, 0.5 and considering the full acquisition duration.

0–590 SECS	BF	BV	MTT	PS	COMBINATION
ŋ2 = 0.1
Normal	331/503 (66%)	236/503 (47%)	377/503 (75%)	448/503 (89%)	247/503 (49%)
	<0.001	0.49	<0.001	<0.001
Tumor	338/497 (68%)	203/497 (41%)	362/497 (73%)	457/497 (92%)	268/497 (54%)
	<0.001	<0.001	<0.001	<0.001
ŋ2 = 0.2
Normal	312/503 (62%)	221/503 (44%)	352/503 (70%)	411/503 (82%)	226/503 (45%)
	<0.001	0.75	<0.001	<0.001
Tumor	313/497 (63%)	189/497 (38%)	337/497 (68%)	442/497 (89%)	238/497 (54%)
	<0.001	0.002	<0.001	<0.001
ŋ2 = 0.5
Normal	261/503 (52%)	150/503 (30%)	277/503 (55%)	352/503 (70%)	211/503 (42%)
	0.002	<0.001	<0.001	<0.001
Tumor	278/497 (56%)	129/497 (26%)	278/497 (56%)	358/497 (72%)	249/497 (50%)
	0.07	<0.001	0.07	<0.001

Note: The P-values of the comparison of each individual CT characteristics and the combination (two-sample test for equality of proportions) are reported under each result.

Discussion

With the advances of the contrast-enhanced functional imaging technology, the development of noninvasive perfusion imaging biomarkers for tissue characterization, cancer prognostication, and detection has emerged as an area of recent focus in clinical cancer research. Moreover, these imaging features have the potential to enhance quantitative evaluations for measuring therapeutic response to antiangiogenic treatment strategies. Future endeavors to further develop and translate this technology should rely on appropriate analytical models for characterizing the multiple sources of variance that are inherent to the acquisition and measurement of temporal changes in contrast enhancement obtained from dynamic CT.

In this article, we presented a Bayesian nonparametric functional analytic approach for the analysis of CT perfusion time courses, which allows for heterogeneity observed across patients and tissues, by clustering the measurements of perfusion characteristics into groups characterized by similar temporal behavior. Our approach takes explicitly into account prior information about the perfusion imaging experiment and the vascular physiology of normal and tumor tissues, by accommodating change-points in the temporal covariance structure of responses over time in order to appropriately describe the temporal heteroskedasticity observed in such measurements. We demonstrate how our approach can lead to improved performance with respect to commonly used classification methods for classifying ROIs in future patients, especially when considering the typical length of acquisition times. Therefore, our approach can assist in discriminating malignant from healthy tissue regions in diagnostic settings.

The clinical study was limited only to metastases to the liver from a specific tumor (metastases arising from neuroendocrine tumors) and consisted of a relatively small number of patients. The extent to which the conclusions of our case study might be generalizable to other tissues and tumors requires future exploration. However, we expect that the proposed methodology for functional data classification is generalizable to any perfusion study. Our data were obtained using a specific acquisition protocol based on relatively high temporal sampling initially (in the first 30s), and more sparsely sampled data in its second phase (out to 590s). While one would not expect more widely spaced temporal sampling in the second phase to have substantial impact on analyses, an ideal dataset might be the one that was acquired at a high temporal sampling in the second phase as well. However, obtaining such data would require high overall radiation exposures and present practical challenges in acquisition because of the constraints of breathing motion and registration.

Furthermore, we should note that our study does not characterize treatment-induced changes in the perfusion values. Using data from a larger randomized clinical trial, the clustering detection could be informed by treatment information, eg, by using an ANOVA-dependent Dirichlet process.62 Finally, our modeling formulation is more computationally demanding than competing approaches. While our extensive simulation studies suggest that computations do not slow down substantially for datasets comprising a few hundred patients’ samples, for larger datasets, computational speed could be improved through the use of variational Bayes approaches.63 The precision of the classification, when considering the full acquisition duration, is dramatically decreased, as presented in Table 3. Scarpa and Dunson64 have recently proposed enriched stick-breaking processes for functional data, which enable incorporation of prior information about attributes of the curves in the classification. A similar feature selection-based strategy could be used to guide the classification of the curves so that specific characteristics of each perfusion parameter are taken into account and/or relevant portions of the time courses can contribute to the inference more than others.

Additionally, while our framework considered ROI-level inference, which is perhaps most prevalent in clinical practice, the methodology could be adapted to accommodate voxelwise analysis following image registration. The methods, when implemented in this context, could be leveraged to create spatial–temporal posterior probability maps, providing diagnostic tools that could potentially be used to enhance existing tumor segmentation approaches.

Supplementary File. MCMC details.