# Partial volume correction for arterial spin labeling using the inherent perfusion information of multiple measurements

## Abstract

### Background

Arterial spin labeling (ASL) provides a noninvasive way to measure cerebral blood flow (CBF). The CBF estimation from ASL is heavily contaminated by noise and the partial volume (PV) effect. The multiple measurements of perfusion signals in the ASL sequence are generally acquired and were averaged to suppress the noise. To correct the PV effect, several methods were proposed, but they were all performed directly on the averaged image, thereby ignoring the inherent perfusion information of mixed tissues that are embedded in multiple measurements. The aim of the present study is to correct the PV effect of ASL sequence using the inherent perfusion information in the multiple measurements.

### Methods

In this study, we first proposed a statistical perfusion model of mixed tissues based on the distribution of multiple measurements. Based on the tissue mixture that was obtained from the high-resolution structural image, a structure-based expectation maximization (sEM) scheme was developed to estimate the perfusion contributions of different tissues in a mixed voxel from its multiple measurements. Finally, the performance of the proposed method was evaluated using both computer simulations and in vivo data.

### Results

Compared to the widely used linear regression (LR) method, the proposed sEM-based method performs better on edge preservation, noise suppression, and lesion detection, and demonstrates a potential to estimate the CBF within a shorter scanning time. For in vivo data, the corrected CBF values of gray matter (GM) were independent of the GM probability, thereby indicating the effectiveness of the sEM-based method for the PV correction of the ASL sequence.

### Conclusions

This study validates the proposed sEM scheme for the statistical perfusion model of mixed tissues and demonstrates the effectiveness of using inherent perfusion information in the multiple measurements for PV correction of the ASL sequence.

## Background

The arterial spin labeling (ASL) sequence provides a noninvasive way to measure the cerebral blood flow (CBF) by utilizing the magnetically labeled arterial blood water as an endogenous tracer to create a “label” image  and subsequently repeats the process to create a “control” image without labeling the arterial blood. The subtraction of the label and control images becomes the perfusion image, which reflects the amount of the arterial blood that is delivered to each voxel after the transit time . Due to fast scan of the label and control images, the perfusion image (label/control difference) is very noisy; therefore, multiple label/control image pairs are commonly acquired and are averaged to impress the noise.

The spatial resolution of the ASL sequence was approximately 3–6 mm. The CBF estimation was contaminated by the partial volume (PV) effect, which results in less accuracy of the CBF quantification . For accurate PV correction, the perfusion contributions of different tissues inside a mixed voxel should be estimated separately . Asllani et al.  proposed a linear regression (LR) method, in which the CBF values of both gray matter (GM) and white matter (WM) are assumed to be constant within an n × n × 1 regression kernel. Under this assumption, the LR method can generate the separate GM’s and WM’s CBF maps, but spatial smoothing may be also introduced into the CBF maps. Then, several methods have been proposed to alleviate the smoothing effect of the LR method [6,7,8]. For multiple inversion-time (TI) ASL data, Chappell et al. reported a PV correction method using a spatially regularized kinetic curve model . To the best of our knowledge, all of the current PV correction methods were performed directly on the averaged image of the multiple label/control pairs, thereby ignoring the inherent perfusion information of the mixed tissues that are embedded in the multiple measurements.

The aim of the present study is to correct the PV effect of the ASL sequence by using the inherent perfusion information of multiple measurements. It was reported that the multiple measurements of the ASL sequence could be regarded as noisy realizations of the original distribution . Therefore, for each voxel composed of mixed tissues, the PV correction problem turns to how to estimate the perfusion contributions of different tissues from multiple noisy measurements. Generally, for magnetic resonance imaging (MRI), the Rician noise model is widely accepted . However, after the label/control difference operation, Gaussian noise is generally considered in the perfusion images of ASL sequence [12, 13]. For the purpose of this study, we first proposed a statistical perfusion model of mixed tissues for the ASL sequence, based on the Gaussian distribution of multiple measurements. With the tissue mixture information obtained from the high-resolution structural image, a structure-based expectation maximization (sEM) scheme was developed to estimate the perfusion contributions of the mixed tissues from multiple measurements.

## Methods

### Perfusion model of a mixed voxel

Considering the low spatial resolution, the GM, WM, and cerebrospinal fluid (CSF) may all contribute to the label/control difference signal, ∆M. No ASL signal typically arises from CSF ; therefore, the perfusion signal ∆M at the spatial position i can be described as

$$\Delta M_{i} = P_{iGM} \Delta M_{iGM} + P_{iWM} \Delta M_{iWM}$$
(1)

where PiGM and PiWM are proportions of GM and WM in the voxel i, respectively. ∆MiGM and ∆MiWM are the difference magnetizations for GM and WM, respectively.

In the current CBF calculation method, the CBF f of a tissue type is obtained by

$$f_{tissue} = \left( {\frac{{\Delta M_{tissue} }}{{M_{0} }}} \right)F_{tissue}$$
(2)

where Ftissue is a tissue-specific parameter, and M0 represents the equilibrium brain tissue magnetization obtained from the M0 image. For a mixed voxel, its CBF comes independently from the GM part ($$f_{GM}^{P}$$) and the WM part ($$f_{WM}^{P}$$) and can be described as

$${\text{CBF}} = f_{GM}^{P} + f_{WM}^{P} = \frac{{P_{iGM} F_{GM} }}{{M_{i0} }}\Delta M_{iGM} + \frac{{P_{iWM} F_{WM} }}{{M_{i0} }}\Delta M_{iWM}$$
(3)

For ASL perfusion studies, PiGM and PiWM can usually be estimated from a high-resolution structural image (e.g., T1 weighted image) of the same subject, and FGM and FWM can be derived from the two-compartment model for the ASL data . Therefore, for a CBF estimation of a mixed voxel, the key problem is to estimate the magnetizations of GM and WM (i.e., ∆MiGM and ∆MiWM) from multiple measurements.

### Statistical perfusion model of mixed tissues

As described in the Introduction section, multiple measurements could be regarded as noisy realizations of the original distribution , and Gaussian noise is generally considered in each measurement of the ASL sequence [12, 13]. Based on the Gaussian distribution of multiple ASL measurements, we first established a statistical perfusion model of mixed tissues.

#### 1. The statistical model of multiple measurements

In the spatial domain, index i is defined to represent the spatial position of a concerned voxel. The intensities of this voxel were acquired by multiple measurements that constitute a column vector Y = {Yit, t = 1, …,T}, where T is the number of multiple measurements collected. Based on the Gaussian assumption, each Yit is a noisy observation of a random variable with a mean of $$\bar{Y}_{i}$$ and a variance of $$\sigma_{i}^{2}$$, i.e.,

$$Y_{it} = \bar{Y}_{i} + n$$
(4)

where n represents statistically independent noise in Yit . Since each measurement of the ASL sequence is independently scanned, the conditional probability of the measurement vector Y can be described as

$$p\left( {{\mathbf{Y}}\left| {\{ \bar{Y}_{i} \} ,\{ \sigma_{i}^{2} \} } \right.} \right) = \prod\limits_{t = 1}^{T} {p\left( {Y_{it} \left| {\bar{Y}} \right.,\sigma_{i}^{2} } \right)}$$
(5)

#### 2. Statistical perfusion model of mixed tissues

The observation Yit contains perfusion contributions from GM and WM. The GM component is denoted by XitGM, with a mean of $$\bar{X}_{iGM}$$ and a variance of $$\sigma_{iGM}^{2}$$. The WM component is denoted by XitWM with a mean of $$\bar{X}_{iWM}$$ and a variance of $$\sigma_{iWM}^{2}$$. Thus, we have

$$p\left( {{\mathbf{X}}\left| {\bar{X}_{iGM} ,\bar{X}_{iWM} ,\sigma_{iGM}^{2} ,\sigma_{iWM}^{2} } \right.} \right) = \prod\limits_{t = 1}^{T} {\left\{ {p\left( {X_{itGM} \left| {\bar{X}_{iGM} ,\sigma_{iGM}^{2} } \right.} \right)p\left( {X_{itWM} \left| {\bar{X}_{iWM} ,\sigma_{iWM}^{2} } \right.} \right)} \right\}}$$
(6)

where X = {XitGM and XitWM, t = 1, …,T} represents a vector of size 2 × T, at position i.

The mean and variance values of each voxel can be calculated by the summation of all contributions at this voxel, i.e.,

$$\bar{Y}_{i} = \bar{X}_{iGM} + \bar{X}_{iWM} \;{\text{and}}\;\sigma_{i}^{2} = \sigma_{iGM}^{2} + \sigma_{iWM}^{2}$$
(7)

By combining the voxel-wise perfusion model in Eq. 3 with the above observation model, we have

$$\bar{X}_{iGM} = P_{iGM} \Delta M_{iGM} \;{\text{and}}\;\bar{X}_{iWM} = P_{iWM} \Delta M_{iWM}$$
(8)
$$\sigma_{iGM}^{2} = P_{iGM} S_{iGM} \;{\text{and}}\;\sigma_{iWM}^{2} = P_{iWM} S_{iWM}$$
(9)

where SiGM and SiWM represent the variance of the GM and the WM signal, respectively. In this study, the PiGM and PiWM, which represent the proportions of GM and WM inside the concerned voxel i, can be estimated from the registered high-resolution structural image, which can be regarded as constants for a concerned voxel.

#### 3. Normal statistical model

For the ASL sequence, the perfusion signal contains GM and WM components. Suppose that each tissue type is independent and follows a Gaussian distribution. Equation 6 becomes

\begin{aligned} & p({\mathbf{X}}\left| {\Delta M_{iGM} ,\Delta M_{iWM} ,S_{iGM} ,S_{iWM} } \right.) \\ & = \prod\limits_{t = 1}^{T} {\left\{ {\left( {\frac{1}{{\sqrt {2\pi P_{iGM} S_{iGM} } }}e^{{ - \frac{{\left( {X_{itGM} - P_{iGM} \Delta M_{iGM} } \right)^{2} }}{{2P_{iGM} S_{iGM} }}}} } \right) \times \left( {\frac{1}{{\sqrt {2\pi P_{iWM} S_{iWM} } }}e^{{ - \frac{{\left( {X_{itWM} - P_{iWM} \Delta M_{iWM} } \right)^{2} }}{{2P_{iWM} S_{iWM} }}}} } \right)} \right\}} \\ \end{aligned}
(10)

The estimation of $$p\left( {{\mathbf{Y}} |\Delta M_{iGM} ,\Delta M_{iWM} ,S_{iGM} ,S_{iWM} } \right)$$ derived from Eq. 5 would generate several nonlinear equations, which are difficult to solve. Given $$\bar{Y}_{i} = \bar{X}_{iGM} + \bar{X}_{iWM}$$ in Eq. 7 and the description in Eq. 10, the EM algorithm may provide an alternative method and effective solution to estimate the model parameters {∆MiGM, ∆MiWM, SiGM, SiWM} based on the structural mixture information derived from a high-resolution image.

### EM algorithm for parameter estimation

In the EM approach [17, 18], the observation Yit is regarded as an incomplete random variable. The XitGM and XitWM are regarded as complete variables, which can reflect the complete perfusion information at each measurement point t for a concerned voxel of position i. The probability distribution of the incomplete data {Yit} can be depicted by the complete data, {XitGM} and {XitWM}, using an integral equation under the condition of {Yit = XitGM + XitWM}:

\begin{aligned} & p\left( {Y_{it} \left| {\Delta M_{iGM} ,\Delta M_{iWM} ,S_{iGM} ,S_{iWM} } \right.} \right) \\ {\kern 1pt} & = \int_{{\left\{ {Y_{it} = X_{itGM} + X_{itWM} } \right\}}} {\left\{ {p\left( {X_{itGM} \left| {\bar{X}_{iGM} ,\sigma_{iGM}^{2} } \right.} \right)p\left( {X_{itWM} \left| {\bar{X}_{iWM} ,\sigma_{iWM}^{2} } \right.} \right)} \right\}dX} \\ \end{aligned}
(11)

In this study, the EM algorithm was used to seek a solution to maximize the conditional expectation of the complete data in Eq. 10. The E-step is to compute the conditional expectation. The M-step subsequently attempts to maximize the expectation of the complete-data log likelihood using the latent variables that were computed in the E-step, given the observations.

E-step This step computes the likelihood p(X|Θ) of the complete data in Eq. 10, given {Yit} and parameter $$\varTheta^{(n)} = \left\{ {\Delta M_{iGM}^{(n)} ,\Delta M_{iWM}^{(n)} ,S_{iGM}^{(n)} ,S_{iWM}^{(n)} } \right\}$$. The conditional expectation is depicted in Eq. 12.

\begin{aligned} Q(\varTheta |\varTheta^{(n)} ) = E_{{Y_{it} = X_{itGM} + X_{itWM} }} [\ln (p({\text{X}}|\varTheta ))|{\text{Y}},\varTheta^{(n)} ] = E_{{_{{Y_{it} = X_{itGM} + X_{itWM} }} }} \left[ { - \frac{1}{2}\sum\limits_{t} {\left\{ {\ln \left( {2\pi P_{iGM} S_{iGM} } \right) + \frac{1}{{P_{iGM} S_{iGM} }}\left[ {X_{itGM}^{2} - 2P_{iGM} \Delta M_{iGM} X_{itGM} + (P_{iGM} \Delta M_{iGM} )^{2} } \right]} \right\}} |Y_{it} ,\varTheta^{(n)} } \right] + E_{{_{{Y_{it} = X_{itGM} + X_{itWM} }} }} \left[ { - \frac{1}{2}\sum\limits_{t} {\left\{ {\ln \left( {2\pi P_{iWM} S_{iWM} } \right) + \frac{1}{{P_{iWM} S_{iWM} }}\left[ {X_{itWM}^{2} - 2P_{iWM} \Delta M_{iWM} X_{itWM} + (P_{iWM} \Delta M_{iWM} )^{2} } \right]} \right\}} |Y_{it} ,\varTheta^{(n)} } \right] = - \frac{1}{2}\sum\limits_{t} {\left\{ \begin{aligned} \ln \left( {2\pi P_{iGM} S_{iGM} } \right) + \frac{1}{{P_{iGM} S_{iGM} }}\left[ {E_{{Y_{it} = X_{itGM} + X_{itWM} }} (X_{itGM}^{2} |Y_{it} ,\varTheta^{(n)} ) - 2P_{iGM} \Delta M_{iGM} E_{{Y_{it} = X_{itGM} + X_{itWM} }} (X_{itGM} |Y_{it} ,\varTheta^{(n)} ) + (P_{iGM} \Delta M_{iGM} )^{2} } \right] + \hfill \\ \ln (2\pi P_{iWM} S_{iWM} ) + \frac{1}{{P_{iWM} S_{iWM} }}\left[ {E_{{_{{Y_{it} = X_{itGM} + X_{itWM} }} }} (X_{itWM}^{2} |Y_{it} ,\varTheta^{(n)} ) - 2P_{iWM} \Delta M_{iWM} E_{{Y_{it} = X_{itGM} + X_{itWM} }} (X_{itWM} |Y_{it} ,\varTheta^{(n)} ) + (P_{iWM} \Delta M_{iWM} )^{2} } \right] \hfill \\ \end{aligned} \right\}} \end{aligned}
(12)

Based on the deduction of the preceding conditional expectation, we have

\begin{aligned} X_{itGM}^{(n)} & = E_{{Y_{it} = X_{itGM} + X_{itWM} }} (X_{itGMt} |Y_{it} ,\varTheta^{(n)} ) \\ & = P_{iGM} \Delta M_{iGM}^{(n)} + \frac{{P_{iGM} S_{iGM}^{(n)} }}{{P_{iGM} S_{iGM}^{(n)} + P_{iWM} S_{iWM}^{(n)} }} \\ & \quad\quad \times \left[ {Y_{it} - (P_{iGM} \Delta M_{iGM}^{(n)} + P_{iWM} \Delta M_{iWM}^{(n)} )} \right] \\ \end{aligned}
(13)
\begin{aligned} X_{itWM}^{(n)} & = E_{{Y_{it} = X_{itGM} + X_{itWM} }} (X_{itWM}^{{}} |Y_{it} ,\varTheta^{(n)} ) \\ & = P_{iWM} \Delta M_{iWM}^{(n)} + \frac{{P_{iWM} S_{iWM}^{(n)} }}{{P_{iGM} S_{iGM}^{(n)} + P_{iWM} S_{iWM}^{(n)} }} \\ & \quad\quad \times {\kern 1pt} \left[ {Y_{it} - (P_{iGM} \Delta M_{iGM}^{(n)} + P_{iWM} \Delta M_{iWM}^{(n)} )} \right] \\ \end{aligned}
(14)
\begin{aligned}(X_{itGM}^{2} )^{(n)} &= E_{{_{{Y_{it} = X_{itGM} + X_{itWM} }} }} [X_{itGMt}^{2} |Y_{it} ,\varTheta^{(n)} ]{\kern 1pt} \\ & = (X_{itGM}^{(n)} )^{2} + \frac{{(P_{iGM} S_{iGM}^{(n)} )\left( {P_{iWM} S_{iWM}^{(n)} } \right)}}{{P_{iGM} S_{iGM}^{(n)} + P_{iWM} S_{iWM}^{(n)} }} \end{aligned}
(15)
\begin{aligned} (X_{itWM}^{2} )^{(n)} &= E_{{Y_{it} = X_{itGM} + X_{itWM} }} [X_{itWM}^{2} |Y_{it} ,\varTheta^{(n)} ]{\kern 1pt} \\ &= (X_{itWM}^{(n)} )^{2} + \frac{{(P_{iGM} S_{iGM}^{(n)} )(P_{iWM} S_{iWM}^{(n)} )}}{{P_{iGM} S_{iGM}^{(n)} + P_{iWM} S_{iWM}^{(n)} }} \end{aligned}
(16)

M-step: This step maximizes the conditional expectation to estimate the next iteration $$\left\{ {\Delta M_{iGM}^{{(n{ + }1)}} ,\Delta M_{iWM}^{{(n{ + }1)}} ,S_{iGM}^{{(n{ + }1)}} ,S_{iWM}^{{(n{ + }1)}} } \right\}$$, which can be described as

$$\frac{\partial Q}{{\partial \Delta M_{iGM} }}|_{{\Delta M_{iGM} = \Delta M_{iGM}^{(n + 1)} }} = 0 \Rightarrow \Delta M_{iGM}^{(n + 1)} = \frac{{\sum\nolimits_{t = 1}^{T} {X_{itGM}^{(n)} } }}{{T \cdot P_{iGM} }}$$
(17)
$$\frac{\partial Q}{{\partial \Delta M_{iWM} }}|_{{\Delta M_{iWM} = \Delta M_{iWM}^{(n + 1)} }} = 0 \Rightarrow \Delta M_{iWM}^{(n + 1)} = \frac{{\sum\nolimits_{t = 1}^{T} {X_{itWM}^{(n)} } }}{{T \cdot P_{iWM} }}$$
(18)
$$S_{iGM}^{(n + 1)} = \frac{{\sum\nolimits_{t = 1}^{T} {\left[ {(X_{itGM}^{2} )^{(n)} - 2X_{itGM}^{(n)} P_{iGM} \Delta M_{iGM}^{(n)} + (P_{iGM} \Delta M_{iGM}^{(n)} )^{2} } \right]} }}{{T \cdot P_{iGM} }}$$
(19)
$$S_{iWM}^{(n + 1)} = \frac{{\sum\nolimits_{t = 1}^{T} {\left[ {(X_{itWM}^{2} )^{(n)} - 2X_{itWM}^{(n)} P_{iWM} \Delta M_{iWM}^{(n)} + (P_{iWM} \Delta M_{iWM}^{(n)} )^{2} } \right]} }}{{T \cdot P_{iWM} }}$$
(20)

Based on the proposed sEM algorithm, we can estimate ∆MiGM and ∆MiWM using the multiple measurements of the ASL sequence.

### Implementation of the sEM scheme for PV correction

The implementation of the proposed sEM scheme for PV correction can be summarized as follows:

1. 1.

Segmentation of high-resolution structural image. The segmented results and ASL data are co-registered. For each mixed voxel at position i, the percentages of GM and WM, PiGM and PiWM, were obtained.

2. 2.

Initialization of the model parameters $$\left\{ {\Delta M_{iGM}^{(0)} ,S_{iGM}^{(0)} ,\Delta M_{iWM}^{(0)} ,S_{iWM}^{(0)} } \right\}$$.

3. 3.

Constitute a column vector with all measurements of the mixed voxel at position i.

4. 4.

Iterative estimation of GM and WM components for the mixed voxel at position i using the column vector in step (3), following Eqs. 1720.

5. 5.

Repeat steps (3) and (4) for the next voxel until all the voxels are corrected.

## Material and evaluation

In this study, the performance of the proposed sEM scheme was evaluated by both digital simulations and clinical data. The two simulations listed below were designed to evaluate its performance quantitatively, especially with regards to noise reduction, lesion detection, and its potential to estimate CBF from fewer measurements. After the simulation studies, the in vivo ASL data were used to evaluate the clinical feasibility.

### Simulation 1

In this simulation, a digital head phantom was generated from a structural MRI brain dataset with a voxel size of 1 × 1 × 1 mm3. After the normalization and segmentation of the MRI data using SPM8 software, the posterior probability images of GM and WM were generated. Next, the images were masked to remove the voxels with probabilities lower than 0.1 [7, 9]. The head phantom was simulated as follows:

1. 1.

The probability images were resampled to a size of 60 × 72 × 60, with a spatial resolution of 3 × 3 × 3 mm3 using SPM8.

2. 2.

Across the whole brain, the WM region was simulated as 20 mL/100 g/min.

3. 3.

The GM was simulated as 60 mL/100 g/min, with a hypo-perfused region (30 mL/100 g/min) and a hyper-perfused region (90 mL/100 g/min). Both of the regions were spherical regions with a radius of 5.

4. 4.

Based on the probability images and the signals of GM and WM, the perfusion signal of each voxel in the 3D perfusion image was generated according to Eq. 1.

5. 5.

It was reported that the noise level of the ASL data ranges from 6.7 to 13.2 according to different labeling schemes and readout sequences . To evaluate the noise impact on PV correction, three different levels of Gaussian noise, with a standard deviation (std) of 5, 10 and 15, respectively, were added into the 3D perfusion image to generate low-, middle-, and high-noise realizations. The highest noise was approximately 25% (15/60) of the GM signal.

6. 6.

Generally, the number of label/control pairs is set as 40–60. To evaluate the proposed method, 40 noisy realizations were generated for each ASL sequence.

### Simulation 2

To evaluate the benefit of PV correction on the lesion detection of small CBF alterations, in this simulation, three regions with different sizes and simulated values, instead of the two regions used in step (3) of Simulation 1, were simulated inside the homogeneous GM tissues: (1) a spherical region of radius 5 with CBF of 75 mL/100 g/min, (2) a 3 × 3 × 3 cubic region with CBF of 45 mL/100 g/min, and (3) a 2 × 2 × 2 cubic region with CBF of 75 mL/100 g/min. The difference between the three regions and the homogeneous GM region were selected from the high std of noise, i.e., 15.

### In vivo data

To test the feasibility of PV correction on in vivo ASL data, the ASL scans were collected from three healthy subjects, which were acquired by a Siemens 3T scanner using the pseudo-continuous ASL perfusion imaging sequence with gradient-echo echoplanar imaging (EPI). The acquisition parameters were TR = 4 s, TE = 11 ms, FOV = 220 × 220 mm2, voxel size = 3.4 × 3.4 × 5 mm3, matrix = 64 × 64 × 20, flip angle = 90°, and post labeling delay = 1.5 s. Forty label/control pairs were acquired. A high-resolution structural image was also acquired with the following parameters: TR = 1900 ms, TE = 2.9 ms, FOV = 250 × 250 mm2, matrix = 256 × 256 × 176, and flip angle = 90°.

The ASL and structural images were preprocessed using SPM8. For each subject, the ASL images was realigned separately for the label and control image series. After realignment, the images were normalized, followed by pair-wise subtraction. The corresponding structural image was normalized and segmented to generate probability images of GM and WM, which were later masked with probabilities lower than 0.1. Finally, the probability images were co-registered with ASL data to obtain PiGM and PiWM at each position i, using a transformation of the structural and ASL coordinates with an MNI coordinate.

### Comparison of PV correction

As is well-known, the EM algorithm is quite sensitive to the initialization. Considering the limited number of measurements and the intensive computation load of the EM algorithm, a relatively accurate initialization from an estimation that uses an uncorrected image or other spatial PV correction method (e.g., the LR method) would lead to accurate estimations and fast convergence. To compare the effect of the PV correction using different methods, the simulated data and the in vivo data were all analyzed using:

1. 1.

No correction. The averaged image was used as the result.

2. 2.

The LR method. The averaged image was used to separately estimate the GM and the WM CBF maps using the LR method with a 5 × 5 × 1 regression kernel, which was suggested to provide the best compromise between smoothing and PV correction [5, 7].

3. 3.

The sEM method, which is the EM algorithm initialized with an estimation from no correction. In this method, $$\Delta M_{iGM}^{(0)}$$ and $$\Delta M_{iWM}^{(0)}$$ were set as the mean value of GM and WM regions from no correction, and $$S_{iGM}^{(0)}$$ and $$S_{iWM}^{(0)}$$ were set as the std of GM and WM. The iteration number was set as 100 to ensure the convergence.

4. 4.

The sEM-LR method, which is the EM algorithm initialized with the LR method. In this method, each 3D difference image was first corrected with the LR method to obtain the initialization of $$\left\{ {\Delta M_{iGM}^{(0)} ,\Delta M_{iWM}^{(0)} ,S_{iGM}^{(0)} ,S_{iWM}^{(0)} } \right\}$$. With this initialization, the GM and WM maps were estimated using the sEM method. The iteration number was also set as 100 to ensure the convergence.

For the simulation data, the root mean square error (RMSE) analysis was performed for a quantitative evaluation of these correction methods.

For the in vivo data, the GM CBF ratio, which is the ratio between the estimated GM CBF and the mean GM CBF of the uncorrected maps, was calculated for each voxel. This index can avoid the bias introduced from a different calibration method in which the CBF value is calculated and permits the assessment of the relative CBF changes after correction .

### The region of interest (ROI) analysis

In this study, the consistency of the mean GM CBF across the whole range of GM probabilities was used to quantitatively evaluate the estimated results from different PV correction methods. To this aim, nine ROIs were automatically defined based on the GM probability images, with the probability range between [10–20%], [20–30%],…, [90–100%], respectively. Next, the mean value of GM CBF in each ROI was calculated. It should be noted that the less independent are the GM CBF values from the GM probability, the better the performance of PV correction is.

## Results

### Simulation results

Figure 1 shows the middle slice of the GM CBF estimation for Simulation 1 using no correction, LR, sEM, and sEM-LR methods. Clearly, the CBF maps derived from the LR, sEM, and sEM-LR methods outperformed those of no correction, with less noise and better restoration. At the edges of hypo- and hyper-CBF regions, the GM map that was estimated by the LR method exhibited a visible smoothing effect.

Figure 2 shows the results of ROI analysis using Simulation 1 when the different PV-corrected methods were performed. It demonstrated that the GM CBF estimation using no correction was underestimated, compared with the ground truth. Corrected by the LR and sEM-LR methods, the GM CBF curves of different GM probabilities were almost consistent with the true line, while that of sEM method was a little underestimated at the relative low GM probability. The performance of the LR and the two sEM-based methods seems to be less affected by the noise level.

To illustrate the effect of different correction methods on the CBF accuracy under different noise levels, the profiles of the lines passing the centers of the hypo- and hyper-CBF regions of the GM CBF maps are shown in Fig. 3, which demonstrates that the sEM and sEM-LR methods provided accurate GM CBF estimations with preserved details and tissue interfaces but are affected by the noise level. Table 1 gives the RMSE values of the estimated CBF maps and the true map, and the differences between them indicated that the sEM-LR method outperformed the LR method at different noise levels.

The effect of the PV correction on lesion detection is shown in Fig. 4. It is obvious that, although the alterations were small, all of the regions with CBF alterations can be detected by using two sEM-based methods, even if the std of the noise was same as the CBF alteration. However, the two small regions (region 2 and region 3 in Fig. 4) were difficult to detect when corrected by the LR method.

Figure 5 demonstrates the GM CBF maps (middle slice) that were estimated from fewer measurements, which indicate that with the increase of measurement numbers, the CBF estimation was more accurate and was less affected by noise. The RMSE values of the CBF maps that were estimated from different numbers of measurements are listed in Table 2, which also illustrate that the restoration was better with the increased number of multiple measurements. In most cases, the RMSEs using the sEM-LR method with fewer measurements (Table 2) were lower than those of the LR method with normal measurements (the corresponding RMSE shown in Table 1).

The computation times of each correction method to correct Simulation 1 were compared using the same computer (Intel CPU E3-1240, RAM of 16G). The computation time of the LR method for the 60 × 72 × 60 averaged image was 19.2 s. With the stopping criterion of 100 iterations, the computation costs for the sEM and sEM-LR methods were 177 s and 982 s, respectively. With the stopping rule of the difference between two adjacent iterations less than 0.001, the time costs of them were 4 s and 792 s, respectively. It should be noticed that the majority time of the sEM-LR was used for the initialization of all spatial label/control difference images using the LR method, which was about 790 s.

### In vivo data

Figure 6 gives the GM CBF ratio of three subjects by using different correction methods. For a better demonstration of the results, the regions enclosed within dotted boxes were zoomed. Compared with the results without correction and estimated from LR method, the proposed sEM and sEM-LR methods reserved more details, especially at the tissue interface.

Figure 7 shows the ROI analysis of the ASL data using different methods. For each subject, the results of the LR and the two sEM-based methods demonstrate less variation (lower standard deviation) than that of the uncorrected data (Table 3), which indicate less independence of the GM CBF values from the GM probability.

## Discussion

The present study proposed a sEM scheme for the PV correction of the ASL sequence. For an accurate estimation of CBF, a statistical perfusion model of mixed tissues was first established. Then, based on the prior tissue mixture obtained from a high-resolution structural image, a structure-based EM algorithm (sEM scheme) was proposed to estimate the perfusion contributions of GM and WM tissues of the mixed voxels from multiple measurements of the ASL sequence. When the contributions of different tissues were estimated, the PV effect embedded in the multiple measurements was naturally resolved.

Different from the previous PV correction studies, the proposed method innovatively utilizes multiple measurements of label/control differences (perfusion images), instead of using the simple averaged image, to estimate the CBF contribution of the GM and WM components in each mixed voxel. The evaluation using computer simulations and the in vivo data demonstrated its superiority in PV correction, especially in the following aspects: (1) Edge preservation. Since the CBF contributions were estimated iteratively from the multiple measurements of a mixed pixel, with less influence from neighboring voxels, the EM estimation was superior in edge preservation and could detect small lesions with a radius of approximately 3.4 mm (calculated from a spherical volume of 2 × 2 × 2 m3 cube). (2) Noise suppression. Unlike the simple averaging of multiple noisy measurements, the sEM scheme restored the GM and WM components from a series of noisy realizations with Gaussian distribution. Thus, the scheme could not only suppress noise, but also could detect small CBF signals effectively, even if strong noise was applied. (3) Fast scan. The CBF estimation using fewer measurements indicated that the proposed method could achieve reasonable imaging quality with fewer label/control pairs and have the potential to shorten the scan time.

Unlike our previous work in which the EM algorithm was used to estimate the tissue mixture inside a mixed voxel [18, 20], in this study, we attempted to integrate the 3D structural image with perfusion series and develop a new sEM scheme for the perfusion estimation of different tissues in a mixed voxel from the multiple measurements of the ASL sequence. Since the contributions of GM and WM to the perfusion signal are independent and different, the proposed sEM scheme could estimate their different contributions effectively. However, if they are correlated or contribute same to the perfusion signal, the sEM method would not help, in which the simple averaging should be good enough.

It is known that the EM algorithm is quite sensitive to the initialization. If the initial values of the model parameters, such as ∆MiGM and ∆MiWM, can be set as close as possible to the true values, better estimations could be obtained with fast convergence. To evaluate the effect of parameter initialization on the CBF estimation, the EM algorithm initialized with parameters estimated without correction and those estimated using the LR method were performed on both simulated and in vivo data. The results indicated that both sEM-based methods (sEM and sEM-LR) outperformed the LR method, while the sEM-LR method performed better than the sEM method only at relatively low GM probabilities (Fig. 2). Following the Markov random field model, the perfusion of a voxel is generally affected by neighboring voxels . Since the proposed sEM method only considers perfusion correction from multiple measurements of the same voxel, a more accurate CBF estimation could be expected if spatial correction is considered further. Therefore, the combination of the proposed sEM with spatial prior obtained from the LR method, i.e., the sEM-LR method, could achieve better performance with the consideration of a spatial neighborhood.

Considering the iterative nature of the EM algorithm, the computation load of different methods was compared. The results indicated that the time cost of the sEM correction was comparable with other methods if a reasonable stopping criterion was used. The major cost of the sEM-LR method came from the initialization of all spatially different images by using the LR method, and not from EM optimization itself. The results also suggest that the use of the difference between two adjacent iterations that were less than 0.001 as the stopping criterion could reduce the computation time remarkably, because most voxels without the tissue mixture could reach the criterion very quickly. If parallel computation was performed, the computation time will be further greatly reduced.

Several limitations of this study should be addressed. Firstly, the proposed method needs multiple measurement information to correct PV effect, thus, this method is more suitable for the ASL sequence with time series, not for 3D ASL sequence. Secondly, the present study assumed that the voxels located at the same 3D spatial position differed only in noise. In practice, the distribution may be affected by temporal CBF variation, which may induce a bias of the CBF estimation for the in vivo data. In this study, we focus on the feasibility to use multiple measurements for an accurate CBF estimation under this assumption, and further studies will be performed to investigate the PV correction by using multiple measurements with consideration of temporal CBF variation. Although further improvement is required, this study validates the proposed statistical perfusion model and demonstrates the effectiveness and necessity of using inherent perfusion information in multiple measurements for PV correction of the ASL sequence.

## Conclusions

In this study, we proposed a statistical perfusion model of mixed tissues for each voxel of the ASL data. Based on this model, the sEM scheme was developed to estimate the contributions of different tissues to the perfusion signal of the mixed voxel with its multiple measurements. Compared to the traditional PV-corrected method, the proposed sEM-based method performs better in edge preservation, noise suppression, and lesion detection while demonstrating the potential to estimate CBF within a shorter scanning time. The results also indicated the effectiveness of using inherent perfusion information in multiple measurements for PV correction of the ASL sequence.

## Abbreviations

ASL:

arterial spin labeling

CBF:

cerebral blood flow

PV:

partial volume

EM:

expectation maximization

MRI:

magnetic resonance imaging

LR:

linear regression

## References

1. 1.

Detre JA, Leigh JS, Williams DS, Koretsky AP. Perfusion imaging. Magn Reson Med. 1992;23(1):37–45.

2. 2.

Alsop DC, Detre JA. Reduced transit-time sensitivity in noninvasive magnetic resonance imaging of human cerebral blood flow. J Cereb Blood Flow Metab. 1996;16(6):1236–49.

3. 3.

Zhao MY, Mezue M, Segerdahl AR, Okell TW, Tracey I, Xiao Y, Chappell MA. A systematic study of the sensitivity of partial volume correction methods for the quantification of perfusion from pseudo-continuous arterial spin labeling MRI. Neuroimage. 2017;162:384–97.

4. 4.

Jezzard P, Chappell MA, Okell TW. Arterial spin labeling for the measurement of cerebral perfusion and angiography. J Cereb Blood Flow Metab. 2018;38(4):603–26.

5. 5.

Asllani I, Borogovac A, Brown TR. Regression algorithm correcting for partial volume effects in arterial spin labeling MRI. Magn Reson Med. 2008;60(6):1362–71.

6. 6.

Ahlgren A, Wirestam R, Petersen ET, Stahlberg F, Knutsson L. Partial volume correction of brain perfusion estimates using the inherent signal data of time-resolved arterial spin labeling. NMR Biomed. 2014;27(9):1112–22.

7. 7.

Liang X, Connelly A, Calamante F. Improved partial volume correction for single inversion time arterial spin labeling data. Magn Reson Med. 2013;69(2):531–7.

8. 8.

Petr J, Schramm G, Hofheinz F, Langner J, van den Hoff J. Partial volume correction in arterial spin labeling using a look-locker sequence. Magn Reson Med. 2013;70(6):1535–43.

9. 9.

Chappell MA, Groves AR, MacIntosh BJ, Donahue MJ, Jezzard P, Woolrich MW. Partial volume correction of multiple inversion time arterial spin labeling MRI data. Magn Reson Med. 2011;65(4):1173–83.

10. 10.

Petr J, Ferré JC, Gauvrit JY, Barillot C. Improving arterial spin labeling data by temporal filtering. In: Proceedings of SPIE. 2010; 7623:76233B.

11. 11.

Gudbjartsson H, Patz S. The Rician distribution of noisy MRI data. Magn Reson Med. 1995;34(6):910–4.

12. 12.

Bibic A, Knutsson L, Stahlberg F, Wirestam R. Denoising of arterial spin labeling data: wavelet-domain filtering compared with Gaussian smoothing. MAGMA. 2010;23(3):125–37.

13. 13.

Wink AM, Roerdink JB. Denoising functional MR images: a comparison of wavelet denoising and Gaussian smoothing. IEEE Trans Med Imaging. 2004;23(3):374–87.

14. 14.

Golay X, Hendrikse J, Lim TC. Perfusion imaging using arterial spin labeling. Top Magn Reson Imaging. 2004;15(1):10–27.

15. 15.

Wang J, Alsop DC, Li L, Listerud J, Gonzalez-At JB, Schnall MD, Detre JA. Comparison of quantitative perfusion imaging using arterial spin labeling at 1.5 and 4.0 Tesla. Magn Reson Med. 2002;48(2):242–54.

16. 16.

Liu Y, Li B, Zhang X, Zhang L, Liang Z, Lu H: Partial volume correction for arterial spin labeling data using spatial-temporal information. In: SPIE Medical Imaging. 2015; 941325–941328.

17. 17.

Dempster AP, Laird NM, Rubin DB. Maximum likelihood from incomplete data via the EM algorithm. J R Stat Soc Ser B (Methodol). 1977;39:1–38.

18. 18.

Liang Z, Wang S. An EM approach to MAP solution of segmenting tissue mixtures: a numerical analysis. IEEE Trans Med Imaging. 2009;28(2):297–310.

19. 19.

Vidorreta M, Wang Z, Rodriguez I, Pastor MA, Detre JA, Fernandez-Seara MA. Comparison of 2D and 3D single-shot ASL perfusion fMRI sequences. Neuroimage. 2013;66:662–71.

20. 20.

Wang S, Li L, Cohen H, Mankes S, Chen JJ, Liang Z. An EM approach to MAP solution of segmenting tissue mixture percentages with application to CT-based virtual colonoscopy. Med Phys. 2008;35(12):5787–98.

21. 21.

Geman S, Geman D. Stochastic relaxation, gibbs distributions, and the bayesian restoration of images. IEEE Trans Pattern Anal Mach Intell. 1984;6(6):721–41.

## Authors’ contributions

Conceived and designed the experiments: YL RL ZL HL. Performed the experiments: YL BL. Analyzed the data: YL ZW. Contributed reagents/materials/analysis tools: YL ZW HL. Wrote the paper: YL. All authors read and approved the final manuscript.

Not applicable.

### Competing interests

The authors declare that they have no competing interests.

### Availability of data and materials

The datasets of the present study are available from the corresponding author upon reasonable request.

Not applicable.

### Ethics approval and consent to participate

The study was conducted according to the principles in the Declaration of Helsinki and approved by the Institutional Board of the Fourth Military Medical University. All subjects received a comprehensive description of the MRI study and gave voluntarily written informed consent before entering the study.

### Funding

This research was supported by grants from National Natural Science Foundation of China (Grant 81871424), the Military Science Foundation (Grant BWS14J038), National Key Research and Development Program of China (Grant 2017YFC0107403).

### Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Author information

Authors

### Corresponding author

Correspondence to Hongbing Lu.

## Rights and permissions
{"ft_lang_label":"__label__en","ft_lang_prob":0.90851814,"math_prob":0.99142057,"size":34889,"snap":"2020-45-2020-50","text_gpt3_token_len":7975,"char_repetition_ratio":0.16812384,"word_repetition_ratio":0.06729893,"special_character_ratio":0.22634642,"punctuation_ratio":0.12292052,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.9972069,"pos_list":[0,1,2],"im_url_duplicate_count":[null,4,null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2020-10-22T01:24:41Z\",\"WARC-Record-ID\":\"<urn:uuid:c8e457ed-56fb-4d07-91a0-5640bd9838d0>\",\"Content-Length\":\"237560\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:7681d3b2-57f9-431a-8346-09fa2bd6b1d2>\",\"WARC-Concurrent-To\":\"<urn:uuid:6cd50421-96e1-4c93-99b1-3a4a56849a75>\",\"WARC-IP-Address\":\"151.101.200.95\",\"WARC-Target-URI\":\"https://biomedical-engineering-online.biomedcentral.com/articles/10.1186/s12938-019-0631-8\",\"WARC-Payload-Digest\":\"sha1:7CH6ZTWZBUPH56IDR6LTPK72KX6NBJOS\",\"WARC-Block-Digest\":\"sha1:5ALOOZV7SZ6TJ4X2XZZ65FSN6VUQ65JI\",\"WARC-Identified-Payload-Type\":\"text/html\",\"warc_filename\":\"/cc_download/warc_2020/CC-MAIN-2020-45/CC-MAIN-2020-45_segments_1603107878662.15_warc_CC-MAIN-20201021235030-20201022025030-00667.warc.gz\"}"}