A Novel Statistical Optimization Algorithm for Estimating Perfusion Curves in Susceptibility Contrast-Enhanced MRI

Dynamic susceptibility contrast-enhanced magnetic resonance imaging is an important tool for evaluating intravascular indicator dynamics, which in turn is valuable for understanding brain physiology and pathophysiology. This procedure usually involves fitting a gamma-variate function to observed concentration-time curves in order to eliminate undesired effects of recirculation and the leakage of contrast agents. Several conventional curve-fitting approaches are routinely applied. The nonlinear optimization methods typically are computationally expensive and require reliable initial values to guarantee success, whereas a logarithmic linear least-squares (LL-LS) method is more stable and efficient, and does not suffer from the initial-value problem, but it can show degraded performance, especially when a few data or outliers are present. In this paper, we demonstrate, that the original perfusion curve-fitting problem can be transformed into a gamma-distribution-fitting problem by treating the concentration-time curves as a random sample from a gamma distribution with time as the random variable. A robust maximum-likelihood estimation (MLE) algorithm can then be readily adopted to solve this problem. The performance of the proposed method is compared with the nonlinear Levenberg-Marquardt (L-M) method and the LL-LS method using both synthetic and real data. The results show that the performance of the proposed approach is far superior to those of the other two methods, while keeping the advantages of the LL-LS method, such as easy implementation, low computational load, and dispensing with the need to guess the initial values. We argue that the proposed method represents an attractive alternative option for assessing intravascular indicator dynamics in clinical applications. Moreover, we also provide valuable suggestions on how to select valid data points and set the initial values in the two traditional approaches (LL-LS and nonlinear L-M methods) to achieve more reliable estimations.

1. Introduction

Dynamic susceptibility contrast-enhanced magnetic resonance imaging (DSC-MRI) is a useful tool for the quantitative assessment of perfusion-related cerebrovascular parameters in clinical applications (Kosior and Frayne, 2010; Emblem et al., 2014; Arzanforoosh et al., 2021; Jin and Cho, 2021). It requires the injection of a paramagnetic contrast agent. A bolus of the intravascular contrast agent passing through the tissue of interest produces local magnetic-field inhomogeneities that lead to a reduction in the transverse relaxation time (T2*) of the bulk tissue (Rosen et al., 1991; Li, 2014). This susceptibility effect is then recorded in a series of rapidly obtained T2*-weighted gradient-echo images. Converting the signal-time curves into concentration-time (i.e., ΔR2^*(t)) curves will yield vascular information associated with tissue perfusion, such as the blood volume, blood flow, and mean transit time (MTT).

However, the concentration-time curve is inevitably contaminated by recirculation of the bolus into the region of interest and residual contrast agent in capillaries. A gamma-variate function has been used to model the first pass of the ΔR2^*(t) curve to eliminate these undesired effects (Norman et al., 1981; Davenport, 1983). The original motivation for using this function as an empirical model was its heuristic resemblance to the expected relaxation-rate time-course during the first bolus passage (Thompson et al., 1964; Axel, 1980). The underlying physical connection of the dilution process with the gamma-variate expression has also been clarified in recent studies (Davenport, 1983; Mischi et al., 2008). The selection of an appropriate fitting method is therefore essential when analyzing intravascular indicator dynamics (Perkiö et al., 2002; Pianykh, 2012; Romain et al., 2017; Quarles et al., 2019).

The nonlinear Levenberg-Marquardt (L-M) optimization method is typically used to solve the fitting problem of gamma-variate function since it provides superior accuracy (Benner et al., 1997; Li et al., 2003). However, the use of a nonlinear method requires reliable initial values to guarantee convergence: a moderate deviation from the true values may cause the fit to diverge, and in extreme cases, these errors might be even higher when applying numerical integration over the sample points (Benner et al., 1997). Moreover, nonlinear methods usually suffer from a high computational load. These reasons explain why nonlinear fitting methods have not played a prominent role in clinical applications. A logarithmic linear least-squares (LL-LS) approach (Madsen, 1992; Chan and Nelson, 2004), on the other hand, is more stable and efficient as well as not requiring the initial values to be guessed. The LL-LS method uses a logarithm operation to transforms the original gamma-variate function into an equivalent, linear regression form. Standard linear least-squares fitting is then applied. However, the logarithm operation complicates the statistics of the data, and means that the noise no longer conforms to a Gaussian distribution. As a result, the LS estimator will be biased and the quality of the fitting will be subject to limits imposed by the statistical nature of the data, especially when there are few data or outliers are present.

In a recent study, we needed to determine the regional cerebral blood volume (CBV) on a finer voxel scale throughout the brain in order to produce a more accurate hemodynamic assimilation (Hu et al., 2012; Zhang et al., 2016). Only five or six valid data points were collected in the first-pass phase for perfusion-curve fitting due to the sampling interval being longer than usual. Since the existing two methods were found to perform poorly, it was found to be necessary to develop an alternative approach that combines the advantages of the two methods and still exhibits satisfactory performance even when there are few data points. In this paper, we present a novel statistical optimization method that is well suited to such a real-world problem. We demonstrate that the original perfusion curve-fitting problem can be transformed into a gamma-distribution-fitting problem, by treating the concentration-time curves as a random sample from a gamma distribution with time as the random variable. A robust maximum-likelihood estimation (MLE) algorithm can then be readily applied to solve this problem. The proposed method was evaluated in experiments using both synthetic and real data. Our new method yields physiological information on cerebrovascular dynamics that are more stable and accurate than those obtained using a nonlinear method, while keeping the advantages of the LL-LS method, such as low computational load, and without requiring the initial values to be guessed.

2. Derivation of the Relative CBV

DSC MR measures the regional CBV by analyzing changes in the signal intensity after the first pass of a paramagnetic contrast agent. A bolus of intravascular paramagnetic contrast agent passing through the tissue of interest, produces local magnetic-field inhomogeneities that lead to a reduction in the transverse relaxation time T2* of the bulk tissue (Figure 1A) (Rosen et al., 1991). This susceptibility effect is then recorded in a series of rapidly obtained T2*-weighted gradient-echo images. There is an exponential relationship between the relative signal attenuation (S(t)/S₀) and the local tissue concentration of contrast agent (C_i(t)) (Rempp et al., 1994):

$\begin{array}{ll}{C}_i(t)=\frac{-k}{TE}·ln\left(\frac{S(t)}{S_0}\right)& (1)\end{array}$

where k is an unknown proportionality factor, TE is the echo time of the imaging sequences, S(t) is the MRI signal intensity at time point t, and S₀ is the baseline MRI signal intensity before administering the contrast agent.

FIGURE 1

Figure 1. (A) The original DSC-MRI signal (S(t)) obtained from one voxel after administering contrast agent. The baseline signal S₀, the first pass bolus and the second pass bolus are illustrated by the arrow, respectively. (B) Concentration-time curves (ΔR2^*) for the same voxel. The ΔR2^* signal measured during the first bolus pass are shown as filled circles, and other sample points are indicated hollow circles. The red area under the fitted ΔR2^* curves (red line) is the relative CBV (rCBV) estimated. a.u., arbitrary unit.

This relationship can be used to convert the observed signal-intensity-vs.-time curves (S(t)) into concentration-time curves (ΔR2^*(t)) on a voxel-by-voxel basis (Figure 1B). According to the indicator dilution theory, the relative cerebral blood volume (rCBV) is proportional to the area under the ΔR2^*(t) curve (Rempp et al., 1994):

$\begin{array}{ll}rCBV=\frac{k}{k_r}·\frac{k_H}{\rho }·\frac{\int ln(S(t)/S_0)dt}{\int ln(S_r(t)/S_{r_0})dt}·\frac{TE}{T{E}_r}& (2)\end{array}$

where S_r(t)/S_r0 is the relative signal reduction in the reference voxel, ρ = 1.04g/mL is the density of brain tissue, and correction factor k_H accounts for the difference in hematocrit between the voxel of interest and the reference voxel. For our calculations of relative concentration values, these correction factors were set as 1.0.

In practice, the concentration-time curves usually do not only reflect the first bolus pass, but are also affected by a second bolus pass starting 15–20 s after the first one due to the recirculation of the contrast agent bolus. Additionally, residual contrast agent often produces a gradually decreasing bias of the concentration-time curve that is not as steep as the assumed hypothetical exponential decay. For these reasons, the original concentration-time curve can not be used directly for the calculation of vascular parameters. The curve of interest is usually modeled with a gamma-variate function to eliminate these undesired effects (Norman et al., 1981; Davenport, 1983):

$\begin{array}{ll}\Delta R{2}^{*}(t)=y=A{(t-t_0)}^{\alpha }{e}^{-(t-t_0)/\beta }& (3)\end{array}$

where t₀ is the time when the contrast agent is applied in the specified area, and A, α, and β are parameters that determine the shape of the function. The rCBV value is proportional to the area (F) under the concentration-time curve, and can be calculated as follows (see Appendix A):

$\begin{array}{ll}F={\displaystyle {\int }_{t_0}^{\infty }}{C}_i(t)dt=A·{\beta }^{\alpha +1}·\Gamma (\alpha +1)& (4)\end{array}$

3. Nonlinear Levenberg-Marquardt Optimization Method

The Levenberg-Marquardt algorithm combines two numerical minimization algorithms: the gradient descent method and the Gauss-Newton method. Specifically, to fit a mode ŷ(t|p) of an independent variable t and a vector of n parameters p to a set of m data points (t_i, y_i), the estimation of parameters p is customary to minimize the sum of the weighted residuals between the known data y_i and the estimated fitting function ŷ(t|p).

$\begin{array}{ll}{\chi }^2(p)={\displaystyle \sum _{i=1}^m}[\frac{y(t_i-ŷ(t|p))}{{\sigma }_{y_i}}{]}^2={(y-\widehat{y}(p))}^T·W·(y-\widehat{y}(p)& (5)\end{array}$

Where σ_yi denotes the measurement error of y(t_i). The weighting matrix W is diagonal with W_ii = 1/σ_yi, or formally, it can be set to the inverse of measurement error covariance matrix. More generally, W can also be treated as an optimization parameter. This goodness-of-fit measure is called the chi-squared error criterion (Gavin, 2019). Immediately, the Levenberg-Marquardt algorithm adaptively varies the parameter updates between the gradient descent update and the Gauss-Newton update,

$\begin{array}{ll}\left[W^T·T·T+\lambda ·I\right]=J^T·W·(y-\widehat{y})& (6)\end{array}$

where J is the Jacobian matrix [∂ŷ/∂p], λ is damping parameter, the parameter update is h. Small values lead to a Gauss-Newton update and large values result in a gradient descent update. In Marquardt's update research, the value of λ are normalized to the values of W^T · T · T, thus,

$\begin{array}{ll}\left[W^T·T·T+\lambda ·diag(W^T·T·T)\right]=J^T·W·(y-\widehat{y})& (7)\end{array}$

and many variations of the Levenberg-Marquardt method also have been developed (Gavin, 2019). In this paper, the nonlinear Levenberg-Marquardt algorithm was provided by MATLAB's own toolbox.

4. Logarithmic Linear Least-Squares Method

Without loss of generality, we assume t₀ = 0. Equation (3) can then be simplified as

$\begin{array}{ll}y(t)=A{t}^{\alpha }{e}^{-t/\beta }& (8)\end{array}$

To find time t_max at which the bolus signal peaks, we takes the first derivative of Equation (8) and set it to 0:

$\begin{array}{l}{y}^{\prime }(t_{max})=A{t}_{max}^{\alpha -1}{e}^{-t_{max}/\beta }\left(\alpha -\frac{t_{max}}{\beta }\right)\\ =0& (9)\end{array}$

which yields

$\begin{array}{ll}{t}_{max}=\alpha \beta \Rightarrow \beta =\frac{t_{max}}{\alpha }.& (10)\end{array}$

Substituting β into Equation (8) and letting t = t_max leads to,

$\begin{array}{ll}{y}_{max}=A{\left(\frac{t_{max}}{e}\right)}^{\alpha }\Rightarrow A=y_{max}{\left(\frac{e}{t_{max}}\right)}^{\alpha }& (11)\end{array}$

Equation (8) can now be rewritten in terms of y_max, t_max, and α:

$\begin{array}{ll}y(t/t_{max})=y_{max}{e}^{\alpha }{\left(\frac{t}{t_{max}}\right)}^{\alpha }{e}^{-\frac{t}{t_{max}}\alpha }& (12)\end{array}$

Taking the natural logarithm on both sides of Equation (12) produces

$\begin{array}{ll}log(y(t^{\prime }))=log(y_{max})+\alpha (1+log(t^{\prime })-t^{\prime })& (13)\end{array}$

where $t^{\prime }=t/t_{max}$. This equation has the form y = C + αx, when t₀ and t_max are known, y_max and α can be determined from the linear regression of the natural logarithm of the observed values, y(t′), with (1 + log(t′) − t′) as the independent variable.

As shown above, the LL-LS method requires additional information about t₀ and t_max to be obtained from the observed data. In practice, t₀ is usually estimated from the time prior to t_max at which S(t) is within one standard deviation (SD) of the initial baseline signal (Chan and Nelson, 2004), and the location of t_max is found by calculating the centroid of the curve (Madsen, 1992). Although the LL-LS method provides a simplified way of assessing rCBV while avoiding the initial-value problem, and it has already been applied in several clinical applications (Chan and Nelson, 2004; Wu et al., 2004), this method still performs worse than the nonlinear method. This can be due to the logarithm operation complicating the statistics of the recorded data, resulting in the noise no longer conforming to a Gaussian distribution. When sufficient data points are available, the noise can still be treated as an approximate Gaussian distribution (according to the Lévy-Lindeberg Central Limit Theory; see Appendix B) (Casella and Berger, 2001). However, when only a few data points are available or the data are contaminated by outliers, the LS estimates become biased and are no longer characterized by minimum variance. Moreover, the estimation of t_max also introduces extra deviation.

5. Proposed Statistical Optimization Method

Assuming t₀ = 0, the first pass of the concentration–time course ΔR2^*(t) is rewritten as follows:

$\begin{array}{ll}y(t)=A_{rCBV}\frac{1}{{\beta }^{\alpha }\Gamma (\alpha )}{t}^{\alpha -1}{e}^{-t/\beta }& (14)\end{array}$

where A_rCBV is the value of rCBV, $\Gamma (\alpha )={\int }_0^{\infty }{x}^{\alpha -1}{e}^{-x}dx$ is the gamma function, t > 0, α > 0, and β > 0. Note that

$\begin{array}{ll}f(t|\alpha ,\beta )=\frac{1}{{\beta }^{\alpha }\Gamma (\alpha )}{t}^{\alpha -1}{e}^{-t/\beta }& (15)\end{array}$

is the probability density function (p.d.f.) of the gamma distribution, and

$\begin{array}{ll}{\displaystyle {\int }_0^{\infty }}f(t|\alpha ,\beta )dt=1.& (16)\end{array}$

Let Y denotes the observed ΔR2^*(t) data of the first bolus pass Y = (y(t₁), y(t₂), …, y(t_k), …), k = 1, …, N that specifics a gamma distribution f(t | α, β) with unknown parameters α, and β. Suppose there is a discrete random sample X = (x₁, x₂, …, x_n) of n independent and identically distributed (i.i.d.) observations, which obey the given gamma distribution f(t | α, β). The random variable x may only take N different values t₁, t₂, …, t_N, and y(t_k) gives the number of variables obtained in time t_k. In other words, the concentration-time curve is considered to be the frequency histogram of the random sample X with a gamma distribution f(t | α, β), where the occurrence frequency of t_k in observation set X is proportional to y(t_k) (Figure 2). This transforms, the original perfusion curve-fitting problem into a gamma-distribution-fitting problem. MLE is a popular choice for estimating parameters of a gamma distribution due to its optimal asymptotic properties (Gupta and Groll, 1961; Harter and Moor, 1965), we therefore applied an MLE approach to solve this problem (Myung, 2003).

FIGURE 2

Figure 2. Schematic of the proposed approach. The concentration curve ΔR2^*(t) is considered the frequency histogram of a random discrete sample X = (x₁, x₂, …, x_n) with a gamma distribution f(t | α, β) associated with the ΔR2^*(t) data. The original perfusion curve fitting problem then can be transformed to a gamma distribution fitting problem.

Denote Y = (y(t₁), y(t₂), …, y(t_k), …). The likelihood function is defined as:

$\begin{array}{ll}L(\alpha ,\beta )|Y)={\displaystyle \prod _{k=1}^n}f(y(t_i)|\alpha ,\beta )& (17)\end{array}$

Maximizing L(α, β|Y) with respect to (α, β) is equivalent to minimizing the negative logarithm of L(α, β|Y), and hence,

$\begin{array}{l}logL(\alpha ,\beta )=-n\alpha \text{ log}\beta -nlog\Gamma (\alpha )\\ +(\alpha -1){\displaystyle \sum }logx-\frac{1}{\beta }{\displaystyle \sum _{k=1}^n}{x}_k& (18)\end{array}$

where the summation is over the n sample values.

Assuming that the log-likelihood function is differentiable, it must satisfy the following partial differential equation

$\begin{array}{ll}\{\begin{array}{c}\frac{\partial (logL(\alpha ,\beta ))}{\partial \alpha }=0\\ \frac{\partial (logL(\alpha ,\beta ))}{\partial \beta }=0\end{array}& (19)\end{array}$

Differentiating Equation (18) we obtain the equations for the MLE:

$\begin{array}{ll}log\widehat{\beta }+\frac{\partial }{\partial \widehat{\alpha }}log\Gamma (\widehat{\alpha })-\frac{1}{n}{\displaystyle \sum _{k=1}^n}log{x}_k=0.& (20)\end{array}$

$\begin{array}{ll}\widehat{\alpha }\widehat{\beta }=\frac{1}{n}{\displaystyle \sum _{k=1}^n}{x}_k& (21)\end{array}$

Substituting $\widehat{\beta }$ in Equation (20) gives

$\begin{array}{ll}log\alpha -\psi (\alpha )=log\overline{x}-\frac{1}{n}{\displaystyle \sum _{k=1}^n}log{x}_k& (22)\end{array}$

where the digamma function $\psi (\widehat{\alpha })=\frac{\partial }{\partial \widehat{\alpha }}log\Gamma (\widehat{\alpha })$ is the logarithmic derivative of Γ(α), and $\overline{x}=1/n\sum _{k=1}^{n}x(t_k)$ denotes the average of the given data. Equation (22) is implicit in $\widehat{\alpha }$, and it is neither possible to find an analytical expression for $(\widehat{\alpha },\widehat{\beta })$ nor to use the Davis tables of ψ-functions directly (Davis, 1933). Masuyama and Kuroiwa prepared tables of $log\widehat{\alpha }-\psi (\widehat{\alpha })$ for the likelihood solutions of the gamma distribution (Masuyama and Kuroiwa, 1951; Greenwood, 1960).

Furthermore, the second derivatives of the log-likelihoods

$\begin{array}{ll}\{\begin{array}{c}\frac{{\partial }^{2}logL(\alpha ,\beta )}{\partial {\alpha }^2}={\psi }_1(\widehat{\alpha })<0\\ \frac{{\partial }^{2}logL(\alpha ,\beta )}{\partial {\beta }^2}=-\frac{\overline{x}}{{\widehat{\beta }}^2}<0\end{array}& (23)\end{array}$

where ${\psi }_1(\widehat{\alpha })=\frac{\partial }{\partial \widehat{\alpha }}\psi (\widehat{\alpha })$, are 1-digamma functions. The second derivatives are negative, and ensure that the log-likelihood function $logL(\widehat{\alpha },\widehat{\beta })$ is a maximum.

In order to obtain a numerical expression, Thom developed an approximation to the solution (Thom, 1958). The digamma function, $\psi (\widehat{\alpha })$ has an asymptotic expansion

$\begin{array}{ll}\psi (\alpha )=log\alpha -1/(2\alpha )-{\displaystyle \sum _{k=1}^m}{(-1)}^{k-1}{B}_k/(2k{\alpha }^{2k})+R_m,& (24)\end{array}$

where B_k is the Bernoulli number for B₁ = 1/6, B₂ = 1/30, …, and R_m is the remainder after m terms.

When α ≥ 1,

$\begin{array}{ll}\left|R_m\right|\le \frac{B_{m+1}}{(2m+2){\alpha }^{2m+2}}& (25)\end{array}$

and can be neglected. The approximation is more accurate when α is larger.

For m = 1 we have

$\begin{array}{ll}\psi (\alpha )=log\alpha -1/(2\alpha )-1/(12{\alpha }^2).& (26)\end{array}$

Substituting in Equation (22) yields

$\begin{array}{ll}12\left(log\overline{x}-\frac{1}{n}{\displaystyle \sum _{k=1}^n}log{x}_k\right){\widehat{\alpha }}^2-6\widehat{\alpha }-1=0& (27)\end{array}$

Simplifying by letting $A=log\overline{x}-\frac{1}{n}\sum _{k=1}^{n}log{x}_k$ produces

$\begin{array}{ll}12A{\widehat{\alpha }}^2-6\widehat{\alpha }-1=0,& (28)\end{array}$

which is a quadratic equation whose only pertinent root is

$\begin{array}{ll}\widehat{\alpha }=\frac{1+\sqrt{1+4A/3}}{4A}.& (29)\end{array}$

Combining this with Equation (21) gives the MLE for the gamma distribution. Thom provided an error table for correcting the estimate obtained from Equation (29) in the case of a small α (Thom, 1958), seen in Table 1.

TABLE 1

Table 1. Errors for correcting the α estimate.

Note that y(t_k) is a scaled version of the occurrence frequency of t_k in observation set X. We have

$\begin{array}{ll}\overline{x}=\frac{{\displaystyle {\sum }_{k=1}^N}{y}_k{t}_k}{{\displaystyle {\sum }_{k=1}^N}{y}_k}.& (30)\end{array}$

To summarize, A, α, and β are expressed in terms of the observation data y_k and associated time t_k as follows:

$\begin{array}{ll}\{\begin{array}{l}A=log\frac{{\displaystyle {\sum }_{k=1}^N}{y}_k{t}_k}{{\displaystyle {\sum }_{k=1}^N}{y}_k}-\frac{{\displaystyle {\sum }_{k=1}^N}{y}_{k}log{t}_k}{{\displaystyle {\sum }_{k=1}^N}{y}_k}\\ \widehat{\alpha }\widehat{\beta }=\frac{{\displaystyle {\sum }_{k=1}^N}{y}_k{t}_k}{{\displaystyle {\sum }_{k=1}^N}{y}_k}\end{array}& (31)\end{array}$

6. Experimental Validation

In this section, we examine the effectiveness of the proposed approach and compare the performances of the three fitting algorithms based on the concentration-time curve. Because the real values corresponding to the DSC MRI data are unknown, numerical simulations were used to investigate the uncertainties of the estimated cerebrovascular parameters (i.e., rCBV or F) for the three fitting algorithms with different signal-to-noise ratios (SNR) and time resolution (Δt) values.

Two subjects participated in this study. Images were acquired with a 1.5-tesla scanner (Signa Excite, GE). The CBV imaging sequence consisted of 30 T2*-weighted images that were collected using an echo-planar gradient-echo (GE) sequence with the following parameters: repetition time = 3, 100 ms, echo time = 80 ms, flip angle = 90°, field of view = 240 × 240mm, matrix size = 128 × 128, and 0.1 mmol/kg Gd-DTPA administered with a power injector. It should be noted that, in order to achieve a sufficient SNR and coverage of the whole brain, the repetition time was increased to be 3.1 s, since the typical value is about 1 s (Zhang et al., 2016).

Since our method is used to estimate the whole-brain rCBV, involving tens of thousands of data points, we presented only a few examples for a more detailed comparison. Thirty voxels were chosen as the region of interest (ROI) from the two subjects. They were picked from various areas, such as the motor area, visual area, and the thalamus. The concentration-time curves were created for each voxel using (Equation 1). Fitting procedures were applied to these curves using three fitting algorithms (nonlinear L-M method, LL-LS method, and the proposed method), which produced a total of 90 parameter sets {α, β}. Figure 3 gives examples of perfusion-curve fitting for one voxel when using the three methods. The figure demonstrates that all of the parameter sets had a reasonable physiological meaning, and so the three methods could be compared fairly. The ideal concentration-time curves as gamma-variate functions were then generated from these sets with the model

$\begin{array}{ll}\Delta R{2}^{*}(t)=\Delta R{2}^{*}(t)+K{\displaystyle {\int }_0^t}\Delta R{2}^{*}(t^{\prime })dt& (32)\end{array}$

where K = 0.02 is the recirculation weighting factor that describes the recirculation and leakage during the recirculation phase of the curve (Weisskoff et al., 1994). Gaussian-distributed noise was added to ΔR2^*(t) at different levels (SNR = 5, 10, 20, 50, and 100) with SNR defined as Chan and Nelson (2004)

$\begin{array}{ll}\text{SNR}=\frac{y_{max}}{\sigma }& (33)\end{array}$

where y_max is the peak of the curves (see Equation 11), and σ is the standard deviation of the added noise. The time resolution Δt was chosen to vary between 0.2 and 3.2 s in steps of 0.2 s. Then the gamma-variate fitting was repeated using the three methods. The onset of the bolus at t₀ was determined by searching from the time point of the maximum back to zero looking for two successive values that were less than a threshold of 10% of the maximum (Benner et al., 1997). The found time points were additionally corrected according to the time resolution and were identical for the three fitting algorithms. Moreover, the estimated parameter {α, β} for the LL-LS method, except where indicated otherwise, was also used as the initial value in the nonlinear fitting procedure.

FIGURE 3

Figure 3. Example of fitting curve obtained using the proposed method (red line, α_mle = 5.8943, β_mle = 1.6095, F_mle = 2.4427), the LL-LS method (green line, α_ls = 2.4364, β_ls = 5.6789, F_ls = 3.0727), and the nonlinear method (blue line, α_non = 6.0289, β_non = 1.7118, F_non = 2.5780) for a single voxel. Filled circles denote valid data points for perfusion-curve fitting; other points are indicate by hollow circles.

Two uncertainty values were calculated for the rCBV (i.e., the area under the concentration-time curves) to evaluate the three methods, Benner et al. (1997)

$\begin{array}{ll}{\mu }_i=\frac{\overline{F}-F_i}{F_i}·100\%& (34)\end{array}$

and

$\begin{array}{ll}{\sigma }_i=\frac{{\sigma }_F}{F_i}·100\%& (35)\end{array}$

where $\overline{F}$ and σ_F are the mean and standard-deviation values of rCBV calculated from 250 synthetic curves with parameter set {α_i, β_i}, and F_i is the real rCBV value calculated from Equation (4). Taken overall 90 parameter sets {α, β}, μ gives the percentage of deviation from the correct value and is a measure of how accurately the parameters calculated from the fit agree with the correct ones, and σ is the coefficient of variation as a percentage of all values calculated for F and represents a measure of the probability that a calculated value lies within a certain range of the correct value (Benner et al., 1997).

Figure 4 shows the calculated uncertainty values for the three fitting methods as functions of the time resolution for the five SNR levels. In general, the quality of the fit depends on the time resolution and the noise level, with the noise having a greater impact on the fitting result when Δt is larger. For all three methods, the uncertainties (μ and σ) of the estimated parameter increased with increasing Δt and with decreasing SNR. The influence of SNR variation was stronger for σ (left panels in Figure 4) than for μ (right panels in Figure 4). The estimated parameter was more stable while Δt decreased and SNR increased (right panels in Figure 4). These observations can be explained by a smaller Δt and larger SNR resulting in more usable sample points with less contamination by noise, thereby achieving more reliable fitting. However, even for the smallest Δt = 0.2 and the lowest noise level (SNR = 100), the estimated parameter deviated from the real values by more than 10% (left panels in Figure 4). This inherent deviation can contribute to the additional effect of contrast agent recirculation in the simulation. In all situations, the nonlinear method (Figures 4E,F) performed better than the LL-LS method (Figures 4C,D), while the proposed approach (Figures 4A,B) exhibited far superior stability and accuracy than the other two methods. While both the conventional LL-LS and nonlinear methods delivered decent estimates when the SNR was large (Figure 4), they led to a failure of the calculated parameters for the lowest SNR (SNR = 5) and the median Δt case (Δt > 1.6 s) where the uncertainties of the estimated parameter were larger than 50% (blue lines in Figures 4C–F). In contrast, the proposed approach generated significantly better results for different time resolutions and noise levels, illustrating the benefits of transforming the original curve-fitting problem into a distribution-fitting problem. Moreover, the fitting procedure for the LL-LS and nonlinear methods failed for the largest time resolution (Δt = 3.2 s) and the highest noise level (SNR = 5). Figure 5 shows the failure rate of the fitting procedure as a function of the time resolution (Δt). Not surprisingly, the number of failures increased with increasing Δt, demonstrating that obtaining more sample points during the first pass of the contrast agent can help to build a more accurate view of the vascular indicator dynamics (Lu and Monahan, 1993).

FIGURE 4

Figure 4. Uncertainty μ (left column) and σ (right column) values calculated using three fitting methods as functions of the time resolution for five SNR values. Note that the ordinate scale is smaller for the proposed method than for other two methods for clarity. The fitting procedure was considered to have failed if the values of the μ and σ uncertainties were larger than 50% of the correct values. From top to bottom: proposed method, LL-LS method, and nonlinear method. (A) Uncertainty μ values calculated using the proposed method. (B) Uncertainty σ values calculated using the proposed method. (C) Uncertainty μ values calculated using LL-LS method. (D) Uncertainty σ values calculated using LL-LS method. (E) Uncertainty μ values calculated using nonlinear method. (F) Uncertainty σ values calculated using nonlinear method.

FIGURE 5

Figure 5. Failures of the fitting procedure with the LL-LS and nonlinear methods as functions of the time resolution (Δt) for SNR = 5. The fitting procedure was considered to have failed if the algorithms did not converge or the calculated values of parameters α and β were imaginary numbers. The estimated parameter {α, β} obtained using the LL-LS method is used as the initial value in the nonlinear fitting procedure. The number of failures of the nonlinear method is thus equal to the number of failures of the LL-LS method plus the number of times that the algorithm for the nonlinear diverged. The number of failures consistently increased as the time resolution decreased. However, there were no failures for the proposed method even for the highest noise level.

7. Discussion and Conclusion

Indicator dilution analysis involves computing perfusion-related parameters (e.g., blood volume, blood flow, and mean transit time) from the observed flow of a contrast agent passing through the vascular system. This analysis method is applied in a diverse range of medical fields. The measurements might yield valuable diagnostic information for use in various applications, such as the assessment of tissue perfusion, the detection of ischemic regions and cancer hypervascularization (Ruediger et al., 1964; Eckersley et al., 2002; Yang et al., 2003). The available measurement methods extend beyond DSC-MRI, to include other bolus-tracking imaging methods, such as contrast-enhanced computed tomography (Pienn et al., 2013), contrast-enhanced ultrasound imaging (Feinstein, 2004), and scintigraphy (Anger, 1964). In this paper, we have proposed a novel statistical optimization strategy for perfusion analysis, with the initial motivation being that the voxel-wise evaluation of the rCBV information needs to be performed from a series of DSC-MRI observations in a real clinical application (Zhang et al., 2016). By designing the concentration-time curves as a random sample from a gamma distribution with time as the random variable, we have transformed the original perfusion curve-fitting problem into a gamma-distribution-fitting problem; MLE is then adopted to solve this problem. Since real values from a real experiment were not available, we used simulation experiments to compare the proposed method with conventional methods. The obtained simulation results have demonstrated that the proposed strategy has the potential to provide a more stable and accurate perfusion analysis than a conventional curve-fitting strategy for different time resolutions and noise levels.

We consider that the improved performance of the proposed method is attributable to the reasons given below. The LL-LS method implicitly assumes that logy(t) conforms to a Gaussian distribution with the right-hand side (RHS) of Equation (13) as its mean and a constant variance. Hence, the small fitting error near t = t₀ (we assume t₀ = 0), that is, at y(t₀) = 0, could cause a large error at logy(t₀). Then, the large change in the left-hand side (LHS) of Equation (13) results a larger error in the estimation of α than for the LS solution based on Equation (13). Unlike the LL-LS method, the proposed method considers that all the observations y(t_k), k = 1, …N are independent samples from a gamma distribution (Equation 15) with unknown parameters α and β. The MLE method estimates those two parameters by maximizing the probability of obtaining the observed data. Moreover, compared with the other two methods, the LL-LS method requires additional information on t_max that also increases the deviation.

Figure 6 provides a comparison of the three methods for different data sets. When the data point near t = t₀ was excluded from the fitting, the LL-LS method obtained estimates similar to those from the other two methods, illustrating the great computational weight for this data point (compare green lines in Figures 3, 6A). In contrast, including the data point far from t = t₀ had little influence on the results (compare green lines in Figures 3, 6B). Because both the nonlinear and proposed methods minimize the estimation error averaged over all observations, they still show relatively stable estimates irrespective of the choice of data points (Figure 6B). This is particularly important for the voxel-wise analysis in the whole brain where it is impossible to check all the data points. The width of the first bolus pass shows a wide variation due to different imaging regions and different molarity/dosages for the contrast agents (Wirestam et al., 2009; Schmiedeskamp et al., 2013). This situation may result in valuable data points being missed or extra data points being included. Our algorithm exhibits more stable and accurate performance in all simulation situations, and hence it has the potential to cope better with undesirable aspects of the image acquisition, such as the inclusion of improper data points, unforeseen data errors, poor image alignment, and the use of new protocols/tracers.

FIGURE 6

Figure 6. (A) Compared with Figure 3, when the first point near to t = t₀ was excluded from the fitting, the three methods produced very similar results, illustrating that this point have greater computational weight than other point in the LL-LS method (α_mle = 6.9087, β_mle = 1.3795, F_mle = 2.4313, α_ls = 5.5731, β_ls = 1.7838, F_ls = 2.5687, α_non = 6.0294, β_non = 1.7116, F_non = 2.5780). (B) In comparison with the first point near to t = t₀, including the additional last point hardly influenced the results. The nonlinear method and the proposed method showed relatively stable performances for the choice of data points due to all observations being treated in an equal manner (α_mle = 5.8765, β_mle = 1.7348, F_mle = 2.6190, α_ls = 4.4527, β_ls = 2.3626, F_ls = 2.6610, α_non = 5.6077, β_non = 1.8767, F_non = 2.6624). The ΔR2^* signal obtained during the first bolus pass is shown as filled circles, while the other sample points are indicated by hollow circles. The area under the fitted ΔR2^* curve corresponds to the estimated rCBV.

Our study also provides suggestions on how to achieve more reliable estimates when using the two conventional methods in practical applications. As discussed above, the data point near t = t₀ (we suggest t < 1 s) should be excluded from the fitting when using the LL-LS method. Moreover, compared with the previous empirical approach for the initial value in the nonlinear method, using estimates from the LL-LS method as the initial values in the nonlinear method can reduce the probability of failure for the fitting procedure by about 10-fold (Benner et al., 1997).

Different from a previous method (Scalzo and Liebeskind, 2016), the parameters were computed using a bolus tracking method based on the deconvolution of the time-density curve on a pixel-by-pixel basis. At first, our method was proposed to be applied to the specific situation of estimating the whole brain rCBV curve (Zhang et al., 2016), which means fitting the perfusion curve by involving only a few data points. As we mentioned above, in the case of a small sample, our proposed strategy can also ensure a more accurate and stable performance. On the other hand, for the parameter estimation, although his paper uses the EM algorithm, which is based on the MLE principle. But in their M-step, they adopt a numerical optimization technique. This iterative strategy can approach the optimal solution well, but it may increase the time cost. Then they used the Bayesian Information Criterion (BIC) to determine the optimal K. As we have known that BIC requires a certain amount of data to ensure accuracy, and insufficient data points will lead to inaccurate estimation. Compared with this numerical optimization, we adopt an approximate asymptotic expansion to determine the optimal parameter. Specifically, according to an approximation that Thom developed to the solution, the digamma function, $\psi (\widehat{\alpha })$ has an asymptotic expansion (Thom, 1958). Therefore, this approximate solution our proposal adopts will allow us to quickly determine the parameters, greatly saving the time cost, and there is no need to manually set the solution procedures.

Furthermore, the proposed strategy could potentially be improved by constructing several bias-correction forms of the MLE estimator to improve the optimal properties of the original MLE in the case of small samples (Choi and Wette, 1969; Giles and Feng, 2009). Also, the strategy could be extended to an estimation scheme containing location parameter t₀ for the three-parameter gamma distribution (Cohen and Whitten, 1982; Anger, 1995). In addition, it is possible to accommodate other mathematical models and include prior knowledge from other modalities so as to further improve the performance (Wu et al., 2004). While it is always possible to transform a distribution-fitting problem into a curve-fitting problem by fitting a curve to a histogram, it is rarely possible to perform the reverse procedure. However, in this paper, we have presented such an example. Despite its success, several questions related to the mathematics underlying the approach still need to be clarified, which will limit our proposal to some extent in explaining the underlying mechanism, such as the optimal asymptotic properties of the MLE estimator for discrete sampling, and the relationship between time resolution and sample size.

In conclusion, the results reported here have demonstrated that the proposed strategy exhibits attractive performance relative to the conventional methods. The algorithm is easy to implement, and it is equivalent to the LL-LS method in terms of computational cost O(n), and it dispenses with the need to guess the initial values. Our method involves transforming the curve-fitting problem into the probability distribution fitting problem, and the perfusion curve can be well fitted by the MLE estimator. Thus, it is not only suitable for contrast-enhanced MRI but also other methods of bolus-tracking perfusion imaging, for example, X-rays, CT, PET, etc. We therefore suggest that the proposed strategy could represent an alternative option for assessing intravascular indicator dynamics from bolus-tracking perfusion imaging. Currently, the method is successfully applied to our study on hemodynamic assimilation (Zhang et al., 2016).

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics Statement

The studies involving human participants were reviewed and approved by College of Science, Zhejiang University of Technology. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

ZH conceived and designed the experiments, performed the experiments and analyzed the data, and wrote the manuscript. ZH, FL, and QL edited and approved the final version of the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work is supported in part by National Key Research and Development Program of China under Grant 2018YFA0701400, in part by the Public Projects of Science Technology Department of Zhejiang Province under Grant LGF20H180015.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Anger, H. O. (1964). Scintillation camera with multichannel collimators. J. Nuclear Med. 5, 515–531.

Google Scholar

Anger, H. O. (1995). Maximum likelihood parameters estimation in the three-parameter gamma distribution. Comput. Stat. Data Anal. 20, 343–354. doi: 10.1016/0167-9473(94)00050-S

CrossRef Full Text | Google Scholar

Arzanforoosh, F., Croal, P. L., van Garderen, K. A., Smits, M., Chappell, M. A., and Warnert, E. A. (2021). Effect of applying leakage correction on rcbv measurement derived from dsc-mri in enhancing and nonenhancing glioma. Front. Oncol. 11:648528. doi: 10.3389/fonc.2021.648528

PubMed Abstract | CrossRef Full Text | Google Scholar

Axel, L. (1980). Cerebral blood flow determination by rapid-sequence computed tomography. Radilogy 137, 679–686. doi: 10.1148/radiology.137.3.7003648

PubMed Abstract | CrossRef Full Text | Google Scholar

Benner, T., Heiland, S., Erb, G., Forsting, M., and Sartor, K. (1997). Accuracy of gamma-variate fits to concentration-time curves from dynamic susceptibility-contrast enhanced MRI: influence of time resolution, maximal signal drop and signal-to-noise. Magn. Reson. Imaging 15, 307–317. doi: 10.1016/S0730-725X(96)00392-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Casella, G., and Berger, R. L. (2001). Statistical Inference, 2nd Edn. Pacific Grove, CA: Duxbury Resource, Center.

Google Scholar

Chan, A. A., and Nelson, S. J. (2004). “Simplified gamma-variate fitting of perfusion curves,” in 2th IEEE International Symposium on Biomedical Imaging (ISBI) (Arlington, VA), 1067–1070.

Google Scholar

Choi, S. C., and Wette, R. (1969). Maximum likelihood estimation of the parameters of the gamma distribution and their bias. Technometrics 11, 683–690. doi: 10.1080/00401706.1969.10490731

PubMed Abstract | CrossRef Full Text | Google Scholar

Cohen, A. C., and Whitten, B. J. (1982). Modified moment and maximum likelihood estimators for parameters of the three-parameter gamma distribution. Commun. Stat. Simul. Comput. 11, 197–216. doi: 10.1080/03610918208812254

CrossRef Full Text

Davenport, R. (1983). The derivation of the gamma-variate relationship for tracer dilution curves. J. Nuclear Med. 24, 945–948.

PubMed Abstract | Google Scholar

Davis, H. T. (1933). Tables of Higher Mathematical Functions, Vol. I, II. Bloomington, IN: Principia Press.

Eckersley, R. J., Sedelaar, J. P., Blomley, M. J., Wijkstra, H., de Souza, N. M., Cosgrove, D. O., et al. (2002). Quantitative microbubble enhanced transrectal ultrasound as a tool formonitoring hormonal treatment of prostate carcinoma. Prostate Cancer Prostatic Dis. 51, 256–267. doi: 10.1002/pros.10065

PubMed Abstract | CrossRef Full Text | Google Scholar

Emblem, K. E., Farrar, C. T., Gerstner, E. R., Batchelor, T. T., Borra, R. J. H., Rosen, B. R., et al. (2014). Vessel calibre a potential MRI biomarker of tumour response in clinical trials. Nat. Rev. Clin. Oncol. 11, 566–584. doi: 10.1038/nrclinonc.2014.126

PubMed Abstract | CrossRef Full Text | Google Scholar

Feinstein, S. B. (2004). The powerful microbubble: from bench to bedside, from intravascular indicator to therapeutic delivery system, and beyond. Am. J. Physiol. Heart Circ. Physiol. 287, 450–457. doi: 10.1152/ajpheart.00134.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Gavin, H. P. (2019). The Levenberg-Marquardt Algorithm for Nonlinear Least Squares Curve-Fitting Problems. Department of Civil and Environmental Engineering, Duke University, 1–19.

Google Scholar

Giles, D. E., and Feng, H. (2009). Bias of the Maximum Likelihood Estimators of the Two-Parameter Gamma Distribution Revisited. Technical report, Department of Economics, University of Victoria.

Google Scholar

Greenwood, J. A. (1960). Aids for fitting the gamma distribution by maximum likelihood. Technometrics 2, 55–65. doi: 10.1080/00401706.1960.10489880

CrossRef Full Text | Google Scholar

Gupta, S. S., and Groll, P. A. (1961). Gamma distribution in acceptance sampling based on life tests. J. Am. Stat. Assoc. 56, 942–970. doi: 10.1080/01621459.1961.10482137

CrossRef Full Text | Google Scholar

Harter, H. L., and Moor, A. H. (1965). Maximum likelihood estimation of the parameters of gamma and weibull populations from censored samples. Technometrics 7, 639–643. doi: 10.1080/00401706.1965.10490304

CrossRef Full Text | Google Scholar

Hu, Z. H., Ni, P. Y., Liu, C., Zhao, Z. H., Liu, H. F., and Shi, P. C. (2012). Quantitative evaluation of activation state in functional brain imaging. Brain Topogr. 25, 362–373. doi: 10.1007/s10548-012-0230-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, S., and Cho, H. J. (2021). Model-free leakage index estimation of the blood-brain barrier using dual dynamic susceptibility contrast mri acquisition. NMR Biomed. 13:e4570. doi: 10.1002/nbm.4570

PubMed Abstract | CrossRef Full Text | Google Scholar

Kosior, J. C., and Frayne, R. (2010). Perfusion parameters derived from bolus-tracking perfusion imaging are immune to tracer recirculation. J. Magn. Reson. Imaging 31, 753–756. doi: 10.1002/jmri.22052

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, K. L., Zhu, X. P., Checkley, D. R., Tessier, J. J. L., Hillier, V. F., Waterton, J. C., et al. (2003). Simultaneous mapping of blood volume and endothelial permeability surface area product in gliomas using iterative analysis of first-pass dynamic contrast enhanced MRI data. Br. J. Radiol. 76, 39–50. doi: 10.1259/bjr/31662734

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X. F. (2014). MRI Perfusion-Weighted Imaging Analysis. Dordrecht: Springer Netherlands.

Google Scholar

Lu, D., and Monahan, W. G. (1993). Effect of sample numbers on the kinetic data analysis of MR contrast agents. Magn. Reson. Med. 30, 131–134. doi: 10.1002/mrm.1910300120

PubMed Abstract | CrossRef Full Text | Google Scholar

Madsen, M. T. (1992). A simplified formulation of the gamma variate function. Phys. Med. Biol. 37, 1597–1600. doi: 10.1088/0031-9155/37/7/010

CrossRef Full Text | Google Scholar

Masuyama, M., and Kuroiwa, Y. (1951). Tabel for the likelihood solution of gamma distributions and its medical applications. Rep. Stat. Appl. Unition Jpn. Sci. Eng. 7, 18–23.

Mischi, M., den Boer, J. A., and Korsten, H. H. M. (2008). On the physical and stochastic representation of an indicator dilution curve as a gamma variate. Physiol. Meas. 29, 281–294. doi: 10.1088/0967-3334/29/3/001

PubMed Abstract | CrossRef Full Text | Google Scholar

Myung, J. (2003). Tutorial on maximum likelihood estimation. J. Math. Psychol. 47, 90–100. doi: 10.1016/S0022-2496(02)00028-7

CrossRef Full Text | Google Scholar

Norman, D., Axel, L., Berninger, W. H., Edwards, M. S., Cann, C. E., Redington, R. W., et al. (1981). Dynamic computed tomography of the brain: techniques, data analysis, and applications. Am. J. Roentgenol. 136, 759–770. doi: 10.2214/ajr.136.4.759

PubMed Abstract | CrossRef Full Text | Google Scholar

Perkiö, J., Aronen, H. J., Kangasmäki, A., Liu, Y., Karonen, J., Savolainen, S., et al. (2002). Evaluation of four postprocessing methods for determination of cerebral blood volume and mean transit time by dynamic susceptibility contrast imaging. Magn. Reson. Med. 47, 973–981. doi: 10.1002/mrm.10126

PubMed Abstract | CrossRef Full Text | Google Scholar

Pianykh, O. S. (2012). Perfusion linearity and its applications in perfusion algorithm analysis. Comput. Med. Imaging Graph. 36, 204–214. doi: 10.1016/j.compmedimag.2011.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Pienn, M., Kovacs, G., Tscherner, M., Johnson, T. R., Kullning, P., Stollberger, R., et al. (2013). Determination of cardiac output with dynamic contrast-enhanced computed tomography. Int. J. Cardiovasc. Imaging 29, 1871–1878. doi: 10.1007/s10554-013-0279-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Quarles, C. C., Bell, L. C., and Stokes, A. M. (2019). Imaging vascular and hemodynamic features of the brain using dynamic susceptibility contrast and dynamic contrast enhanced mri. Neuroimage 187, 32–55. doi: 10.1016/j.neuroimage.2018.04.069

PubMed Abstract | CrossRef Full Text | Google Scholar

Rempp, K. A., Brix, G., Wenz, F., Becker, C. R., and Lorenz, F. G. W. J. (1994). Quantification of regional cerebral blood flow and volume with dynamic susceptibility contrast-enhanced MR imaging. Radiology 193, 637–641. doi: 10.1148/radiology.193.3.7972800

PubMed Abstract | CrossRef Full Text | Google Scholar

Romain, B., Rouet, L., Ohayon, D., Lucidarme, O., d'Alché Buc, F., and Letort, V. (2017). Parameter estimation of perfusion models in dynamic contrast-enhanced imaging: a unified framework for model comparison. Med. Image Anal. 35, 360–374. doi: 10.1016/j.media.2016.07.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosen, B. R., Belliveau, J. W., Buchbinder, B. R., McKinstry, R. C., Porkka, L. M., Kennedy, D. N., et al. (1991). Contrast agents and cerebral hemodynamics. Magn. Reson. Med. 19, 285–292. doi: 10.1002/mrm.1910190216

PubMed Abstract | CrossRef Full Text | Google Scholar

Ruediger, E. P., Knopp, M. V., Hoffmann, U., Milker, S. Z., and Brix, G. (1964). Multicompartment analysis of gadolinium chelate kinetics: blood-tissue exchange in mammary tumors by dynamic MR imaging. Circ. Res. 14, 502–515.

PubMed Abstract

Scalzo, F., and Liebeskind, D. S. (2016). Perfusion angiography in acute ischemic stroke. Comput. Math. Methods Med. 2016:2478324. doi: 10.1155/2016/2478324

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmiedeskamp, H., Andre, J. B., Straka, M., Christen, T., Nagpal, S., Recht, L., et al. (2013). Simultaneous perfusion and permeability measurements using combined spin- and gradient-echo MRI. J. Cereb. Blood Flow Metab. 33, 732–743. doi: 10.1038/jcbfm.2013.10

PubMed Abstract | CrossRef Full Text | Google Scholar

Thom, H. C. S. (1958). A note on the gamma distribution. Mon. Weather Rev. 86, 117–122. doi: 10.1175/1520-0493(1958)086<0117:ANOTGD>2.0.CO;2

CrossRef Full Text | Google Scholar

Thompson, H. K., Starmer, C. F., Whalen, R. E., and McIntosh, H. D. (1964). Indicator transit time considered as a gamma variate. Circ. Res. 14, 502–515. doi: 10.1161/01.RES.14.6.502

PubMed Abstract | CrossRef Full Text | Google Scholar

Weisskoff, R. M., Boxerman, J. L., Sorenson, A. G., Kulke, S. M., Campbell, T. A., and Rosen, B. R. (1994). “Simultaneous blood volume and permeability mapping using a single Gd-based contrast injection,” in Proceedings of International Society for Magnetic Resonance in Medicine (ISMRM) (San Francisco, CA), 279.

Wirestam, R., Engvall, C., Ryding, E., Holtås, S., Stånhlberg, F., and Reinstrup, P. (2009). Changes in cerebral perfusion detected by dynamic susceptibility contrast magnetic resonance imaging: normal volunteers examined during normal breathing and hyperventilation. J. Biomed. Sci. Eng. 2, 210–215. doi: 10.4236/jbise.2009.24034

CrossRef Full Text | Google Scholar

Wu, Y., An, H. Y., Krim, H., and Lin, W. L. (2004). An independent component analysis approach for minimizing effects of recirculation in dynamic susceptibility contrast magentic resonance imaging. J. Cereb. Blood Flow Metab. 27, 632–645. doi: 10.1038/sj.jcbfm.9600374

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, S., Law, M., Zagzag, D., Wu, H. H., Cha, S., Golfinos, J. G., et al. (2003). Dynamic contrast-enhanced perfusion MR imaging measurements of endothelial permeability: differentiation between atypical and typical meningiomas. Am. J. Neuroradiol. 24, 1554–1559.

PubMed Abstract | Google Scholar

Zhang, Y., Wang, Z., Cai, Z., Lin, Q., and Hu, Z. (2016). Nonlinear estimation of bold signals with the aid of cerebral blood volume imaging. Biomed. Eng. Online 15, 1–11. doi: 10.1186/s12938-016-0137-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Appendix

A. Calculating rCBV From Γ-Variate Function

The area F under the concentration time curve is proportional to the blood volume, and can be calculated as follows:

$\begin{array}{l}F={\displaystyle {\int }_{t_0}^{\infty }}{C}_i(t)dt\\ ={\displaystyle {\int }_{t_0}^{\infty }}A·{(t-t_0)}^{\alpha }{e}^{-(t-t_0)/\beta }dt\\ ={\displaystyle {\int }_0^{\infty }}A{t^{\prime }}^{\alpha }{e}^{-t^{\prime }/\beta }d{t}^{\prime }(t^{\prime }=t-t_0)\\ =A·{\beta }^{\alpha +1}{\displaystyle {\int }_0^{\infty }}(\frac{t^{\prime }}{\beta }{)}^{\alpha }{e}^{-t^{\prime }/\beta }d(\frac{t^{\prime }}{\beta })\\ =A·{\beta }^{\alpha +1}·\Gamma (\alpha +1)& (\mathrm{A1})\end{array}$

B. Lévy-Lindberg Central Limit Theory

Let X₁, X₂, …be a sequence of iid random variables with EX_i = μ and 0 < Var $X_i={\sigma }^2<\infty$. Define ${\overline{X}}_n=(1/n)\sum _{i=1}^n{X}_i$. Let G_n(x) denote the cumulative distribution function (cdf) of $\sqrt{n}({\overline{X}}_n-\mu )/\sigma$. Then, for any x, −∞ < x < ∞,

$\begin{array}{ll}{\displaystyle \underset{n\to \infty }{lim}}{G}_n(x)={\displaystyle {\int }_{-\infty }^x}\frac{1}{\sqrt{2\pi }}{e}^{-y^2/2}dy;& (\mathrm{B1})\end{array}$

that is, $\sqrt{n}({\overline{X}}_n-\mu )/\sigma$ has a limiting standard normal distribution (Casella and Berger, 2001).