Inferring CT perfusion parameters and uncertainties using a Bayesian approach

Introduction

Stroke is the third leading cause of death in the US, killing about 140,000 Americans each year (1). There are two types of strokes, ischemic and hemorrhagic, for which the medical management can be quite different. About 87% of all cases are ischemic strokes (2), meaning there is a complete or partial blockage of the blood supply from a cerebral artery. Standard treatment options for ischemic strokes are to either dissolve or remove the blood clot (i.e., thrombolysis or thrombectomy). Early treatment decisions are crucial for the survival and recovery of stroke patients. To aid the diagnosis and selection of the appropriate treatment, the location and volume of the infarct core and penumbra are often calculated—the infarct core indicates the brain region that is not salvageable, and the penumbra indicates the brain region that is at risk of progression to infarction. The penumbra usually surrounds the ischemic core and may still be salvageable. Classification of the infarct core and penumbra is often derived by thresholding certain brain hemodynamic parameters, e.g., cerebral blood volume (CBV, mL/100 g), cerebral blood flow (CBF, mL/100 g/min), mean transit time (MTT, s), and time to peak (TTP, s). For example, in (3), the ischemic core was defined as regions with relative MTT ≥145% of the normal tissue and absolute CBV <2.0 mL/100 g.

Computed tomography perfusion and dynamic contrast-enhanced magnetic resonance imaging are the primary imaging techniques used to obtain the hemodynamic parameters of the brain. In this study we are interested in CT perfusion, a procedure performed after the administration of intravenous iodine contrast to monitor the first pass of a contrast bolus through the cerebral circulation. Based on the well-known indicator dilution theory (4), from a CT perfusion scan one can obtain the hemodynamic functional parameters related to the blood passage in the tissue, including CBF, CBV, MTT, and TTP. Various approaches, including singular value decomposition (SVD)-based methods, have been applied to solve this inverse problem. In the clinical setting, perfusion analysis is usually done with SVD-based methods, the implementation of which varies across vendors.

Unfortunately, SVD-based methods have some known problems, which make their derived hemodynamic parameters unreliable. First, they make assumptions about the ideal underlying physiological model. For example, the residual function is assumed to be an impulse function, which is not valid for a real scan. Because of this, the estimated parameters can be biased. As shown in (5,6), the SVD-based method can produce biased CBF and MTT estimates. Several authors have proposed advanced deconvolution methods, such as Tikhonov regularization, to address this problem (7-9). Second, these methods often cannot well tolerate imperfect measurements, including relatively long scan intervals (10,11) and truncated measurements in time (12,13). Third, these methods are deterministic, meaning they provide an estimate of the mean of the parameter of interest but not its probability distribution. Factors that affect the reliability of the measurement in a CT perfusion scan, such as the protocol design and perturbance of data, therefore cannot be accounted for in the estimation. For example, a too-long scan interval reduces the available samples for estimation and hence the accuracy. In other cases, patient movement could produce inconsistent voxel attenuation measurements that would inevitably increase the uncertainty of the estimation. Post-processing (e.g., spatial and temporal filtering) may introduce further uncertainty of estimation. For quantification, it is desirable to account for such randomness with proper probability theory (14). Therefore, a probabilistic inference method seems warranted to infer both the hemodynamic parameters and their associated uncertainties.

In this paper, we propose such an accurate and probabilistic method and show its potential value in clinical application. In Section “Methods”, we formulate the estimation problem as a Bayesian inference problem. A variational technique, mean-field approximation, is applied, which makes the problem more tractable as a classic expectation-maximization problem. The posterior distributions (with mean estimate and variance) of CBF and MTT can be obtained. Their coefficients of variation (CoV) are used to indicate the uncertainty in inference, which reflects the reliability of a given measurement. In Section “Experiments and Results”, we first demonstrate the superiority of the proposed method over the conventional method in simulation studies where ground truth is available. We also demonstrate the potential usage of the proposed method in evaluating different suboptimal scan conditions. Finally, patient scans were processed with the proposed method. The results are discussed in Section “Discussion”.

Methods

We first present the implementation of the SVD-based methods. Then the forward maximum log-likelihood model of CT perfusion is presented. We then detail the proposed Bayesian inference algorithm and the quantitative analysis design.

SVD-based methods

In a CT perfusion scan, a mask image is first acquired. The enhanced image at a single time point can be obtained by subtracting the mask image from an image acquired later in the scan. The enhancement of tissue depends on the iodine concentration, which at one voxel can be represented as:

g(t)=h(t)∗Ca(t)[1]

where h(t) is the flow-scaled residual function, indicating the fraction of contrast that remains in the voxel at time t after its arrival; Ca(t) is the arterial input function (AIF), indicating the contrast change in the arteries; and * represents the convolution operation. Eq. [1] describes the well-known indicator dilution theory (4) (Figure 1). Based on this theory, one can obtain h(t) from a CT perfusion scan, which can be used to derive the hemodynamic functional parameters related to the blood passage in the tissue, including CBF, CBV, MTT, and TTP which are defined as:

CBF=1ρ⋅max(h(t))[2]

CBV=1ρ∫0∞h(τ)dτ[3]

MTT=CBV/CBF[4]

TTP=argmaxtg(t)[5]

Figure 1 The indicator dilution theory model (A) describes the blood supply to the tissue at one voxel in the brain. After contrast injection, the bolus enters the voxel via an arterial inlet Ca(t), is dilated in the capillary bed, and leaves the voxel via the venous outlet. Due to the limited spatial resolution of the CT image, Y(t) contain numerous capillary beds as well as arterioles and venules. (B) The residual function h(t) that describes the hemodynamic functions in the physiological model. CBF, CBV, and MTT can be derived from h(t). CBF, cerebral blood flow; MTT, mean transit time.

where ρ is the mean density in the capillary bed, which is set to 1.04 g/mL in this study (15). A detailed explanation of how Eqs. [2-5] is derived can be found in (16).

A straightforward way to obtain the residual function h(t) is to perform deconvolution on Eq. [1]. In practice, g(t) is sampled at discrete time points. Thus, the deconvolution problem is an inverse algebraic problem, which can be solved by SVD. The deconvolution is performed as illustrated in (17). To suppress noise in the resulting residual function, the least significant eigenvectors are consecutively removed until h(t) has an oscillation index below a certain threshold. In this study, a threshold of 0.095 is used as suggested in (17). To further suppress the spatial noise, an alternate deconvolution method with Tikhonov regularization is also implemented using the Regularization Tools package (18). The parameter selection and regularization setup are as suggested in (7).

Forward maximum log-likelihood model

Starting from Eq. [1], the measurement process in a CT perfusion scan can be represented as

Y(t)=g(t)+ε(t)=h(μ,t)∗Ca(t)+ε(t)[6]

where Y(t) is the enhanced tissue time curve at the voxel of interest, µ represents the hemodynamic parameters to be estimated, h(µ,t) is the flow-scaled residual function, and ε(t) is the noise. The flow-scaled residual function h(µ,t) can be further written as a function of the CBF, tissue density ρ, and residual function r(t):

h(μ,t)=CBF⋅ρ⋅r(t)[7]

According to (19–21), r(t) can be modelled as a box shaped function with exponential decay. Hence, h(µ,t) becomes:

{h(μ,t)=CBF⋅ρ⋅exp(−(t−T0)MTT−T0),t≥T0h(μ,t)=CBF⋅ρ,t<T0[8]

where T₀=0.632×MTT. In this study, we consider µ to contain two parameters, CBF and MTT. Other parameters can be calculated; for example, CBV can be calculated as CBV = CBF × MTT, and TTP can be derived directly from the time-attenuation curve (TAC) as the time when the peak enhancement is reached, as shown in Eqs. [2-5]. There is no bolus arrival delay modelled here. The convolution in Eq. [6] can be expanded by sampling the integral of the arterial enhancement curve:

Y(t)=∫0tCa(τ)⋅h(μ,t−τ)dτ+ε(t)[9]

As explained in the introduction, Y(t) in practice is sampled at N discrete time points. Without knowledge of the measurement errors, ε is then modelled as a zero-mean Gaussian distributed in each frame. Since the noise in each frame is temporally independent, we can express the discrete version of Eq. [9] as:

yi=gi(μ)+εiεi∼Ni(0,ωis2)[10]

where the index of the time frame is i, and ω_i is proportional to the number of the detected photon counts (duration of the ith frame as a surrogate). s² is the variance at a given frame. From Eq. [10] we can build the data mismatch term, and the inference problem becomes one of finding the estimate that maximizes the log-likelihood

P(Y|μ)∝∏i=1Nexp(−[yi−gi(μ)]22ωis2)[11]

logP(Y|μ)=∑iN−12logωis2−12(yi−gi(μ))T1ωis2(yi−gi(μ))+C0[12]

μ^=argmaxμlogP(Y|μ)[13]

where C₀ is a constant term. In the case where the length of each frame is the same, one can set ω_i=1. The above optimization problem is essentially equivalent to maximum likelihood estimation, which can be solved by nonlinear least squares fitting, and the way to perform the inference is similar to that of previous perfusion estimation methods (22).

Bayesian inference algorithm

Now we extend the above maximum likelihood estimation to a full Bayesian inference, which provides the probability distribution rather than just the maximum likelihood estimate of a given parameter in Eq. [13]. The derived distribution will inherently contain the measurement uncertainty of each parameter. According to Bayes’ theorem, the posterior of the estimate µ over measurement Y is:

P(μ|Y)=P(Y|μ)⋅P(μ)P(Y)[14]

where P(Y|µ) is the likelihood in Eq. [12], and P(µ) is prior. The fundamental difference here from Eq. [12] is that we are no longer interested in the maximum likelihood estimate but the full distribution of µ. Solving an inference problem within the Bayesian framework involves maximizing the posterior probability P(µ|Y) of µ given the observed data Y. In practice, computing the posterior P(µ|Y) in a closed form is often intractable. An approximation is therefore often used, for which the Markov chain Monte Carlo (MCMC) and variational Bayesian are two popular choices. To achieve a reasonable computational load, the approximate but efficient variational Bayesian method is used here. The basic idea is to find a simple analytical distribution q(µ|Y) to approximate the intractable posterior P(µ|Y) such that the Kullback-Leibler (KL) divergence of these two is minimized. Given the probabilistic model Eq. [14], the log evidence can be written as:

logp(Y)=∫q(μ|Y)logp(Y)dμ=∫q(μ|Y)logP(Y|μ)P(μ)P(μ|Y)dμ=∫q(μ|Y)logP(Y,μ)q(μ|Y)P(μ|Y)q(μ|Y)dμ=∫q(μ|Y)logq(μ|Y)P(μ|Y)dμ+∫q(μ|Y)logP(Y,μ)q(μ|Y)dμ=E∗[logq(μ|Y)P(μ|Y)]+E∗[logP(Y,μ)q(μ|Y)]=KL+ELBO[15]

where P(Y,µ) is the joint probability, E* denotes the expectation with respect to q(µ|Y), ELBO is the lower bound of the evidence, and KL is the divergence between approximate q(µ|Y) and P(µ|Y). Because the KL divergence is always positive, ELBO provides a lower bound on the log-likelihood evidence. Therefore, one can maximize ELBO by finding the distribution q(µ) that approximates the posterior close to the true posterior (or minimizes KL):

q(^μ)=argmaxq(μ)ELBO=argmaxq(μ)E∗[logq(μ|Y)logP(Y,μ)][16]

now we have a second KL divergence in the square brackets, which is always positive, and the optimal solution can be found by equating the numerator and denominator.

Considering the numerator, we first need to ensure it is tractable. Variational Bayesian factorizes q(µ|Y) (in physics this is known as mean-field approximation) to sort the parameters into separate groups, each with their own approximate posterior (23). Note that this approximation does not imply that the parameters are uncorrelated for a given measurement. A simple choice of the distribution of the approximated true posterior is multivariate Gaussian. The mean-field approximation theory assumes that the posterior distribution of each parameter is Gaussian. Supposing that we have a number, k, of parameters to be estimated, the approximate posterior has the following analytical expression:

q(μ|Y)=∏kN(μk,mk,σk)[17]

logq(μ|Y)=∑k−12μk1σkμk+μk1σkmk+C1=−12μTσμ+μTσμ+C1[18]

where σ is a matrix with 1/σ_k as the diagonal elements, m_k is the mean estimate of each parameter, and C₁ is a constant term. Then, for the denominator, we can insert the likelihood from Eq. [12]. As for the prior, we choose the conjugate prior—the prior is said to be conjugate to the likelihood if and only if the factorized posterior has the same parametric form (24). Therefore, the conjugate prior has a multivariable Gaussian distribution with mean m₀ and variance σ02. The numerator thus becomes:

logP(Y,μ)=logP(Y|μ)P(μ)=∑i[−12(yi−gi(μ))T1ωis2(yi−gi(μ))]−12(μ−m0)T1σ0(μ−m0)+C2[19]

where C₂ is a constant term. To ensure tractability and to allow the method to generalize to any nonlinear model, g(µ) is approximated by a first-order Taylor expansion about the mode (mean) of the posterior:

gi(μ)≈gi(m)+Ji(μ−m)(J)j,k=d(g(μ)j)dμk|μ=m[20]

where J is the Jacobian. With such linearization applied to g(µ), Eq. [19] becomes:

logP(Y,μ)≈∑i[−12(yi−gi(m)+Ji(μ−m))T1ωiS2(yi−gi(m)+Ji(μ−m))]−12(μ−m0)T1σ0(μ−m0)[21]

By equalizing Eq. [21] and Eq. [17], we have updated equations for the mean m and σ:

σnew=σ0+∑i1ωisJiTJimnew=σ0m0+∑i1ωisJi(yi−gi(m)+Jim)σ0+∑i1ωisJiTJi[22]

The above derivation matches the description in (24,25): under certain conditions, a Bayesian inference problem reduces to an optimization problem, and the optimal solution can be found by an algorithm similar to expectation maximization. Specifically, the expectation term in Eq. [16] is an expectation step, and the updates in Eq. [22] are a maximization step. The values of the parameters are calculated based on the current values, and these values then used for the next iteration and so on until convergence. Because the independence of the estimated parameters is assumed, each parameter in µ (CBF and MTT) can be updated simultaneously to iteratively maximize the KL. Convergence is determined by observing the trend of changes in KL divergence between the approximate and true posteriors. When the change is less than a threshold, the iterating process terminated and the parameters at the current iteration serve as the final estimates.

Hence, the distribution q(μk)~N(mk,σk2) that approximates q(µ|Y) to P(µ|Y) can be found, where σ_k contains the uncertainty information of the kth parameter estimated from a given scan. Instead of directly using the variance of each parameter, we use the coefficient of variation, which has the advantage of being independent of the intensity scale. It can simply be represented as a ratio between the standard deviation and the mean of a given parameter: CoV=σ/m, which is unitless. One can conclude that the higher the CoV, the more uncertainty there is associated with the parameter.

Experiments

Simulation studies

We simulated 3D (2D+t) CT perfusion scans using a modification of the brain phantom described in (26) (downloadable at https://www5.cs.fau.de/research/data/digital-brain-perfusion-phantom/index.html). We defined two annotated tissue classes to mimic a scan with ischemia: healthy tissue and tissue with reduced CBF and increased MTT (Figure 2A). Perfusion parameters in Table 1, in which was assigned according to (27), were assigned to the annotated tissues. From the perfusion phantom, we could generate dynamic projection data based on the forward model defined in Eqs. [8,9]. A total of 50 one-second frames were simulated, for which TACs are shown in Figure 2B. Data were generated by forward projecting each frame with an acquisition protocol (Table 2) based on a simulated scanner model. Poisson noise was added to the projections by setting the number of detected counts in a blank scan to 1e5. Noise-free and noisy projections were reconstructed to generate two sets of acquired dynamic images, from which perfusion parameters were inferred. In all simulations, reconstruction was performed with a helical Feldkamp-Davis-Kress algorithm (28,29). A total of four simulations were performed to evaluate the proposed method:

Figure 2 Simulation phantom setup. (A) CBF and MTT perfusion phantoms from which data were generated. Two lesions with reduced CBF and increased MTT are labeled with arrows. (B) Time attenuation curves of voxels in different regions. AIF is scaled by 0.1 for visual inspection. A region-of-interest for quantification was drawn inside the lesion as indicated by the red circle. CBF, cerebral blood function; MTT, mean transit time; AIF, arterial input function.

Comparison with conventional methods

In the first simulation, we compared the proposed method with conventional methods with ideal simulated measurements. The noise-free dynamic images were used to compute the perfusion maps with both methods, and the bias was assessed by comparing image profiles across the lesion in perfusion maps generated from the SVD-based methods, the proposed method, and the ground truth. Note that only this simulation used the noise-free images, the rest using noisy simulated images.

Effect of scan interval

In the second simulation, we evaluated scan protocols with different scan intervals using the proposed method. Here, scan interval indicates the time between the acquisition of the two frames. It was shown previously that scan interval variation could have an impact on perfusion imaging by introducing bias in the perfusion parameters (10,11). We varied the scan intervals by undersampling the simulated noisy dynamic frames. For a normal setup, the scan interval was zero; for long-interval setups, the scan intervals were either 1 s or 3 s, reducing the number of frames to 25 and 12, respectively. The hypothesis was that the Bayesian method can tolerate large intervals better than the SVD-based methods. Moreover, the measurement would be less reliable when fewer data were available, that is, the proposed Bayesian method would identify larger uncertainties in parameters. We also calculated the predicted enhanced TAC for a region-of-interest (ROI) in the lesion for Bayesian methods by applying the estimated distribution of the perfusion parameters to the forward model in Eq. [10]. This enhanced TAC was then compared with the true one.

Effect of truncation

In the third simulation, we evaluated the effect of truncated measurements on parameter inference with the proposed method. Truncation commonly occurs in CT perfusion when the acquisition is too short to capture the complete attenuation change due to poor cardiac output or large vessel occlusion (30). It has been shown that incomplete acquisition falsely reduces perfusion estimates and hence may alter patient management (12,13). To simulate such an effect, we made the last 20 or 30 s of the measurement unavailable for parameter inference. We hypothesized that the proposed Bayesian method would handle missing data better than the SVD-based method. Moreover, we hypothesized that the measurement would be less reliable because of the missing data, that is, the proposed Bayesian method would identify larger uncertainties in parameters. Again, we calculated the predicted enhanced TAC for an ROI in the lesion for the Bayesian method and compared it to the true one.

Effect of patient motion

In the fourth simulation, we evaluated the effect of motion on parameter inference with the proposed method. Involuntary patient motion can create inconsistencies in temporal measurements and hence have a negative impact on perfusion map estimation (31,32). We examined the effect of motion on the uncertainty of the parameter estimates.

For this investigation, we generated simulated motion-corrupted projections of the phantom by applying previously measured human motion (33) during forward projection and reconstructing the motion-corrupted images. Perfusion maps derived from the motion-free and motion-corrupted dynamic images were compared. The hypothesis was that the measurement would be less reliable in the motion case because of inconsistent data, that is, the proposed Bayesian method would identify larger uncertainties in parameters estimated from motion-corrupted data.

Patient study

The anonymized raw data of two patients who had previously undergone a head CT scan in the Department of Radiology at Westmead Hospital in Sydney, Australia, were collected with the approval of the Human Research Ethics Committee of the Western Sydney Local Health District. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). All human studies were performed with the approval of the Human Research Ethics Committee of the Western Sydney Local Health District (approval number 2019/ETH02047). The scans were acquired with a Definition Force scanner (Siemens Healthineers, Forchheim, Germany). Before each scan, 45 mL of iodine contrast agent was injected intravenously at a rate of 6 mL/s, followed by a 40 mL saline flush. After that, a total of 33 frames (26×1.5 s, 7×3 s) were acquired in axial (shuttle) mode. The first acquired frame was labeled as the baseline frame (mask). The reconstructions were initially performed with the vendor’s filtered back-projection algorithm with a smoothing kernel of H40f. The pixel dimensions of the reconstructed image were 0.4×0.4×5 mm³. DICOM volumes were exported without further post-processing.

One patient dataset was motion-free and the other presented visible motion artifacts. For the former, the proposed Bayesian method was directly applied to derive the perfusion parameters, while, for the dataset contaminated by the motion, we identified both inter- and intra-frame motion artifacts in the exported dynamic reconstructed frames. These images were labeled as NMC. Patient motion would be expected to result in greater uncertainty of parameter estimates, whereas a reduction in motion artifacts and blurring following the application of a motion correction algorithm would be expected to recover the reliability of the parameter estimates. As a test of the proposed method, we applied it to evaluate the effect of motion correction methods in terms of the perfusion parameters. Two motion correction methods were implemented in this study. The first method only compensated for inter-frame motion by directly operating on the exported images (34,35). All dynamic frames were rigidly registered to the first (baseline) frame to remove the inter-frame motion. The corrected images were labelled as MC1. The second method included an additional step before inter-frame motion compensation, which compensated for intra-frame motion using the method described in (36,37). This iterative data-driven method jointly estimated the six degrees-of-freedom rigid motion of the head at each projection and the motion-corrected reconstruction from the raw projection data. At each iteration, head pose estimates were updated analytically at multiple projection angles by adjusting the head pose to reduce the differences between the forward and measured projections. The motion-corrected image was updated by taking the estimated motion into account during reconstruction with a fully 3D maximum likelihood expectation maximization or Feldkamp-Davis-Kress algorithm. Motion correction was applied to each frame independently. This was followed by inter-frame motion correction as in the first correction method. The corrected images were labeled as MC2. After motion correction, perfusion analysis was performed by applying the proposed Bayesian inference algorithms to use NMC, MC1, and MC2 images.

In contrast to the simulation study, AIF was unknown here and was derived directly from the images as follows. We first applied a threshold to the summation of all frames and drew a 2D ROI over the internal carotid region. We then corrected the AIF for the partial volume effect due to the small diameter of the target arteries—a venous output function was derived similarly in a straight sinus region, and the AIF was scaled to have the same area under the curve as the venous output function to account for the underestimation (38). After inference, hematocrit correction was performed on CBV and CBF [multiplying with a constant 0.733 (15)] to account for the different hematocrit values between arteries and capillaries.

Quantification

All generated perfusion parameter images were compared with the truth in Figure 2A by calculating the mean squared error and the bias. Denoting IM0 as the true image and IM1 as the image for comparison, we calculated the mean squared error MSE:

MSE=1N∑jN(IM1j-IM0j)2[23]

where N is the total number of voxels, and j is the index of the voxel. The bias between two images was represented as:

bias(%)=|mean(IM1ROI)−mean(IM0ROI)mean(IM0ROI)|×100%[24]

This calculation was only performed in an ROI drawn inside the lesion, as shown in Figure 2A. The CBF and MTT intensity distributions were also calculated for that ROI, which were also Gaussian distributed. The mean of the distribution was the mean at all voxels in the ROI; the standard deviation of the distribution was the mean of the square root of the variance at all voxels in the ROI. The wider the distribution, the greater the uncertainty of the estimated parameters in the ROI.

Results

Simulation studies

Comparison with conventional methods

Figure 3A shows the accuracy (MSE and bias) of CBF and MTT maps derived with the proposed method and with the deconvolution methods. Bias was calculated from the voxels of an ROI within a lesion (Figure 2A). Figure 3B shows profiles from the MTT images in the lower row of Figure 3A, at the location indicated by the red dotted line. Compared with the SVD-based method, the proposed Bayesian method resulted in lower MSE and bias for CBF and MTT, and, therefore also, lower bias for CBV. The MTT image profile obtained with the Bayesian method corresponded well to the reference profile.

Figure 3 A comparison of the results with SVD and the proposed method without noise in Simulation 1. (A) Compared with the true images (REF), the proposed method shows markedly less bias in lesions than the SVD-based methods for both CBF and MTT. (B) Profile plots (red dotted line in a) of MTT images show the reduction in bias with the proposed method.

Effect of scan interval

In the second simulation, we assessed the impact of varying the scan interval on the MSE and bias of parameter estimates. As shown in Figure 4A, the SVD-based methods resulted in markedly increased MSE and bias at large scan intervals, while the proposed method was much less affected by increased scan interval. At larger scan intervals, we observed that the mean CBF and MTT maps were noisier than those obtained with normal sampling (Figure 4). The CoV maps and the CBF distribution plots indicated differences in CBF mean and variance with different scan intervals. The number of iterations used for full data, and scan intervals of 1 s and 3 s were 6, 8 and 12, respectively. We also calculated the predicted enhanced TAC sampled in a lesion ROI for all methods in Figure 5.

Figure 4 The effect of varying the scan interval in Simulation 2. For the full data, the SNR of the estimated CBF map (from the proposed method) was 3.2. (A) Mean CBF estimate and its CoV for three different scan intervals. The CBF and CoV values calculated with the SVD-based method have larger bias than those calculated with the Bayesian method, and bias increases with the interval. (B) The CBF intensity distribution of a sampled ROI (yellow circle in A) with different intervals. The wider the distribution, the greater the uncertainty in the parameters. Similar results were obtained for MTT (not shown). SNR, signal to noise; CBF, cerebral blood flow; CoV, coefficient of variation; SVD, singular value decomposition; MTT, mean transit time.

Figure 5 Time-attenuation curves for the lesion ROI from (A) simulation 2 and (B) simulation 3. With the inferred Gaussian distributed CBF and MTT (mean ± SD), we can generate the lower and upper bounds of the time-attenuation curves. The predicted enhanced TAC with the proposed method and true enhanced TAC are well correlated. With more deterioration in data, the predicted enhanced TAC is more biased and with more uncertainty. HU, Hounsfield unit.

Effect of truncation

In the third simulation, our observations on the effect of truncation in Figure 6 were similar to those regarding scan interval in Figure 4. MSE and bias in MTT estimates were markedly degraded as truncation was increased with the SVD-based method while the proposed method was relatively unaffected. With severe truncation, the MTT estimate was biased and its variance higher, as indicated by the CoV images, compared to the case with full measurement data. The number of iterations of the Bayesian method for the full data and truncations of last 25 s and 30 s of data were 6, 6 and 9, respectively. Figure 5B shows the predicted enhanced TAC sampled in a lesion ROI for all methods.

Figure 6 The effect of truncating the measured data in Simulation 3. For the full data, the SNR of the estimated MTT map (from the proposed method) was 2.8. (A) Mean MTT estimate and CoV for full data and two different degrees of truncation. The MTT and CoV values calculated with the SVD-based method have larger bias, and the bias increases with the truncation. (B) The MTT intensity distribution of a sampled white matter region with different amounts of missing data. Similar results were obtained for CBF (not shown). SNR, signal to noise; CBF, cerebral blood flow; CoV, coefficient of variation; SVD, singular value decomposition; MTT, mean transit time.

Effect of patient motion

In the fourth simulation with simulated motion, we observed that the mean CBF and MTT map had noticeable visual artifacts (Figure 7). The CBF_CoV images indicated that motion made the scan less reliable. High CoV values corresponded to regions with artifacts in the mean estimate maps, and CoV values were generally lower in the motion-free image. This suggests that CoV mapping is useful for identifying parts of the brain where parameter estimates are problematic due to motion-induced uncertainty. The numbers of iterations of the Bayesian method used for MC and NMC images were 6 and 7, respectively.

Figure 7 The effect of motion in Simulation 4. (A) Mean CBF estimate and CoV with and without simulated motion. Severe motion creates increased uncertainty in CBF estimates, as revealed in the CoV map. (B) The CBF intensity distribution of a sampled white matter region. Similar results were obtained for MTT (not shown). (C) A segment of realistic rigid motion that is applied. Note that only in-plane motion is simulated (one rotation, two translations) in this 2D+t simulation. CBF, cerebral blood flow; CoV, coefficient of variation.

Patient study

For the motion-free dataset, the proposed Bayesian method was directly applied to derive the perfusion parameters. According to the clinical report, this patient had an ischemic lesion due to a left middle cerebral artery infarct with increased MTT and a corresponding drop in CBF. Figure 8 shows the CBF, MTT and CBV maps inferred using our method. The lesion and the normal tissue had clearly different hemodynamic parameters, indicated by the arrows in Figure 8. The average CBF value within a region in the lesion was 13.2±4.1 mL/100 g/min, while for a region in the normal tissue it was 22.6±6.2 mL/100 g/min. The average MTT value in the lesion was 22.3±4.8 s, while in the normal tissue it was 11.4±3.4 s. The average CBV values in the lesion and normal tissues were comparable. In addition, the uncertainty of the estimates can be assessed alongside the mean estimates as indicated in the middle column of Figure 8.

Figure 8 Mean and CoV maps inferred from a patient scan using the proposed method, in which a lesion due to the left middle cerebral infarct was indicated by the arrow. The average CBF value in the lesion was over 40% lower than the normal tissue and the average MTT value in the lesion was over 96% higher than the normal tissue. Nine iterations of the Bayesian method were performed.

For the dataset with motion artifacts, MC1 and MC2 were generated as described in the Methods section. The goal was to apply the proposed method to evaluate the effectiveness of two motion correction methods in terms of recovering the perfusion parameters. For both CBF and MTT, mean estimates and CoV were derived with MC1, MC2, and noMC. In Figure 9 severe artifacts are visible with noMC. With MC1, the quality of the dynamic scan was successfully recovered, and the inferred parameters were more reliable than those obtained with noMC, although there were still some residual artifacts. With MC2, the images were further improved, and the reliability of the estimation was the highest among the three. The results correspond well to the hypothesis of progressive recovery of the correction methods (from noMC to MC1 to MC2).

Figure 9 Mean and CoV maps of CBF inferred from a motion-corrected clinical CT perfusion scan using the proposed method. (A) MC2 yields fewer visible artifacts and lower overall CoV than MC1 or noMC. (B) CBF distributions of a selected ROI (yellow circle in A). MC2 provided the narrowest distribution, indicating that it provides the lowest uncertainty. Six, eight and twelve iterations of the Bayesian method were used for noMC, MC1 and MC2, respectively. Similar results were obtained for MTT (not shown).

Discussion

CT perfusion imaging results can be affected by many factors in practice, including the selection of scanner acquisition protocols, the choice of the post-processing technique, data truncation, and data perturbance caused, for example, by patient movement. All may affect the reliability of hemodynamic parameter estimation, which is crucial for accurate the diagnosis and appropriate treatment of stroke patients. Conventional inference methods cannot tolerate such imperfection well and cannot assess the reliability of the data, although such information could be useful diagnostically. We proposed a Bayesian method that can calculate CT perfusion maps and estimate the uncertainty at each voxel. The proposed method assumes better residual function hence the underlying perfusion model. Compared with conventional SVD-based methods, this approach has the advantages of (I) estimating parameters with less bias than conventional methods, especially when the measurements are imperfect and (II) providing the probability of an estimate, which can be used to guide the design of scan protocols or to compare the performance of different data processing or artifact reduction methods. We applied the proposed method to three simulated abnormal scan scenarios involving increased frame intervals, truncated measurement, and patient movement. The results indicated that the uncertainty measure provided by the method could allow scan protocols of varying reliability to be differentiated. In one patient study, the estimate CBF and MTT maps indicated an ischemic lesion. In a second patient study, we evaluated the proposed method on a stroke scan in which the patient exhibited severe motion. Two motion correction techniques were compared in terms of how much they improved the scan reliability, as assessed by the proposed Bayesian method.

The running time of the proposed method was longer than the conventional SVD-based methods, as Bayesian inference is an iterative process. For example, to produce the results in Figure 3, SVD took 2 s per slice while Tikhonov-SVD took 6 s and the Bayesian method took 10 s. Although the running time was longer for the Bayesian method, more accurate estimates, and additional information about the distribution of the parameters of interest were obtained.

The proposed method successfully provided mean and uncertainty estimates for CBF and MTT. However, it is also possible to infer other parameters, such as bolus arrival time delay and permeability, if the associated acquisition model is assumed. For example, a delay exists between the injection time and the contrast agent arrival time in the tissue of interest, which has been shown to affect the CBF estimates, especially for negative time shifts and at high SNR (17). To include this effect, the convolution model in Eq. [1] could be modified, and the inference would proceed with minor change in future work. Permeability is an important indicator of the degree of blood-brain barrier damage. The Patlak model is commonly used to infer this parameter and could be incorporated with the proposed Bayesian method. However, a longer scan time (~200 s) is usually needed to provide enough information for the Patlak model, and such a scan protocol is not available at this time in our center.

Another possible way to infer the distribution of perfusion parameters is the frequentist method, which often requires repeatable measurements (which is not feasible for patients) or simulated measurements using techniques such as bootstrapping. We chose a Bayesian method in this study instead since we believe it is more appropriate for our application where the measured data are fixed and parameters are random. In addition, the Bayesian method has the advantage that prior knowledge can easily be included in the Bayesian inference process, whereas it is not obvious how to do this for bootstrapping (39). In the future work, spatial-temporal regularization could be included in the proposed method to reduce speckle noise in estimation and enable CT perfusion scans to be performed with lower tube current and radiation dose. A Variational Bayesian approach was chosen in this study. As mentioned earlier in the Introduction, MCMC is another popular Bayesian method shown to be efficient in inferring kinetic parameters in dynamic contrast-enhanced cardiac MR applications (40). MCMC can also incorporate the spatial prior information to reduce the speckle noise in estimation map (41).

We inferred the parameter distributions from dynamic reconstructed images. In principle, it is possible to infer the desired parameters directly from the measured raw projections without reconstruction. Applications of direct kinetic parameters estimation exist in CT (42), PET (43,44), and MRI (40,45,46). However, these approaches can only provide maximum likelihood estimates of parameters, not the full distribution. Another concern in CT perfusion is that it is difficult to separate the coupling between the attenuation and kinetic distribution since a perfect forward model can often not be guaranteed. For example, if there is patient motion during a scan, simultaneous motion correction and perfusion map inference can be problematic because errors from one side will likely propagate to the other.

A limitation of our measurement model is that the observation noise is assumed to be an independently and identically distributed (i.i.d.) Gaussian variable. Within each frame, the noise level is scaled according to the duration of that frame. This assumption leads to great convenience for computation. If the CT scan dose is very low, the actual noise model can be complicated (47,48) and the one assumed here may yield poor inference results. Other inference techniques may be required to consider the correct noise model. For example, techniques approximating the Poisson to Gaussian distribution may be required before applying Bayesian inference. Another limitation within the proposed measurement model is the bolus arrival delay was not considered and hence inferred. As stated earlier, the model can be modified to account for this effect. A further clinically relevant limitation is that we did not evaluate the estimated parameter maps for identifying the infarct core and penumbra region, which could be relevant to a clinical decision. A more detailed investigation of the clinical impact of the method is needed.

Conclusions

We proposed and applied a Bayesian approach that infers the hemodynamic parameters in CT perfusion imaging. The method can yield mean parameter estimates and uncertainty estimates of such parameters under imperfect scan scenarios. It can be used to assess the confidence of parameter estimates under imperfect scan conditions to evaluate protocol designs. Further studies are required to fully evaluate the clinical utility of the proposed method.

Acknowledgments

The authors wish to thank Johan Nuyts from KU Leuven for the discussion; Krystal Moore from the Department of Radiology at Westmead Hospital for assistance with data collection; and Karl Stierstorfer, Siemens Healthcare, Forchheim, Germany, for his help with reading the raw CT data.

Funding: This work was partly funded by Key Laboratory for Magnetic Resonance and Multimodality Imaging of Guangdong Province (2020B1212060051).

Footnote

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://dx.doi.org/10.21037/qims-21-338). Dr. DL serves as an unpaid editorial board member of Quantitative Imaging in Medicine and Surgery. The other authors have no conflict of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). This study was performed with the approval of the Human Research Ethics Committee of the Western Sydney Local Health District (approval number 2019/ETH02047).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Yang Q, Tong X, Schieb L, Vaughan A, Gillespie C, Wiltz JL, King SC, Odom E, Merritt R, Hong Y, George MG. Vital Signs: Recent Trends in Stroke Death Rates - United States, 2000-2015. MMWR Morb Mortal Wkly Rep 2017;66:933-9. [Crossref] [PubMed]
2. Benjamin EJ, Blaha MJ, Chiuve SE, Cushman M, Das SR, Deo R, et al. Heart Disease and Stroke Statistics-2017 Update: A Report From the American Heart Association. Circulation 2017;135:e146-603. [Crossref] [PubMed]
3. Wintermark M, Flanders AE, Velthuis B, Meuli R, van Leeuwen M, Goldsher D, et al. Perfusion-CT assessment of infarct core and penumbra: receiver operating characteristic curve analysis in 130 patients suspected of acute hemispheric stroke. Stroke 2006;37:979-85. [Crossref] [PubMed]
4. Meier P, Zierler KL. On the theory of the indicator-dilution method for measurement of blood flow and volume. J Appl Physiol 1954;6:731-44. [Crossref] [PubMed]
5. Riordan AJ, Prokop M, Viergever MA, Dankbaar JW, Smit EJ, de Jong HW. Validation of CT brain perfusion methods using a realistic dynamic head phantom. Med Phys 2011;38:3212-21. [Crossref] [PubMed]
6. Ibaraki M, Ohmura T, Matsubara K, Kinoshita T. Reliability of CT perfusion-derived CBF in relation to hemodynamic compromise in patients with cerebrovascular steno-occlusive disease: a comparative study with 15O PET. J Cereb Blood Flow Metab 2015;35:1280-8. [Crossref] [PubMed]
7. Calamante F, Gadian DG, Connelly A. Quantification of bolus-tracking MRI: Improved characterization of the tissue residue function using Tikhonov regularization. Magn Reson Med 2003;50:1237-47. [Crossref] [PubMed]
8. Frindel C, Robini MC, Rousseau D A. 3-D spatio-temporal deconvolution approach for MR perfusion in the brain. Med Image Anal 2014;18:144-60. [Crossref] [PubMed]
9. Giacalone M, Frindel C, Robini M, Cervenansky F, Grenier E, Rousseau D. Robustness of spatio-temporal regularization in perfusion MRI deconvolution: An application to acute ischemic stroke. Magn Reson Med 2017;78:1981-90. [Crossref] [PubMed]
10. Wintermark M, Smith WS, Ko NU, Quist M, Schnyder P, Dillon WP. Dynamic perfusion CT: optimizing the temporal resolution and contrast volume for calculation of perfusion CT parameters in stroke patients. AJNR Am J Neuroradiol 2004;25:720-9. [PubMed]
11. van Ommen F, Kauw F, Bennink E, Dankbaar JW, Viergever MA, de Jong HWAM. Effect of prolonged acquisition intervals for CT-perfusion analysis methods in patients with ischemic stroke. Med Phys 2019;46:3156-64. [Crossref] [PubMed]
12. Borst J, Marquering HA, Beenen LF, Berkhemer OA, Dankbaar JW, Riordan AJ, Majoie CB. MR CLEAN investigators. Effect of extended CT perfusion acquisition time on ischemic core and penumbra volume estimation in patients with acute ischemic stroke due to a large vessel occlusion. PLoS One 2015;10:e0119409 [Crossref] [PubMed]
13. Copen WA, Deipolyi AR, Schaefer PW, Schwamm LH, González RG, Wu O. Exposing hidden truncation-related errors in acute stroke perfusion imaging. AJNR Am J Neuroradiol 2015;36:638-45. [Crossref] [PubMed]
14. Barrett HH, Myers KJ. Foundations of image science. Hoboken, NJ: Wiley-Interscience, 2004.
15. Kudo K, Terae S, Katoh C, Oka M, Shiga T, Tamaki N, Miyasaka K. Quantitative cerebral blood flow measurement with dynamic perfusion CT using the vascular-pixel elimination method: comparison with H2(15)O positron emission tomography. AJNR Am J Neuroradiol 2003;24:419-26. [PubMed]
16. Fieselmann A, Kowarschik M, Ganguly A, Hornegger J, Fahrig R. Deconvolution-Based CT and MR Brain Perfusion Measurement: Theoretical Model Revisited and Practical Implementation Details. Int J Biomed Imaging 2011;2011:467563 [Crossref] [PubMed]
17. Ostergaard L, Weisskoff RM, Chesler DA, Gyldensted C, Rosen BR. High resolution measurement of cerebral blood flow using intravascular tracer bolus passages. Part I: Mathematical approach and statistical analysis. Magn Reson Med 1996;36:715-25. [Crossref] [PubMed]
18. Hansen PC. Regularization tools: A Matlab package for analysis and solution of discrete ill-posed problems. Numer Algorithms 1994;6:1-35. [Crossref]
19. Bredno J, Olszewski ME, Wintermark M. Simulation model for contrast agent dynamics in brain perfusion scans. Magn Reson Med 2010;64:280-90. [Crossref] [PubMed]
20. Fieselmann A, Ganguly A, Deuerling-Zheng Y, Zellerhoff M, Rohkohl C, Boese J, Hornegger J, Fahrig R. Interventional 4-D C-arm CT perfusion imaging using interleaved scanning and partial reconstruction interpolation. IEEE Trans Med Imaging 2012;31:892-906. [Crossref] [PubMed]
21. Pisana F, Henzler T, Schönberg S, Klotz E, Schmidt B, Kachelrieß M. Noise reduction and functional maps image quality improvement in dynamic CT perfusion using a new k-means clustering guided bilateral filter (KMGB). Med Phys 2017;44:3464-82. [Crossref] [PubMed]
22. Bennink E, Riordan AJ, Horsch AD, Dankbaar JW, Velthuis BK, de Jong HW. A fast nonlinear regression method for estimating permeability in CT perfusion imaging. J Cereb Blood Flow Metab 2013;33:1743-51. [Crossref] [PubMed]
23. Robert CP, Casella G. Monte Carlo statistical methods. Berlin: Springer, 2004.
24. Blei DM, Kucukelbir A, McAuliffe JD. Variational Inference: A Review for Statisticians. J Am Stat Assoc 2017;112:859-77. [Crossref]
25. Gelman A, Carlin JB, Stern HS, et al. Bayesian data analysis. Boca Raton, FL: CRC Press, 2013.
26. Manhart MT, Kowarschik M, Fieselmann A, Deuerling-Zheng Y, Royalty K, Maier AK, Hornegger J. Dynamic iterative reconstruction for interventional 4-D C-arm CT perfusion imaging. IEEE Trans Med Imaging 2013;32:1336-48. [Crossref] [PubMed]
27. Parkes LM, Rashid W, Chard DT, Tofts PS. Normal cerebral perfusion measurements using arterial spin labeling: reproducibility, stability, and age and gender effects. Magn Reson Med 2004;51:736-43. [Crossref] [PubMed]
28. Danielsson P, Edholm P, Eriksson J, et al. Towards exact reconstruction for helical cone-beam scanning of long objects. A new detector arrangement and a new completeness condition. Pittsburgh, PA: 1997 Meeting on Fully 3D Image Reconstruction in Radiology and Nuclear Medicine, 1997:141-4.
29. Tam KC, Samarasekera S, Sauer F. Exact cone beam CT with a spiral scan. Phys Med Biol 1998;43:1015-24. [Crossref] [PubMed]
30. Campbell BC, Christensen S, Levi CR, Desmond PM, Donnan GA, Davis SM, Parsons MW. Cerebral blood flow is the optimal CT perfusion parameter for assessing infarct core. Stroke 2011;42:3435-40. [Crossref] [PubMed]
31. Konstas AA, Goldmakher GV, Lee TY, Lev MH. Theoretic basis and technical implementations of CT perfusion in acute ischemic stroke, part 1: Theoretic basis. AJNR Am J Neuroradiol 2009;30:662-8. [Crossref] [PubMed]
32. Fahmi F, Beenen LF, Streekstra GJ, Janssen NY, de Jong HW, Riordan A, Roos YB, Majoie CB, Vanbavel E, Marquering HA. Head movement during CT brain perfusion acquisition of patients with suspected acute ischemic stroke. Eur J Radiol 2013;82:2334-41. [Crossref] [PubMed]
33. Kim JH, Sun T, Alcheikh AR, Kuncic Z, Nuyts J, Fulton R. Correction for human head motion in helical x-ray CT. Phys Med Biol 2016;61:1416-38. [Crossref] [PubMed]
34. de Lucas EM, Sánchez E, Gutiérrez A, Mandly AG, Ruiz E, Flórez AF, Izquierdo J, Arnáiz J, Piedra T, Valle N, Bañales I, Quintana F. CT protocol for acute stroke: tips and tricks for general radiologists. Radiographics 2008;28:1673-87. [Crossref] [PubMed]
35. Fahmi F, Marquering HA, Borst J, Streekstra GJ, Beenen LF, Niesten JM, Velthuis BK, Majoie CB, vanBavel E. DUST study. 3D movement correction of CT brain perfusion image data of patients with acute ischemic stroke. Neuroradiology 2014;56:445-52. [Crossref] [PubMed]
36. Sun T, Kim JH, Fulton R, Nuyts J. An iterative projection-based motion estimation and compensation scheme for head x-ray CT. Med Phys 2016;43:5705. [Crossref] [PubMed]
37. Sun T, Jacobs R, Pauwels R, Tijskens E, Fulton R, Nuyts J. A motion correction approach for oral and maxillofacial cone-beam CT imaging. Phys Med Biol 2021;66:125008 [Crossref] [PubMed]
38. van der Schaaf I, Vonken EJ, Waaijer A, Velthuis B, Quist M, van Osch T. Influence of partial volume on venous output and arterial input function. AJNR Am J Neuroradiol 2006;27:46-50. [PubMed]
39. Alfaro ME, Zoller S, Lutzoni F. Bayes or bootstrap? A simulation study comparing the performance of Bayesian Markov chain Monte Carlo sampling and bootstrapping in assessing phylogenetic confidence. Mol Biol Evol 2003;20:255-66. [Crossref] [PubMed]
40. Dikaios N. Stochastic Gradient Langevin dynamics for joint parameterization of tracer kinetic models, input functions, and T1 relaxation-times from undersampled k-space DCE-MRI. Med Image Anal 2020;62:101690 [Crossref] [PubMed]
41. Scannell CM, Chiribiri A, Villa ADM, Breeuwer M, Lee J. Hierarchical Bayesian myocardial perfusion quantification. Med Image Anal 2020;60:101611 [Crossref] [PubMed]
42. Fang R, Gupta A, Sanelli PC. Direct estimation of permeability maps for low-dose CT perfusion. 2016 IEEE 13th International Symposium on Biomedical Imaging (ISBI) 2016:739-42.
43. Kamasak ME, Bouman CA, Morris ED, Sauer K. Direct reconstruction of kinetic parameter images from dynamic PET data. IEEE Trans Med Imaging 2005;24:636-50. [Crossref] [PubMed]
44. Wang G, Qi J. Direct estimation of kinetic parametric images for dynamic PET. Theranostics 2013;3:802-15. [Crossref] [PubMed]
45. Dikaios N, Arridge S, Hamy V, Punwani S, Atkinson D. Direct parametric reconstruction from undersampled (k, t)-space data in dynamic contrast enhanced MRI. Med Image Anal 2014;18:989-1001. [Crossref] [PubMed]
46. Guo Y, Lingala SG, Zhu Y, Lebel RM, Nayak KS. Direct estimation of tracer-kinetic parameter maps from highly undersampled brain dynamic contrast enhanced MRI. Magn Reson Med 2017;78:1566-78. [Crossref] [PubMed]
47. Whiting BR, Massoumzadeh P, Earl OA, O'Sullivan JA, Snyder DL, Williamson JF. Properties of preprocessed sinogram data in x-ray computed tomography. Med Phys 2006;33:3290-303. [Crossref] [PubMed]
48. Ma J, Liang Z, Fan Y, Liu Y, Huang J, Chen W, Lu H. Variance analysis of x-ray CT sinograms in the presence of electronic noise background. Med Phys 2012;39:4051-65. [Crossref] [PubMed]