Modeling and Reconstruction of Mixed Functional and Molecular Patterns

Functional medical imaging promises powerful tools for the visualization and elucidation of important disease-causing biological processes in living tissue. Recent research aims to dissect the distribution or expression of multiple biomarkers associated with disease progression or response, where the signals often represent a composite of more than one distinct source independent of spatial resolution. Formulating the task as a blind source separation or composite signal factorization problem, we report here a statistically principled method for modeling and reconstruction of mixed functional or molecular patterns. The computational algorithm is based on a latent variable model whose parameters are estimated using clustered component analysis. We demonstrate the principle and performance of the approaches on the breast cancer data sets acquired by dynamic contrast-enhanced magnetic resonance imaging.


INTRODUCTION
Functional imagingtechnologies are providing researchers and physicians with exciting new tools to study important disease-causing biological processes in living tissue [1,2]. Dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) uses various molecular weight contrast agents to assess tumor vascular permeability and quantify cellular and molecular abnormalities in blood vessel walls [3]. DCE-MRI can characterize vascular heterogeneity and elucidate features that distinguish angiogenic blood vessels from their normal counterparts, and has potential utility in assessing the efficacy of angiogenesis inhibitors in cancer treatment [2][3][4]. Although DCE-MRI can provide a meaningful estimation of vascular permeability when a tumor is homogeneous, many malignant tumors show markedly heterogeneous areas of permeability and vascular endothelial growth factor (VEGF) expression; thus the signal of each pixel often reflects multiple microenvironments in a tumor representing a complex summation of vascular permeability with various diffusion rates [2,3].
Several quantitative methods based on parametric compartment modeling (CM) have been developed to dissect the spatial distribution of vascular heterogeneity associated with tumor angiogenesis [3,5,6]. These methods estimate the fundamental kinetics components (called factors) and the associated factor weights (called factor images) [6][7][8]. Each factor is interpreted as the time course of a compartment, whereas each factor image is interpreted as local weights representing the spatial distribution of vascular permeability with different diffusion rates [2,9]. The parametric model chosen may not fit the data obtained, and each model makes a number of assumptions that may not be valid for every tissue or tumor type. The causes for modeling failures are complex and often not well understood [6,10]. Key reasons include multiple tissue compartments, an incorrect arterial input function, and numerical nonidentifiability of the parametric model [3,6,[9][10][11]. This motivates the consideration of clustered component analysis (CCA) that can be based on a flexible compartment latent variable model [12,13]. The objective is to factorize the underlying angiogenic permeability distributions (APD) and time activity 2 International Journal of Biomedical Imaging curves (TAC) from dynamically mixed DCE-MRI image sequences [14,15].

THEORY AND METHODS
We first introduce a simple form of compartment latent variable model for DCE-MRI. Without loss of generality, we initially focus on the two-tissue compartment model shown in Figure 1. The tracer characterization within a region of interest can be approximated by a set of first-order differential equations [15]:ċ where c f (t) and c s (t) are the tissue activity in the fast turnover and slow turnover pools, respectively, at time t; c p (t) is the tracer concentration in plasma (the input function); c(t) is the measured total tissue activity; k 1 f and k 1s are the unidirectional transport constants from plasma to tissue (permeability in ml/min/g: spatially varying); k 2 f and k 2s are the rate constants for efflux (diffusion in /min: spatially invariant) in the fast flow and slow flow pools, respectively [5][6][7].
Mathematical consideration based on a latent variable model suggests a simple method to convert temporal kinetics (1) to spatial information [7]. Let ] T be the observed tracer activities of pixel i measured at L time points. Now consider a source vector of spatial permeability distributions k(i) = [k f (i), k s (i), k p (i)] T together with an L×3 mixing matrix A(t) which maps the latent space into the data space: where the TACs associated with different APDs are and ⊗ denotes the convolution operation. Relationship (2) describes how the observed multivariate data are generated by a process of mixing the latent components, as illustrated in Figure 2.
Since the APDs k(i) and TACs A(t) are both unknown, what we seek in the above model is an algorithm that can perform blind source separation to recover the source patterns from their observed mixtures. Based on the realistic assumption that the APDs are spatially heterogeneous (e.g., piecewise stationary with insignificant spatial overlap) [2,3,16], CCA on x(i) over the time domain aims to perform a nonparametric multivariate clustering of pixel TACs similarly to the successful application in functional MRI analysis [13].
Intuitively, when there are only pure-volume pixels, a oneto-one association between pixel TAC x(i, t) and one of the source TACs a j (t) exists-except for a local scaling by k j (i) and some additive statistical variation where ε(i) is the noise term being both temporally and spatially white Gaussian distributed with zero mean and unknown variance σ 2 x(i) and a j = [a j (t 1 ), a j (t 2 ), . . . , a j (t L )] T [13]. However, note that the forms of compartment TACs in (4) are no longer necessarily parametric as in (3) and are much more flexible; this should help reduce the potential for modeling failures. To perform a top-down CCA on x(i) to estimate a j based on (4), the shape rather than the magnitude k j (i) of the pixel TAC is of the interest [7,11,17]. By performing both "centering" and "normalization" over time, given by each pixel TAC can be transformed to a constant scale with mean zero independent of amplitude variations, denoted by the normalized x n (i). There has been considerable success in using the standard finite normal mixture (SFNM) distribution to model clustered data sets, taking a sum of the following general form [18]: where π is the mixing factor and g is the Gaussian kernel with mean a j and covariance matrix C j .
Finding an estimate of the mixing matrix A comes down to performing a maximum likelihood estimation of the SFNM model (6), where the joint log-likelihood is given by where N is the number of the pixels. This clustered component analysis task can be, fortunately, solved by the expectation-maximization (EM) algorithm that maximizes the joint log-likelihood [13,15,18] π j , a j , Yue Wang et al.

3
Fast flow pool c f (t) Slow flow pool c s (t) Vascular space plasma c p (t)  where the "soft" splits of a pixel TAC allow x n (i) to contribute simultaneously to multiple source TACs. Specifically, in order to compute the expectation step of the EM algorithm, we must first estimate the posterior probability that each pixel TAC x n (i) is of source TAC a j , namely, the maximization step of the EM algorithm [13]. Such estimated posterior Bayes probabilities of pixel TAC x n (i) associated with one of the source TACs are given by (9) and the compartment TACs are the normalized and weighted sample averages of pixel TACs in the light of their compartment memberships estimated by (9), computed via the expectation step of the EM algorithm [13]: for j = f , s, p. Having determined the mixing matrix A = [a f , a s , a p ] representing the compartment TACs, the APDs can be reconstructed using a least squares fit according to (2), resulting in where T denotes matrix transpose.
To perform CCA, we use the visual statistical data analyzer (VISDA) algorithm [18]. The main function of VISDA is cluster modeling, discovery, and visualization. In addition to the multivariate soft clustering by the EM algorithm, VISDA also includes model selection by minimum description length (MDL) criterion and cluster initialization by hierarchical clustering. To capture all of the hidden clusters, VISDA is both statistically principled and visually insightful that incorporates both the power of statistical methods and the human gift for pattern recognition. VISDA uses an adaptive boosting of discriminatory subspaces involving hierarchical mixture modeling, selected optimally by the MDL criterion, and allows the complete data set to be visualized at the top level and so partitions data set, with clusters and subclusters of data points visualized at deeper levels. Complementary to (9) and (10), the M step in EM algorithm also involves the update rules for cluster factors and covariance matrices for j = f , s, p. The E step involves assigning to the clusters, probabilistically, contributions from the data points and the M step involves re-estimating the parameters of the clusters in the light of this assignment. The algorithm cycles back and forth until the joint likelihood function is maximized. When there are multiple compartment mixture regions, one remaining issue in pixel TAC clustering is the model selection that refers to the detection of cluster number K 0 . The EM model fitting cannot be used to estimate K 0 since the ML is a nondecreasing function of K 0 , thereby making it useless as a model selection criterion. This problem can be, fortunately, solved by using MDL criterion in conjunction with EM clustering. MDL is a proven information-theoretic criterion for model selection and has proven asymptotically consistent. The major thrust of MDL-based cluster validation has been the formulation of a model fitting procedure in which an optimal model is selected from the several competing candidates such that the selected model best fits the observed data. Specifically, the optimal value of K 0 is selected by minimizing where the first term on the right is the approximation error and the second term on the right is the estimation error whose role is to penalize large value of K 0 . For K 0 = K min , . . . , K max , the values of MDL are calculated and a model with K 0 clusters is selected that will correspond to the minimum MDL value.

EXPERIMENT AND RESULTS
In this section, we first demonstrate the performance of clustered component analysis when applied to real DCE-MRI data sets. The data was acquired at the NIH Clinical Center using gadolinium DTPA as the contrast agent. The threedimensional DCE-MRI scans were performed every 30 seconds for a total of 11 minutes after the injection. For the purpose of comparison, Figure 3 shows the estimated TACs associated with the input function as well as the fast and slow flows obtained by the advanced parametric compartment modeling method. The corresponding reconstructed APDs are given in Figure 4. This represents an advanced breast tumor case where active angiogenesis occurs often in the peripheral area (i.e., boundary with fast flow), while the inner core reflects hypoxia (dominated by slow flow) [2,3]. We then apply CCA to the same data set. The DCE-MRI sequence contains a total of 18 images taken at different times, of which, we remove the first few images that do not show sufficient contrast accumulation and use the remaining 12-15 images in the experiment. After an analysis by VISDA, MDL criterion determines that there is clearly more than one pixel TAC cluster. By targeting the two major compartment sources, the corresponding source TACs are estimated and the hidden APD images are subsequently reconstructed. Figure 5 shows the extracted source images (i.e., vascular permeability) carrying out fast and slow diffusions. The projected distribution of pixel TACs clearly re-Yue Wang et al. veals a multicluster data structure, and the scatter plot of the source images shows the expected globally negative correlation dependence.In addition, the ability of estimating the input function has aprofound impact in multivariate quan-tifications, since otherwise blood samples will be taken invasively at the radial artery or from an arterialized vein which, however, poses health risks and is not compatible with clinical practice. The outcome of CCA on further decomposing the mixture into the three underlying compartments is given in Figure 6 that presents a very consistent result with the one obtained by the independent method based on the compartment modeling shown in Figures 3 and 4. To test the stability of the performance by CCA, we have applied this method to a series of realistically simulated data sets in which various realistic TACs are numerically synthesized and the mixed observations are generated by weighting the real APDs (e.g., given in Figures 4 and 6) by these synthesized TACs. The results show that CCA can successfully reconstruct the hidden clustered components under various mixtures and TAC conditions.
As an example of more challenging problems, with significant practical utility, we report the preliminary application of CCA method in a longitudinal study to monitor a breast tumor's response to anti-angiogenic therapy. Defective endothelial barrier function due to vascular endothelial growth factor (VEGF) expression is one of the best-documented abnormalities of tumor vessels, resulting in spatially heterogeneous high microvascular permeability to macromolecules. Initial results suggest that changes in vascular permeability and volume fraction can be detected in a responsive tumor soon after therapy begins. Vascular permeability has been reported to correlate closely with VEGF expression in tumors, and decrease significantly after anti-VEGF antibody treatment and after the administration of other inhibitors of VEGF signaling. In breast and cervical cancers, a decrease in transendothelial permeability often accompanies tumor's response to chemotherapy and an early increase or no change in permeability can predict non-responsiveness or poorer prognosis.
Three sets of DCE-MRI data were acquired before and during the treatment period, each with three-months apart. Figure 7 shows the DCE-MRI images as a potential endpoint in assessing the response to therapy. The introduction of imaging cancer therapies by DCE-MRI has posed new challenges to traditional anatomic imaging approaches, because the vascularity of a tumor can change without a corresponding change in tumor size and vice versa. Our preliminary experiment shows promising results on the application of CCA to this problem, see Figure 8. For example, the extracted source TACs closely resemble the expected compartmental kinetics of the contrast agent, and both the APD images and TACs show the expected changes of the patterns over time, consistent with clinical assessment of a responsive case. Of particular scientific value, our results show that tumorinduced vascular activities were significantly reduced after a positive response to anti-angiogenesis chemotherapy, despite a noticeable increase in tumor volume during the initial treatment period.

DISCUSSION AND FUTURE WORK
In CCA approach, the significant overlap between the source image boundaries in space causes a potential partial volume effect (PVE) [16]. It can be shown that PVE will lead to a biased estimation of the compartment TACs. We can incorporate such PVE into the SFNM model that can be, fortunately, estimated by a constrained EM algorithm [16]. Specifically, we will first apply MDL criterion (13) to estimate the most appropriate number of distinctive temporal clusters that will also include the so-called composite boundary clusters [16]; we will then estimate the PVE-SFNM model by only updating the parameters of pure-volume clusters usingz i j , followed by the assignment of the parameter values for the partial volume clusters based on the PVE model [16]. Alternatively, we can simply consider those composite boundary clusters as intrinsic compartment TACs where the emphasis will be on . The pattern changes are consistent with the clinical foundings that, even as a responsive case, most tumors' volume will grow initially but shrink lately with much reduced vascular activities.
interpreting such compartments in relation to biological or clinical parameters. In longitudinal studies, special care should be taken in assessing the response to therapy. Inour present experiment, we have considered longitudinal samplings separately and performed blind source separation for each of the samplings. It would be more meaningful to consider the estimated TAC before the therapy starts as a baseline reference, and then estimate the source images in the follow-up studies to see whether the spatial distribution of fast and slow permeabilities changes. We can also use the estimated baseline source image as the reference, and subsequently recover the TACs in the follow-up studies to detect the changes of diffusion rates. The multivariate quantifications that reflect the efficacy of angiogenesis inhibitors has great potential but is at an early stage. Part of the challenge stems from incomplete knowledge of how blood vessels are affected. For example, angiogenesis inhibitors can block the growth of new blood vessels from existing vessels, but may also eliminate certain existing vessels, such as tumor vessels.
We believe that our comparative studies provide useful information on the utility of the proposed methods for computed simultaneous imaging of multiple functional or molecular biomarkers. Given the difficulty of the task, while the optimality of these methods may be data or modality dependent, we wouldexpect them to be important toolsin dynamic image formation and analysis. For example, since angiogenesis is complex process critical to growth and metastasis of malignant tumors, the clinical value and promise of DCE-MRI in imaging tumor angiogenesis before and during therapy provides strong incentive for advancing the imaging formation method [2][3][4]. Specifically, with the prior information on the nonnegativity of the mixing matrix and sources, new principle and perhaps improved methods may yet become possible [19]. Here we wish to propose a nonnegative least-correlated component analysis (nLCA) when the hidden sources and mixing matrix are known to be nonnegative [20]. This concept has powerful features which are of considerable universal applicability since it eliminates the condition of source independence and non-Gaussianity required by independent component analysis [19]. It can be shown that when the mixing matrix is nonnegative, the correlation between the mixtures is always positively increased, namely, the correlation increase theorem [20]. Such a positive increase in correlation after nonnegative mixing immediately suggests a possible recovering mechanism for blind source separation of dependent sources. With the encouraging preliminary success tested on real data sets, we are currently investigating the existence and uniqueness of nLCA solution that exploits the nonnegativity constraint on correlated yet well-grounded sources in the light of correlation increase theorem [19,20].