Hemodynamic Segmentation of Brain Perfusion Images with Delay and Dispersion Effects Using an Expectation-Maximization Algorithm

Automatic identification of various perfusion compartments from dynamic susceptibility contrast magnetic resonance brain images can assist in clinical diagnosis and treatment of cerebrovascular diseases. The principle of segmentation methods was based on the clustering of bolus transit-time profiles to discern areas of different tissues. However, the cerebrovascular diseases may result in a delayed and dispersed local perfusion and therefore alter the hemodynamic signal profiles. Assessing the accuracy of the segmentation technique under delayed/dispersed circumstance is critical to accurately evaluate the severity of the vascular disease. In this study, we improved the segmentation method of expectation-maximization algorithm by using the results of hierarchical clustering on whitened perfusion data as initial parameters for a mixture of multivariate Gaussians model. In addition, Monte Carlo simulations were conducted to evaluate the performance of proposed method under different levels of delay, dispersion, and noise of signal profiles in tissue segmentation. The proposed method was used to classify brain tissue types using perfusion data from five normal participants, a patient with unilateral stenosis of the internal carotid artery, and a patient with moyamoya disease. Our results showed that the normal, delayed or dispersed hemodynamics can be well differentiated for patients, and therefore the local arterial input function for impaired tissues can be recognized to minimize the error when estimating the cerebral blood flow. Furthermore, the tissue in the risk of infarct and the tissue with or without the complementary blood supply from the communicating arteries can be identified.


Introduction
The advent of modern imaging modalities has enabled studies to unravel temporal hemodynamic patterns in different regions of the brain to assist in the assessment of cerebrovascular diseases. Dynamic-susceptibility-contrast MR (DSC-MR) imaging records signal changes associated with different blood supply patterns following the intravenous injection of a bolus of contrast agent. Based on the bolus profile in the arterial compartment and on indicator dilution theory [1,2], computations of cerebral hemodynamic parameters, such as relative cerebral blood volume (rCBV), relative cerebral blood flow (rCBF), mean transit time (MTT), and time to peak (TTP), are possible. Hemodynamic parameter maps have been intensively used in clinical applications, including the assessment of brain tumors [3][4][5], brain ischemia [4][5][6][7], occlusive cerebrovascular disease [8], and radiation necrosis [3,5].
Classification of brain hemodynamics into different components is essential to facilitate the analysis and assessment of brain perfusion. In our previous studies, fast independent component analysis (FastICA) [9,10], the noiseless independent factor analysis (NIFA) [11], and analysis based on expectation-maximization (EM) of a mixture of multivariate Gaussians (MoMG) [12] were developed to discern brain tissues. In the EM-MoMG method, the data was similarly zero-mean normalized and reduced by principal component analysis (PCA) [13]. Under the assumptions that the distribution of the reduced data associated with each tissue type was multivariate Gaussian distributed and that the overall distribution of the reduced data was a MoMG [14], alternating iterative expectation (E) and maximization (M) steps [15] were performed. The E-steps calculated the expected log-likelihood conditions using the observed data along with the currently estimated model parameters. The M-steps updated the model parameters by maximizing the log-likelihood values. During this EM process, the posterior probabilities in relation to each tissue type were obtained at each pixel, which in turn provided sufficient statistical validity, i.e., maximal probability, to determine the tissue type of each pixel. As a consequence, all tissue compartments could be segmented out simultaneously using the EM-MoMG method, instead of the one-by-one segmentation process used in the FastICA method. However, the EM-MoMG method required a good initial guess of the model parameters to achieve satisfactory results. In this study, we improved the EM method by using the results of hierarchical clustering (HC) on whitened perfusion data as initial parameters for a MoMG.
The segmentation of brain tissues using perfusion images relies on the clustering of bolus transit-time profiles to discern areas of different tissue types. However, cerebrovascular diseases would change the local perfusion of the impaired tissues with different levels of delay and dispersion, and therefore complicate the overall distribution of the signal profiles. The delay and dispersion phenomena were resulted from the defect of the blood vessel structures, such as the internal carotid artery (ICA) stenosis, moyamoya disease, and arteriovenous malformation. Assessing the accuracy of the segmentation technique under delay/dispersed circumstance is critical to accurately evaluate the severity of the vascular disease and locate the impaired tissues. However, previous literatures have focused on the effects of delay and dispersion in the deconvolution algorithm for estimating the rCBF [16][17][18], but not the effects in the tissue segmentation and impaired tissue recognition. In this study, Monte Carlo simulations were conducted to evaluate the performance of proposed method under different levels of delay, dispersion, and noise of signal profiles in tissue segmentation. Finally, we analyzed the perfusion data from five normal participants, a patient with unilateral ICA stenosis, and a patient with moyamoya disease using the proposed method.

Materials and Methods
Tissue Segmentation using the HC-EM-MoMG Method Brain regions were extracted from perfusion images using Otsu's method [19], followed by an erosion and dilation operation [20]. The brain region for each image was assumed to comprise N pixels, and the observation of 65 temporal images was represented by a 65|N matrix. The dimension of each data set was reduced by PCA [13] from 65|N to r|N. During PCA, at least 99% of the data variance was retained.
Similar to PCA, the data whitening process first calculates eigenvalues and corresponding eigenvectors for the covariance matrix of the zero-mean normalized data. The zero-mean data were then transformed via a whitening matrix, which was constructed of eigenvectors and eigenvalues, so that the covariance matrix became an identity matrix [13]. The compressed, transformed data resulting from whitening dimension reductions were subsequently used in HC.
During HC, the compressed data are denoted by a matrix X with size r|N and each r|1 column vector in X is referred to as a feature vector, x. Initially, HC was carried out in the data set X based on a created cluster tree with a multilevel hierarchy in which any two nearby clusters at one level become merged into one cluster at the next higher level. The HC algorithm comprised three major parts. First, each group contains a single feature vector of X, i.e., N groups in X, and the element in the dissimilarity matrix is the squared Euclidean distance between any two featured column vectors in X, which was defined by where d 2 rs is the Euclidean distance between two featured column vectors x r and x s . Second, in order to construct the cluster tree, we employed the Ward's method [21] to measure the distance between two clusters which was defined by where d 2 (p,q) is the distance between two clusters of featured vectors, n p and n q are the numbers of featured vectors in cluster p and q, respectively, kk 2 is the Euclidean distance, x x p and x x q are the centroids of clusters p and q, respectively, which are defined by where x pi is the ith featured vector in cluster p and x qj is the jth objects in cluster q. The two clusters with the smallest betweengroup distance are grouped together to form a new group. The algorithm proceeds until all of the feature vectors fall within a single group, thus forming a hierarchical clustering tree. Third, the level or scale of clustering is determined by cutting the hierarchical cluster tree. The HC results for the condensed data, X, were used to be the initial guess of the mean vectors, the covariance matrices of the MoMG, and the proportion of each tissue class for the subsequent EM estimation [12]. The E-step of EM estimation is to compute the posterior probabilities of tissue classes p(i|x n ,h j21 ) based on the estimation from previous j-1 th iteration: Where i = 1,…,K represent labels of tissue classes, g xn m i ,S i ½ represents the multivariate Gaussian density with the mean vector m i and covariance matrix S i of tissue class i, p i denotes the proportion of each class i, and h j{1 denotes the parameters fm j{1 i ,S j{1 i ,p j{1 i g at the j-1 th iteration. The M-step consists of estimating parameters h j as follows (please see Appendix in [12] for details): Once all tissues of interest were identified, averaged signal-time curves for each tissue type were computed and the arterial input function (AIF) was modeled from the averaged concentration-time curve of artery in order to compute the rCBV, rCBF, and MTT of each segmented tissue types [22][23][24]. Specifically, we first computed the concentration-time curve C t (t) for each pixel using the formula: where k is a 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. 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 of the contrast agent under the concentration-time curve, C t (t), to that under the AIF, C a (t), rCBV~Ð first pass The rCBF can be computed based on the relationship with concentration-time curve for each pixel: where 6 denotes convolution, ? denotes multiplication, and R(t) is the residue function for the pixel. The rCBF : R(t) curve for each pixel 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 [25]. Finally, the MTT of contrast-agent particles passing through a pixel was defined as

Design of Hypothetical Compartments for Monte Carlo Simulations
Simulated dynamic images based on a region of interest (ROI) in one of the raw data sets were adopted. In those images each simulated data set accommodated the possibilities of having four to nine hypothetical clusters (Fig. 1). Pixels within each hypothetical tissue area were assumed to be spatially independent, and pixel intensities across time were assumed to be multivariate Gaussiandistributed. The averaged signal-time curves, each of which was a 65|1 vector, and the covariance matrices of intensities across time, each of which had a dimension of 65|65, within ROIs in the raw data were employed as hypothetical parameters for the MoMG model to create sets of noise-free simulated dynamic images using random number generators. In addition, Gaussian noise was added to each set of noise-free simulated dynamic images to produce the signal to noise ratio (SNR) levels of 40 and 70. For each of the hypothetical clusters and each noise level, a Monte Carlo simulation comprising 1000 runs was performed. Accordingly, the total number of simulated image sets was 6 (number of hypothetical clusters) |3 (noise levels) |1000(repe-(repetitions for each combination) = 18000. Each set of simulated dynamic images was rearranged into a 65|N matrix before classification, where N (N~5502 in the simulation) was the number of pixels representing the brain area.

Monte Carlo Simulations for EM-MoMG with Differenet Initializations
The parameters of MoMG for EM segmentation were initialized by 2 different means: 1) HC results from whitened data; 2) random sampling of columns from the matrix representing the PCA data, which was used in our previous EM-MoMG method [12]. The number of tissue classes to be tested ranged from 4 to 9. After comparing posterior probability values at each pixel, the cluster resulting in the maximum probability was identified, that pixel was then labeled accordingly. The classification rate for each cluster, denoted by r i , was defined by the observed proportion of agreements, i.e., r i~fi =n i , where f i was the number of agreements for the i th cluster and n i was the total number of pixels in the i th class. The classification rate for a set of simulated dynamic images, denoted by R, was computed by R~100| , where K is the number of clusters. No delay and dispersion effects were added to these simulated perfusion images.

Monte Carlo Simulations with Delay and Dispersion
Delay and dispersion are two popular phenomena in the hemodynamics of subjects with cerebrovascular diseases. Thus, Monte Carlo simulations were conducted to investigate the influences of these two factors on the proposed HC-EM-MoMG method. The concentration-time curves, C tissue (t), for anomalous tissue types were created using the following formula [22] where the effective residue function R eff ,tissue (t) [18] is given by in which t delay,tissue and b tissue denote time delays relative to AIF and dispersion, respectively. The values of MTTtissue, rCBF tissue , t delay,tissue , and b tissue were estimated to simulate the concentrationtime curves of the delayed/dispersed artery (dArtery), delayed/ dispersed gray matter (dGM), and delayed/dispersed white matter (dWM). Specifically, we calculated the values of rCBF and MTT for artery, GM, and WM, based on the raw data of hypothetical compartments. The normal rCBF and MTT values of artery, GM, and WM were also used for dArtery, dGM and dWM, respectively. The derived simulated concentration-time curves were then converted into signal-time curves using the following equation where S 0 is the baseline signal and TE is echo time. In addition to the signal-time curves for dArtery, dGM, and dWM, signal-time curves for artery, GM, WM, vein combined with sinus (vein+sinus), sinus, as well as for CSF including the corpus plexus (CSF+cp) tissues were also used in these simulations. Hence, there were 9 hypothetical tissue types in the test of delay and dispersion effects.
With t delay values of 0, 1, 2, 3, 4, and 5 s and dispersion b values of 0, 1, 2, 3, 4, and 5 s, overall 36 combinations of abnormal conditions were simulated (Fig. 2). The number of voxels in the impaired components, i.e., dArtery, dGM and dWM, were designated to be 50%, 25%, or 13% of the total voxels in the artery, GM and WM components. Noise-free and SNR levels of 40 and 70 were created. The total number of simulated data sets was 6 (different t delay 's) |6 (different b's) |3 (different numbers of voxels in abnormal components) |3 (noise levels) |1000 (repetitions for each combination) = 324000. Fig. 3 illustrates the 9 simulated mean signal-time curves in which the time settings of t delay and b for the impaired compartments were both 5 s.

Participants and Data Recording
This study received prior approval from the Institutional Review Board of Taipei Veterans General Hospital. Each participant provided written informed consent before participating in this study. Five healthy volunteers (3 males and 2 females) aged from 18 to 47 years, one 78 year-old male with a unilateral ICA stenosis, and one 57 year-old female with moyamoya disease participated in this study. A multi-slice gradient-echo echo-planarimaging pulse sequence on a 1.5 Tesla scanner (Signa CV/i, GE Medical Systems, Milwaukee, WI, USA) was used to acquire dynamic perfusion images. For the healthy subjects, trans-axial imaging was used with TE/TR = 60/1000 ms, flip angle = 90 degrees, FOV = 24 cm 624 cm, matrix = 1286128, slice thickness/gap = 5/5 mm for 7 slices, one acquisition, and 70 images per slice location with a one second temporal resolution. Twenty ml of Gd-DTPA-BMA (OmniscanH, 0.5 mmol/ml, Nycomed Imaging, Oslo, Norway) followed by 20 ml of normal saline were delivered administratively using a power injector (SpectrisH, Medrad, Indianola, PA, USA) at a flow rate of 3-4 ml/s in the antecubital vein. Some imaging protocol settings were changed for the ICA patient, i.e., TE/TR = 40/1000 ms, flip angle = 60 degrees, slice thickness/gap = 7/14 mm for 7 slices. All routines were implemented using MATLAB (MathWorks, Inc., Natick, MA, USA) codes and carried out on a 2.66 GHz Intel-Core2-based personal computer. Table 1 presents the average classification rates+standard deviations (in percentages) for Monte Carlo simulations obtained from the EM-MoMG method initialized by 2 different approaches. The EM-MoMG method initialized by the results from HC on whitened data outperformed the initialization approach using the random sampling on PCA data. The average classification rates were higher using HC on whitened data (between 99.0+0.39% for noise-free conditions with four hypothetical tissue clusters and 93.6+1.42% for the condition of SNR = 40 with 9 hypothetical tissue clusters, Table 1). The standard deviations of classification rates were also smaller (always less than 1.42%) using the HC results on whitened data as initialization, thereby implying the HC-EM-MoMG method can produce accurate and reliable segmentations.

Simulation Results for Different Initializations of EM-MoMG
The mean classification rates from all initialization methods decreased slightly with the change from the noise-free condition to SNR = 40. In the case of the proposed HC-EM-MoMG method, the minimal mean classification rate and maximal standard deviation were 93.6% and 1.42%, respectively, in the 9 hypothetical tissue clusters with SNR = 40 ( Table 1). The classification performances degraded as the number of hypothetical tissue clusters increased. In the worst case, i.e., SNR = 40, the results obtained from the HC-EM-MoMG method were from 98.9+0.40% to 93.6+1.42%.
In addition, we computed minimal description length (MDL) [26] values when different numbers of tissue classes were used in  Simulation Results with Various Delays, Dispersions, and Percentages of Abnormal Voxels Table 2 shows the simulation results with different delays and dispersions under various noise conditions. The results showed that the classification rates increased as either delay or dispersion increased. For example, the classification rates under noise-free conditions increased from 71.0+5.32% to 89.9+5.92% when the delay increased from 1 s to 5 s, and increased from 76.4+5.87% to 87.3+4.36% as dispersion increased from 1 s to 5 s. In general, the classification rates were satisfactory when the differences in either delay or dispersion between normal and corresponding impaired compartments were longer than one second. As examples, the classification rates were mostly over 80% when (b, t delay ) = (0, 2), (1, 1), (2, 0), and were over 85% when (b, t delay ) = (0, 3), (1, 2), (2, 1), (3,0). The results also show that as noise increased from the noise-free condition to SNR = 40, the mean classification rates also decreased. As examples, the mean classification rates decreased from 87.0+7.94% with noise-free condition to 86.6+6.55% with SNR = 40 at (b, t delay ) = (5,5). In addition, the mean classification rates obtained from the six normal tissues and the three impaired compartments (dArtery, dGM, dWM) under different levels of delay and dispersion were 4% to 20% lower than the rates from the 9 normal tissues shown in Table 1. This result suggested that some of the perfusion profiles of the impaired compartments may be equivocal and more difficult to distinguish. However, as delay and dispersion were both longer than one second, the mean classification rate was greater than 85%.
The classification rates that were determined using different percentages of abnormal voxels are presented in Table 3. The results demonstrated that a higher abnormal proportion produces a lower mean classification rate and a larger standard deviation. The mean classification rates decreased from a range between 78.2 and 90.1% with a 13% abnormal percentage to a range between 56.3 and 88.9% with a 50% abnormal percentage. Moreover, the standard deviations increased from at most 6.11% with a 13% abnormal percentage to 9.65% with a 50% abnormal percentage.

Classification Results from Normal Participants
Data from the five normal participants were processed and their perfusion images segmented into different compartments based on minimizing MDL values. Various tissue types, such as artery, GM, WM, vein+sinus, sinus, CSF+cp, CSF, vein with noise, sinus with noise, as well as artifacts, were segmented from normal data sets. The ratios of GM to WM for rCBV, rCBF, and MTT were 2.19660.097, 2.25960.119, and 0.96860.023, respectively. These ratios were in agreement with other published reports [3,27,28].

Results from a Patient with Right Internal Carotid Artery (ICA) Stenosis
The stenosis subject was a 78 year-old male subject with a rightside ICA 99% stenosis. A right-lateral neck angiogram of the subject showed the high-degree stenosis on the right ICA (white arrows in Fig. 4(b)). Another cerebral angiogram with an anteriorposterior projection at the early arterial phase, and with contrast injection from the aortic arch, was provided to illustrate the delayed perfusion on the right ICA. The right hemisphere of the brain exhibited delayed circulation (white arrow in Fig. 4(c)). This delayed perfusion was observed on consecutive trans-axial DSC-MR images (see Fig. 4(a)). The subject's right hemisphere (abnormal side; left side of the images) had a signal drop and recovery between the 14 th and 28 th images, whereas the left hemisphere (normal side; right side of images) had a signal drop and recovery between the 12 th and 22 nd images, i.e., the signal drop occurred 2 s earlier and recovery occurred 6 s earlier on the normal side than on the stenosis-affected side. In addition, a trans-axial T 2 -weighted MR image ( Fig. 4(d)) and a diffusion-weighted image (Fig. 4(e)) at a similar slice location as the DSC-MR image confirmed that there was no evidence of infarct or stroke on the slice. Segmented results (Fig. 4(f)) from the DSC-MR image using HC-EM-MoMG shows nine compartments, namely artery (red), GM (green), WM (brown), CSF+cp (gray), right ICA stenosis induced delayed/dispersed artery (pink), delayed/dispersed GM (yellow), delayed/dispersed WM (light brown), vein+sinus (blue), and artifact (dark blue). In Fig. 4 (f), it is worth noting that the GM and WM component maps cross the midline from the normal side to the stenosis side, especially in the anterior and posterior regions. Those results support the clinical findings that the subject's anterior and posterior communicating arteries in the circle of Willis were intact. The averaged signal-time curves and corresponding spatial maps ( Fig. 4(g) and (f)) demonstrate that tissues with similar characteristic perfusion profiles can be grouped by the EM algorithm. The results suggest that the impaired parts (dArtery, dGM and dWM) can be discerned from their contralateral counterparts (artery, GM and WM) by applying the HC-EM-MoMG method. One should note that the delayed or dispersed arterial input function (dAIF) obtained from the signal-time curve of dArtery was used to calculate the rCBF for dGM and dWM for minimizing the error induced by the delay and dispersion whereas the AIF obtained from the normal artery was used for other normal tissues such as GM, WM, CSF+cp, and vein+sinus. Fig. 4(h)-(j) exhibit hemodynamic-parametric maps (rCBV, rCBF and MTT) and the details for each of the classified tissue types are given in

Results from a Patient with Moyamoya Disease
Another case was a 57 year-old female patient with moyamoya disease. She suffered from sudden vision loss on her left side for 3 days. Her right cerebral angiogram at the early arterial phase with lateral projection showed stenosis of the distal internal carotid artery, proximal segments of the anterior and middle cerebral arteries (red arrows in Fig. 5(b)). Intracranial collateral arteries, socalled moyamoya vessels, are seen in the basal ganglia region (blue arrow in Fig. 5(b)). The prominent branches of external carotid artery (arrowheads in Fig. 5(b)) form the extracranial collaterals. Another arterial phase angiography taken at 1.5 seconds later than The abnormal percentage was the ratio of number of abnormal voxels to total number of voxels in a tissue type.
In the case of abnormal percentage = 50%: numbers of voxels in dArtery, dGM and dWM are 551*50%&275, 1741*50%&870 and 1636*50% = 818, respectively. Other two cases can be calculated in the same manner by replacing percentage of 50% with 25% and 13%, respectively. doi:10.1371/journal.pone.0068986.t003 Fig. 5(b) shows the late arrival blood supply (arrows in Fig. 5(c)) via the collaterals. Her right parietal region is less blood-irrigated (arrowheads in Fig. 5(c)). The full perfusion could be shown by the consecutive trans-axial DSC-MRIs (see Fig. 5(a)). In addition, a trans-axial T2-weighted MR image shows abnormal high signals in her parietal-occipital region (arrows in Fig. 5(d)) and a diffusionweighted image confirms an acute infarct with diffusion restriction (arrows in Fig. 5(e)). Segmented results (Fig. 5(f)) from DSC-MRI shows eight hemodynamic components, namely artery (red), gray matter (green), white matter (brown), CSF (gray), ischemic or infarct area (purple), areas supplied by collateral arteries with risk of infarct (cyan), vein and sinus (blue), and artifacts (dark blue). The averaged signal-time curves (Figs. 5(g)) demonstrated that the artery arrived earliest with maximum signal drop followed by the GM or vein+sinus, WM or CSF, risk of infarct, and infarct. Infarct area, mainly the right parietal-occipital brain, was the last to receive blood flow and shows much smaller signal drops compared to other areas. The resulting right occipital infarct as seen in Figs. 5(d) and (e) may be explained by the hemodynamic aberration shown in the infarct area (purple). Whereas the risk of infarct area presented the averaged signal-time curve with delay and dispersion effects compared to that of the GM tissue, and the spatial map of the risk of infarct resided near to the infarct area. In addition, using the classified regions ( Fig. 5(f)) as masks on parametric images (Figs. 5(h)-(j)), the hemodynamic parameters, (rCBV, rCBF, MTT), for compartments artery, GM, WM, CSF, infarct, risk of infarct, and vein+sinus, can be easily calculated, respectively (see Table 5). The infarct area presented much lower rCBV (1.46+0.93 ml?100 g 21 ) and rCBF (9.95+5.37 ml?100 g 21 ?min 21 ) but higher MTT (8.50+2.63 s) in comparison with that of GM and WM areas. In contrast, the risk of infarct area performs relatively toward normal rCBV (4.49+1.83 ml?100 g 21 ) and rCBF (29.986 12.68 ml?100 g 21 ?min 21 ), which implies this area is worth in therapy to recover. These values can be re-computed after treatment to assess therapeutic effects.

Discussion
In this study, we improved the EM method by using the results of HC on whitened data to initialize the model parameters used in the MoMG model. It is worth noting that performing HC on whitened data is equivalent to conducting clustering on the independent component images resulting from FastICA [10]. To clarify, let the zero-mean-normalized and whitened data be X. In a matrix of such data, each row represents an image and each column encodes the temporal information for each voxel, which can be regarded as a feature vector in the clustering process. In addition, let the optimal rotation matrix in the FastICA method be denoted by G. The resultant independent component images that manifest one major tissue type in each image are therefore denoted by GX. Since any rotation matrix holds the property G T G~I, where the superscript 'T' stands for the matrix transpose, it follows that the Euclidean distance between any pair of feature vectors, namely, x r and x s , remains the same after rotational transformation, i.e., (Gx r {Gx s ) T . This shows that clustering on x r 's is the same a clustering on Gx r 's.
Monte Carlo simulations were carried out to compare the performance using HC results from whitened data and the random sampling of columns from the matrix representing the PCA data. The results of HC on the whitened data produced superior data clustering in terms of accuracy and variance than that obtained from random sampling on PCA data. That superiority was due to the similarity of scales in the corresponding components in any pair of column (feature) vectors from the whitened data. Through that similarity, each component in the subtraction of each paired vectors contributed equally to the calculation of Euclidean distance. On the contrary, in PCA data, the first components were the dominant terms in any pair of column vectors, and they overwhelm the contributions from the remaining components, producing results in which the features from the rest of components of each column vector have been suppressed. Thus, clusters resolved from the dissimilarity matrix can be better discriminated based on the whitened data. In brief, the simulation results suggest that the HC-EM-MoMG method performed well, particularly when the SNR was higher than 40.
Various tissue types, such as artery, GM, WM, vein+sinus, sinus, CSF+cp, CSF, vein with noise, sinus with noise, as well as artifacts, were segmented from normal data sets. Additionally, impaired tissue types, namely, dArtery, dGM, dWM, were segmented from an ICA subject's data, and the infarct and risk of infarct were distinguished from others normal tissues in the moyamoya data. Since the slice locations from each of the normal subjects were not exactly the same, the segmented tissue types and the optimal number of clusters (between 7 and 8) varied slightly  among subjects. Nevertheless, the artery, GM, WM, CSF+cp (or CSF) could be consistently segmented from the five normal data sets and the delivery of the contrast agent appeared in the following order: artery, GM, WM, and CSF+cp (or CSF). The proposed HC-EM-MoMG method successfully showed spatial and temporal hemodynamic patterns for each of the dissected tissue compartments. Three advantages of using the proposed method are as follows. First, the averaged signal-time curve for arterial tissue can be used as an AIF; in particular, the dAIF can be identified in the subject with an ICA stenosis. That allows calculation of the hemodynamic parameters of abnormal tissues, which can be used to compensate for delay and dispersion effects. In the ICA stenosis subject, for example, the rCBF values of the dGM and dWM obtained from the normal AIF were 46.21+17.75 and 22.90+8.06, respectively, whereas from the dAIF, the values were 39.51+15.32 and 19.54+6.60, respectively. Second, in this method, voxels with similarly abnormal perfusion can be grouped into the same cluster, the corresponding signal-time curves and the hemodynamic parameters, namely rCBV, rCBF, MTT, and TTP, can be quantified for each tissue type, which can aid in diagnosis and therapeutic assessment. As presented in Fig. 4(f), the spatial location of each impaired compartment (dArtery, dGM, dWM) can be identified, and theirs hemodynamic parameters showed longer TTP and MTT, but lower rCBF compared to the normal tissues (see Table 4). In the case of moyamoya, the infarct area exhibited smaller signal drop and higher mean signal intensity of averaged signal-time curve which was similar to that of CSF, implying insufficient blood supply and death of brain tissue ( Fig. 5(g)). This also can be observed with low rCBV and rCBF in infarct area. However, the hemodynamic parameters of risk of infarct, which were better than that of infarct area, presented lower rCBV, rCBF but longer TTP and MTT compared to GM. Although the risk of infarct exhibited the ''toward-abnormal'' states, it was potentially salvageable by clinical intervention for preventing enlargement of the infraction. The value of the recognition of risk of infarct area using HC-EM-MoMG method may facilitate the prognosis for treatment, which is similar to using MR perfusion-diffusion mismatch in identifying patients with acute ischaemic stroke for thrombolysis [29,30]. Third, the segmented spatial maps derived from the method are useful when evaluating the hemodynamic compensation mechanism from collateral circulation through the circles of Willis; in particular, the integrity of the anterior communicating artery (AcoA), both side posterior communicating arteries (PcoA), and the basilar artery (BA) can be assessed. This compensation mechanism is illustrated in Fig. 4 (f) where the spatial distribution of the normal GM crosses the midline of two hemispheres within the anterior carotid artery (ACA) and posterior carotid artery (PCA) territories. As a result, the hemodynamics of these two areas on the stenotic side and in the normal GMs were grouped into one cluster, thereby suggesting the presence of perfusion compensation via the AcoA and BA. This result was supported by an angiogram assessment (red arrows in Fig.4 (c)) in which AcoA and BA patency was evident.
In future studies, the concept of the bargaining problem can be implemented for the optimisation of the parameters to assess the delay and dispersion effects. From the perspective of Game theory [31], the delay and dispersion parameters can be served as two players in the game, and the optimisation of these parameters can be achieved by the notion of an equilibrium point [32][33][34]. The optimisation of this two-player benefit in a non-zero-sum game is based on the optimisation of payoffs in the terms of game theory, where we expect to observe equilibrium in a set of alternatives in a spatial game theoretical framework [35]. Target recognition techniques, such as the cross-plot method [36,37], can also be integrated for the automatic identification of tisse type after hemodynamic clustering using our proposed HC-EM-MoMG algorithm. Cross-plots of binary patterns (i.e., the binary maps of tissue templates) are exploered as image signatures for the observed target (i.e., the binary maps of tissue clusters) by capture of the spatial features with minimal computation [36,37]. Some other segmentation mehods can be applied to facilitate the preprocessing of brain perfusion images, for example, an active contour model can automatically define the region of interest to remove the nontarget image pixels, such as skull and scalp, or to locate the brain lesion, such as brain tumors [38].
In summary, we improved the EM method by using the results of HC on whitened data to initialize the parameters of MoMG models, termed as HC-EM-MoMG method. The results of Monte Carlo simulations confirmed that the mean classification rates of the Figure 5. Segmentation results using the HC-EM-MoMG method for a patient with moyamoya disease. (a) Consecutive trans-axial DSC-MRI (15 th to 38 th second). (b) Right cerebral angiography at early arterial phase on lateral projection shows stenosis of the distal internal carotid artery, proximal segments of the anterior and middle cerebral arteries (red arrows). The moyamoya vessels are seen in the basal ganglia region (blue yellow). The prominent branches of external carotid artery (arrowheads) form the extracranial collaterals. (c) Another arterial phase angiogram taken at 1.5 seconds after (a) shows the late arrival blood supply (arrows). Her right parietal region is less blood-irrigated (arrowheads). (d) A trans-axial T2weighted MR image shows abnormal high signals in her parietal-occipital region (arrows). (e) Diffusion-weighted image confirms an acute infarct with diffusion restriction (arrows). (f) Segmented results shows eight hemodynamic components, namely artery (red), gray matter (green), white matter (brown), CSF (gray), ischemic or infarct area (purple), areas supplied by collateral arteries with risk of infarct (cyan), vein and sinus (blue), and artifacts (dark blue). (g) The averaged signal-time curves of corresponding tissues demonstrated that the artery arrived earliest with maximum T2* signal drop followed by the GM or vein+sinus, WM or CSF, risk of infarct, and infarct. (h) rCBV map (scale unit is ml?100 g 21 ). (i) rCBF map (scale unit is ml?100 g 21 ?min 21 ). (j) MTT map (scale unit is second). doi:10.1371/journal.pone.0068986.g005 Table 5. The hemodynamic parameters of segmented areas for the patient with moyamoya disease.

Artery
GM WM CSF infarct risk of infarct vein+sinus proposed method can achieve over 93.6% and small standard deviations, which was superior to that using the random sampling as initializations. Moreover, this is the first study to assess the delay and dispersion effects on the hemodynamic segmentation, which is critical to the cerebrovascular diseases. Our results of Monte Carlo simulations showed that most classification rates of the HC-EM-MoMG method were greater than 80% whenever the differences in delay or dispersion between normal and corresponding impaired compartments were longer than one second. Finally, the analysis of data from the patients with unilateral ICA stenosis and moyamoya disease illustrated the effectiveness of the method for identifying impaired compartments. The proposed method can quantify impaired hemodynamics and serve as an aid in diagnosis and therapeutic assessment, such as evaluation of the integrity of the circle of Willis and identification of the risk of infarct area.