Blind source separation of hemodynamics from magnetic resonance perfusion brain images using independent factor analysis.

Perfusion magnetic resonance brain imaging induces temporal signal changes on brain tissues, manifesting distinct blood-supply patterns for the profound analysis of cerebral hemodynamics. We employed independent factor analysis to blindly separate such dynamic images into different maps, that is, artery, gray matter, white matter, vein and sinus, and choroid plexus, in conjunction with corresponding signal-time curves. The averaged signal-time curve on the segmented arterial area was further used to calculate the relative cerebral blood volume (rCBV), relative cerebral blood flow (rCBF), and mean transit time (MTT). The averaged ratios for rCBV, rCBF, and MTT between gray and white matters for normal subjects were congruent with those in the literature.

1. Introduction

Dynamic susceptibility-contrast perfusion magnetic resonance (MR) imaging is a commonly used method for the noninvasive assessment of cerebral blood perfusion because of the availability of MR imaging units and lack of exposure to ionizing radiation [1-10]. After a bolus injection of intravascular contrast agent, the passage of bolus induces the susceptibility inhomogeneity, which in turn causes a relative decrease of image intensities of brain tissues from the baseline. Various tissues manifest distinct blood-supply patterns since the contrast agent arrives consecutively, leading to temporal intensity drops at different time instants. Based on the perfusion data, we can calculate the hemodynamic parametric maps, that is, relative cerebral blood volume (rCBV), relative cerebral blood flow (rCBF), and mean transit time (MTT), by employing the indicator dilution theory [11, 12]. This has been used in the assessment of many brain disorders such as tumors [2, 7, 13], stroke [14], infection [15], and moyamoya disease [16, 17].

The estimation of rCBV and rCBF, however, requires arterial concentration of the contrast agent as an arterial input function (AIF). This is a demanding task, and many methods have been proposed to address the issue. Additionally, classification of blood-supply patterns for various tissue types in brain based on bolus transit profiles is essential for the assessment of brain perfusion. Wiart et al. [18] manually selected single or multiple pixels for the tissue of interest and used the averaged temporal profile as a reference in producing a similarity map for segmenting gray matter (GM) and white matter (WM) from perfusion images. This method is advantageous for easy implementation but limited to extraction of a single tissue type per similarity map. If one attempts to segment multiple tissue types, associated similarity maps need to be created one by one after cumbersome selections of reference pixels, which may be prone to operator influence. Alternatively, Martel et al. [19] used the conventional factor analysis, under the hypothesis that the observable signal-time curve for each pixel is a weighted sum of pure physiological factors, to classify the signal-time curves of arterial and venous structures from MR dynamic perfusion images. This method has also been applied in nuclear medicine to separate cardiac components and extract left ventricular input function from dynamic H2 15O PET images [20-22]. It should be noted that, although the factor-analysis-based methods are attractive, additional assumptions of a priori knowledge are required to obliquely rotate the eigenvectors to yield meaningful physiological factors [23-25]. Besides, only one or two physiological factors were resolved in past studies. To our knowledge, the automatic simultaneous segmentation of multiple perfusion compartments has been less explored. A recent related work was to use the expectation-maximization algorithm with a mixture of multivariate Gaussian models for fitting the perfusion MR images so that each pixel can be labeled by the resulting maximal posterior probability [26].

The aim of this study is to automatically classify the spatiotemporal blood-supply patterns which enables us to extract the arterial compartment for modeling the AIF and segment other tissue regions as well. To this end, we employ an independent factor analysis (IFA) [27], a data-driven method, allowing us to blindly separate mixed signals into independent-factor (or independent-source) components for multitissue hemodynamic classification. That is, the hemodynamics of each tissue type can be dissected without making a priori spatial and temporal assumptions of physiology. The factors in the IFA, in contrast to the conventional factor analysis, are modeled by a finite mixture of Gaussian functions that can be used as a constraint to remove rotational factors [27]. This method has been applied to successfully separate the background factors and noise artifacts from the stimulus-evoked MEG and EEG sensor data contaminated by large background brain activity [28]. In this study, the use of IFA is based on two inherent assumptions: (1) signal intensity of each pixel is a linear mixture contributed from different tissues, referred to as the partial volume mixing, which is a well-known phenomenon due to finite resolution of MR scan; (2) the anatomic structures of pure tissue types are spatially independent (nonoverlapped with each other). The classification of multitissue hemodynamics consists of two steps. The first step is to identify a dominant tissue type on each independent-factor (IF) image resolved from the IFA based on the arrival order of corresponding signal-time curve and a priori knowledge of brain anatomy. The second step is to automatically extract regions of the dominant tissue from each selected IF image by means of an optimal threshold [29].

This paper is organized as follows. The protocol of MR imaging and data recordings is first described followed by the introduction of IFA method and calculation of the pixel-by-pixel rCBV, rCBF, and MTT maps. Computer simulations are presented to validate the application of IFA method on two-dimensional independent factors. Resultant five IF images in conjunction with corresponding signal-time curves from a data set are exhibited. We then calculate and display the rCBV, rCBF, and MTT maps based on the extracted arterial compartment. The averaged ratios for rCBV, rCBF, and MTT between GM and WM from five normal subjects are also computed. Finally, we discuss and conclude this study.

2. Subjects and Data Recording

Five healthy volunteers (three males and two females) aged from 18 to 47 were recruited to participate in this study. Written informed consent was obtained from each volunteer before this study. A multislice gradient-echo EPI pulse sequence on a 1.5-Tesla scanner (Signa CV/i; GE Medical Systems, Milwaukee, WI, USA) was used and the imaging parameters were transaxial imaging, TE/TR  =  60/1000 milliseconds, flip angle  =  90 degrees, FOV = 24 cm × 24 cm, matrix = 128 × 128, slice thickness/gap  =  5/5 mm for 7 slices, one acquisition, and 100 images per slice location. Twenty milliliters of Gd-DTPA-BMA (Omniscan, 0.5 mmol/mL; Nycomed Imaging, Oslo, Norway) followed by 20 mL of normal saline were delivered administratively by a power injector (Spectris; Medrad, Indianola, PA, USA) at a flow rate of 3-4 mL/sec in the antecubital vein. The first thirteen and last thirty-seven images were removed from 100 images and fifty images, which exhibited stable baseline and discernible temporal signal changes, were kept in a slice location for analysis. The temporal resolution is one second. Figure 1 displays part of dynamic perfusion images at an upper slice location encompassing the first circulation from a 39-year-old volunteer. With 16384 (=128 × 128) pixels for each image, the observation of 50 temporal images can be represented by a 50 × 16384 matrix. The signal-to-noise ratios (SNRs) of these five data sets were 51, 55, 77, 83, and 88, respectively.

Figure 1

Perfusion MR images from a 39-year-old volunteer. Images encompass the first circulation (7–30, from left to right, top to bottom) at an upper slice location.

3. Independent Factor Analysis

In this study, we assume that there are N pure tissue types presented on perfusion images and that signal intensity of each pixel is a linear combination of contributions from N pure tissues due to the effect of partial volume mixing in MR scanning. Let the observed noisy mixtures be denoted by a matrix y of size 50 × 16384, the linear mixing by a matrix H of size 50 × N, the independent factors by a matrix x of size N × 16384, and the zero-mean Gaussian noise by a matrix u of size 50 × 16384, where each row in x represents an image for a tissue type, each column of H represents a signal-time curve for an associated tissue type, and covariance matrix of the zero-mean Gaussian noise is denoted by Λ. Under the assumption that the sources x are mutually statistically independent, the IFA-based blind source separation (BSS) technique attempts to recover the unknown mixing matrix H and hidden IF images x from the observed noisy mixtures [27]:

In order to recover the IF images, we assume that each row of the matrix x is a realization of a random variable whose probability density p(x ) is in a form of mixture of Gaussians (MoG) given by

where g stands for the Gaussian function, q is a variable with value running over number of Gaussians with means μ , variances ν used in MoG model for the source x , j = 1,…, N, and ω 's are the mixing proportions satisfying ∑ ω = 1. Let q = (q 1,…, q ) denote all possible combinations of the individual q . The joint probability density p(x) in N-dimensional space can be formed by the product of the one-dimensional MoG's in (2) due to the mutual independence of x , which is itself an MoG

where

In our implementation, n was determined to be 2 for each source x . To estimate the IF model parameters, the Kullback-Leibler distance function [27] is employed to measure the difference between the probability density of the observed signals given the model parameters W, that is, p(y | W), and the observed one, that is, p°(y), which is given by

where W = {H, Λ, } denotes the unknown mixing matrix H, noise covariance Λ, and MoG parameters . The operator E denotes the average over the observed y. The first term in (5) is the negative log-likelihood of the observed signals given the model parameters W, and the second term is the entropy of observed signals, which is independent of W and will henceforth be dropped. Note that, based on p(q, x, y | W) = p(q)p(x | q)p(y | x), p(x | q) = g(x − , V ), and p(y | x) = g(y − H x, Λ), the probability density p(y | W) can be expressed in a closed form:

where p(y | q) = g(y − H , H V H + Λ). By minimizing ε(W) with respect to W based on the expectation-maximization (EM) method and denoting W′ the parameters obtained from the previous iteration, the iterative algorithm can be summarized as follows [27].

Initialize all the unknown parameters, that is, W 0 = {H 0, Λ0, 0}.

Compute the conditional mean of

Update the mixing matrix and noise covariance:

Update the source MoG parameters:

Rescale the parameters to cancel the extra freedom of scaling for facilitating convergence:

Repeat (2)–(5) until the error function ε(W) converges.

Finally, the IF images can be reconstructed from the estimated model parameters W and measured data y using the least-mean-square estimator:

where p(q | y) = p(q)p(y | q)/∑ p(q′)p(y | q′), A = (H Λ−1 H + V −1)−1 H Λ−1, and b = (H Λ−1 H + V −1)−1 V −1 .

4. Calculation of the rCBV, rCBF, and MTT

Prior to the calculation of pixel-by-pixel rCBV, rCBF, and MTT maps, we computed the concentration-time curve C (t) for each pixel using the formula

where k is an unknown constant, TE is the echo time, and S(t) and S 0 are the signal intensities of each pixel at time t and at the baseline, respectively [1, 4, 5, 30]. By using the indicator dilution theory, one can determine the rCBV for each pixel as a ratio of the area integrating over the first pass (e.g., from 1st to 22nd images in Figure 1) of the contrast agent under the concentration-time curve, C (t), to that under the AIF, C (t) [11, 12],

where C (t) is the concentration-time curve for the arterial region. The rCBF can be computed based on the relationship with concentration-time curve for each pixel [31]:

where ⊗ denotes convolution, · denotes multiplication, and R(t) is the residue function for the pixel. The rCBF · R(t) curve for each pixel in (14) can be resolved using the singular value decomposition (SVD) method and the value of rCBF at each pixel was determined by the maximum value of rCBF · R(t) curve [1]. Finally, the MTT of contrast-agent particles passing through a pixel was defined to be [1, 4, 5]

5. Results

In order to validate the implementation of IFA algorithm, a computer simulation was conducted before perfusion data were processed. Since the purpose of this simulation was to verify applicability of the IFA method on the separation of two-dimensional IF images, rather than validate the theory of IFA method (readers are referred to [27] for detailed computer simulations), we simply created four hypothetical IF images (128 × 128) with zero background and foreground related to four major tissue types, that is, artery (IF1: 603 pixels), GM (IF2: 1683 pixels), WM (IF3: 1514 pixels), and others (IF4: 924 pixels including vein, sinus, choroid plexus, and cerebral spinal fluid) (the left four plots in Figure 2(a)). In practice, the four tissue-related areas were segmented from the perfusion data and copied to the corresponding locations. The corresponding averaged temporal profile of each area during the first pass was simplified into five points to generate the hypothetical mixing matrix H (5 × 4) (also see the right most plot in Figure 2(a)) We arranged the four IF images into a matrix x (4 × 16384) which was multiplied by the H to create a matrix y = H x (5 × 16384). Additional multivariate Gaussian random noises with zero mean and diagonal covariance matrix (diagonal terms were 10) were added to the noise-free simulated data y (y is displayed in Figure 2(b)). The SNR was 240.

Figure 2

Computer simulation for the evaluation of the IFA method. (a): The four hypothetical IF images and columns of the hypothetical mixing matrix (the right most plot). (b): The simulated images. (c): The convergence of the IFA algorithm.

Now the task was to estimate four IF images and mixing matrix from the matrix y, to compare them with the hypothetical ones. The number n of Gaussian functions in the MoG model was chosen to be 2 for modeling the probability density function of each IF image and minimizing the computational cost since we have experienced that the use of larger n did not improve the results in either simulated or perfusion data. The values of cost function (the first term in (5)) reduced quickly and converged after around 100 iterations (Figure 2(c)). The correlation value between each pair of hypothetical and estimated IF images and that between each pair of hypothetical and estimated mixing weightings were higher than 0.9999. We further computed a matrix , where was the estimated mixing matrix, and the result was

which was very close to the identity matrix if the estimate was correct. We repeated the simulation when the SNR was reduced to 40, which was lower than that in the normal data. Results showed that the correlation value between each pair of hypothetical and estimated mixing weightings remained as high as 0.9998 and that between each pair of hypothetical and estimated IF image only degraded slight to between 0.9883 and 0.9997. The values of matrix J were

The high correlation values and good approximation of J matrices not only validated the correctness and convergence of our implementation of the IFA method, but also promised the suitability of the IFA method in segregating various tissue types from perfusion MRI data.

The number of IF images (N), that is, number of tissue types, needs to be determined prior to the IFA process. Various numbers, ranging from 4 up to 6, have been assumed in the calculation and we have found that N = 5 elucidated distinctly discernible tissue types, appearing to agree with our knowledge of brain anatomy and physiology. Results from five normal image data sets have been confirmed by a neuroradiologist who is one of the coauthors in this study with expertise in perfusion imaging. One of the results was shown in Figure 3, where the upper panel depicts the resultant five IF images, which are artery (ar), GM, WM, vein and sinus (vs), and choroid plexus (cp), respectively. Figure 3(b) displays the corresponding signal-time curves (columns of H ), respectively, which were all normalized to unit variance and their baselines were shifted to 1.0 for the comparison. The tissue types of IF images were identified based on their anatomical structures and the arrival order of contrast agent, that is, artery follows by GM, WM, vs, or cp. Among all IF images, the IF image with brighter pixels representing artery can be easily recognized because of its signal-time curve presenting the fastest signal drop (red curve in Figure 3(b)). The pixels of each major tissue type can be further segmented out from each IF image using an automatically optimal thresholding technique, that is, the Otsu's method in this study. The final segmentation result was depicted by a color-coded composite image (left panel in Figure 3(c)) and the averaged signal-time curve within each colored area, that is, each tissue type, was computed (right panel in Figure 3(c)).

Figure 3

Classification results of a normal subject (whose perfusion images are shown in Figure 1). (a) Five major tissue types exhibited in IF images are ar, GM, WM, vs, and cp, respectively. (b) The normalized signal-time curves correspond to five IF images, where the one (red) with the fastest signal drop is the artery. (c) Left: a color-coded composite map is used to represent the final segmentation result for different tissue types. Right: five signal-time curves correspond to the averaged intensities of color-coded areas. (d) The rCBV map (left), the rCBF map (middle), and the MTT map (right). (Scale unit for rCBV and rCBF: arbitrary unit; scale unit for MTT: second.)

To create the hemodynamic parametric maps, namely, rCBV, rCBF, and MTT shown in the left, middle, and right panels in Figure 3(d), respectively, we first converted the averaged arterial signal-time curve and temporal intensity profile at each pixel into the concentration-time curves using (12) followed by (13), (14), and (15).

From the segmentation results of the five normal data sets, we have observed that arterial areas were all reliably segmented and the contrast agent consistently arrived first at the artery, followed by GM, WM, vs, and cp. The averaged ratios for rCBV, rCBF, and MTT between GM and WM were 2.139 ± 0.190, 2.598 ± 0.184, and 0.789 ± 0.098, respectively, which were congruent with those in the literature [3, 4, 9, 13].

6. Discussion

This study describes an IFA-based method to classify pixels of the same tissue type on perfusion images based on bolus transit profiles and the assumptions of spatial independence as well as the partial-volume mixing effect. The IFA method is flexible in learning the source densities from observed data so that sources can be more accurately modeled by the mixture of Gaussians for the facilitation of subsequent separation. The IFA technique is related to projection pursuit [32, 33], where “interesting” projections of multidimensional data are pursued for optimal visualization of data and exploratory data analysis. Projection pursuit is usually performed by finding directions in which the data is least Gaussian distributed. Since the independent factors in the IFA are modeled by mixture of Gaussian functions, it is interesting to investigate whether the independent factors present maximal non-Gaussian clustering structures in future work. Our results indicated that the brighter pixels in a cluster were homogeneous in most tissues and can be segmented out using Otsu's method for automatic determination of optimal thresholds [29]. Otsu's method is simple in terms of implementation and computation and is robust to histogram irregularities caused by noise.

The determination of the number N of IF images is an important issue to be addressed. As suggested by Attias [27], one can determine this number using a comprehensive method, for example, the model comparison method [34], or a simpler but less precise method, for example, the number of significant eigenvalues of data covariance matrix. In this study, we decided N based on a priori knowledge of the brain anatomy and arrival order of blood flows for different tissues. When N was chosen to be five, each IF image exhibited one dominating tissue cluster with brighter intensities, allowing consistent and reliable segregation of the artery, GM, WM, vs, and cp from the normal subjects. When N was less than five, two major tissue types, for examples, artery with GM or vs with cp, can appear at one of IF images. On the other hand, if N was larger than five, either one major tissue type was split into two IF images, for example, part of sinus was separated from vs, or repeatedly exhibited at two IF images.

It is worthy to note that the number of unknown parameters depends on the number of dynamic MR images. The more images were used, the more parameters needed to be estimated in the noise covariance Λ and in the mixing matrix H in which columns corresponded to the signal-time curves. We have found that the data points encompassed the duration of the first and second circulations, for example, 50 in the illustrative example (Figure 3), were adequate to resolve the IF images. A simple way to determine data length is to average each perfusion image and plot the averaged images with respect to time. The curve delineating the first and second passes is readily seen and the duration can be easily decided.

7. Conclusions

The IFA-based method provides several advantages for cerebral hemodynamic studies. First, various tissue compartments on perfusion images can be classified systematically. Second, the arterial compartment can be modeled consistently on the same slice location for calculation of rCBV, rCBF, and MTT maps. Third, bolus transit profiles of these tissues can be well separated, providing information in addition to rCBV, rCBF, and MTT maps, such as temporal scenarios and recirculation of contrast agent. This method also promises differentiation of pathological and nonpathological blood-supply patterns in future clinical applications.