|Home | About | Journals | Submit | Contact Us | Français|
The aim was to investigate the feasibility of making relative cerebral blood flow (rCBF) maps from MR images acquired with short TR by measuring the initial arrival amount of Gd-DTPA evaluated within a time window before any contrast agent has a chance to leave the tissue. We named this rCBF measurement technique utilizing the early data points of the Gd-DTPA bolus the “early time points” method (ET), based on the hypothesis that early time point signals were proportional to rCBF. Simulation data were used successfully to examine the ideal behavior of ET while monkey’s MRI results offered encouraging support to the utility of ET for rCBF calculation. A better brain coverage for ET could be obtained by applying the Simultaneous Echo Refocusing (SER) EPI technique. A recipe to run ET was presented, with attention paid to the noise problem around the time of arrival (TOA) of the contrast agent.
In susceptibility contrast-enhanced MRI of Gd-DTPA bolus injection, relative cerebral blood flow (rCBF) has commonly been analyzed using the deconvolution approach reported by Ostergaard et al. (Ostergaard et al., 1996b). Another strategy had been suggested by Buxton et al. (Buxton, 2002) to calculate rCBF by measuring the amount of Gd-DTPA arriving within a time window before any contrast agent has a chance to leave the tissue. That time window is in the order of the tissue mean transit time τ, a quantity known to be a few seconds for gray matter in humans. Given the basic assumption that the contrast agent has not had the chance to leave from the tissue, the quantity of Gd-DTPA present in the tissue will be proportional to local CBF. We name this rCBF measurement technique which utilizes the early contrast-enhanced MR data points the “early time points” method (ET). The basic assumption, which will be cited repeatedly, will be called ET’s basic assumption for the rest of the manuscript.
A perfusion analysis technique like the deconvolution method utilizes the arterial input function (AIF) from an artery to evaluate flow given by the expression below:
where C(t) is the contrast agent concentration curve, a is the perfusion term, AIF is arterial input function, is the convolution symbol and R(t) is the residue function which is the probability that a molecule of the agent that entered the capillary bed at t=0 is still there at time t (Buxton, 2002; Ostergaard et al., 1996b). Eq (1) can equally be used to explain the properties of ET. In ET, the residue function R(t) = 1 when ET’s basic assumption holds (Buxton, 2002), leading to Eq (2) below which shows that perfusion result evaluated by ET is independent of AIF.
The independence of ET results from AIF is associated with an important observation of the early time points. If the shape of AIF changes due to Gd-DTPA injection rate or some other factors, the time course C(t) of contrast concentration (ΔR2*) would change accordingly, but Eq (2) demonstrates that perfusion value estimated using the ET method does not change at all throughout the range of the early time points. A box-shaped or rectangle function rect(t) (Buxton, 2002) is a good model of a residue function that meets the demand of ET’s basic assumption and will be presented later in the Methods section.
One advantage of ET results being independent of AIF is skipping altogether the long known problem of possible inaccuracies in estimating the “true” AIF (Ellinger et al., 2000; Kiselev, 2001; Yamada et al., 2002), a problem that continues to be receiving a fair amount attention in the perfusion publications. Eq (2) is valid, however, only if AIF is the same for all pixels. The assumption of a single global AIF is applied throughout the manuscript for our investigation of the properties of ET. In the Discussion section, we will present a discussion on the problem of different AIF’s for different pixels due to delay and dispersion and why we chose to develop our understanding of ET using the global AIF assumption.
When it comes to transit delays of the bolus to the vascular bed, ET inherently requires an accurate correction of the transit delay to be made before the initiation of rCBF evaluation. In this manuscript, whenever the “early time points” or the “bolus time course” is mentioned, the transit delay is expected to have been corrected for the time course. The accuracy of transit delay correction depends, of course, on the accuracy of the estimation of the transit delay, a question that will be explored below.
While ET measures only rCBF, it can also derive absolute CBF the same way the deconvolution method does. It should be pointed out that the deconvolution outputs (the height of the maximum amplitude of the deconvolved residue function) using the method of Ostergaard et al. also yield only rCBF with respect to MR data analysis. One can get absolute CBF only if a proportionality constant, required for absolute CBF estimation, can be obtained independently and separately from a different technique like PET. One compares average PET CBF with the average MR blood flow derived from the deconvolution analysis to obtain the proportionality constant and then calculate absolute CBF. It is as easy for ET as for the deconvolution method to make absolute CBF maps by obtaining a similar proportionality constant using a calibration from PET. In this manuscript we focus only on the rCBF result of the ET analysis of MR data.
Since we are interested in the early time points of the rising bolus time course, a TR acquisition much shorter than the commonly used parameter of TR of ~1.5 sec in clinical perfusion imaging would be preferable. We proposed collecting MR data with echo planar imaging (EPI) acquisition using a short TR of 300ms. While TR=300ms is by no means an absolute requirement, it should cover the low range of τ, allowing a less ambiguous investigation of ET’s properties.
The tradeoff of short TR acquisition is the limited number of slices which reduces the clinical acceptance of ET. Therefore, in addition to applying standard gradient echo EPI (GE EPI) for our main MRI investigation, we briefly introduced at the end of this report a recently developed fast sequence, the Simultaneous Echo Refocusing (SER) EPI technique (Feinberg et al., 2002; Reese et al., 2002) which can acquire at least two slices in a single readout. SER could therefore be applied to obtain more slices than standard GE EPI in the same short TR. The utility of SER EPI had been reported in diffusion (Reese et al., 2009) and fMRI studies (Chan et al., 2008).
In order to learn how ET can be used for rCBF calculation, we need to identify the conditions and properties associated with the part of the bolus concentration time course which is compatible with ET’s basic assumption. The effort to clarify those conditions and properties of ET led to the investigation of four technical topics to be presented below. We supported our investigation using contrast-enhanced experimental MRI data and simulated perfusion data created using the algorithm of Wu et al. (Wu et al., 2003b).
The four technical topics (Topic I – IV below) of ET we chose to study highlight several of ET’s properties and limitations. Topics I clarifies ET’s limitations due to noise. Topic II deals with a number of different possible techniques used by ET to calculate rCBF. Topic III considers the question of which part of early time points would be useful for rCBF calculation. Topic IV was chosen to demonstrate whether rCBF results by ET would match perfusion results in the literature obtained previously by PET and MRI. All four topics reinforce each other in demonstrating the applicability of ET for rCBF evaluation. With regard to the blood flow terminology we employed in those four topics, we used flow ratio to describe the ratio between different flow levels of simulation data. With experimental contrast-enhanced MRI data, “flow ratio” and “rCBF” were used interchangeably.
The time of arrival (TOA) of the bolus is a critical ET parameter which has a different value for each pixel. The value of TOA is used for the correction of the transit delay of the bolus for each pixel. Once TOA is found and the transit delay corrected, any single early time point on the rising bolus time curve could theoretically be used for rCBF calculation provided that the same time location is used for every pixel. Unfortunately, the presence of MRI noise makes it difficult to identify the precise location of TOA. The small amount of contrast agent entering the pixel at TOA gives a very low contrast-to-noise ratio (CNR) of the MR signal of the bolus. A misestimated TOA would give an inaccurate evaluation of rCBF by ET.
We must notice there is also a different TOA/noise problem completely separate from the difficulty of identifying the true TOA. It is important to recognize that, even if the true TOA is known, as is the case in noiseless simulation data, low CNR could undermine significantly the accuracy of calculated rCBF results near TOA. In other words, the knowledge of true TOA is not sufficient to secure accurate rCBF results.
Therefore, there are two separate ways for noise around TOA to affect the accuracy of calculated rCBF. First, one may not be able to find the correct TOA due to low CNR at the beginning of the rising concentration time course. Secondly, even if the correct TOA value could somehow be found, low CNR at TOA would still affect the accuracy of calculated rCBF. It is important to analyze these two noise related questions on TOA separately in order not to mix up the interpretation of their different effects.
The usefulness of a TOA searching algorithm to estimate TOA can be assessed first with noiseless simulation data, then with noisy simulation data generated with a predetermined signal to noise ratio (SNR) close to what can be encountered in contrast-enhanced MRI settings. Given that there is currently no accepted optimal way to find TOA, the TOA searching algorithm we selected was just one of many possible options and we did not conduct a comparison of different competing algorithms. Using both noiseless and noisy simulation data (SNR= 20), we proposed to assess the reasonableness of the calculated TOA obtained with our algorithm.
What procedure can be used to calculate rCBF by ET? Buxton (Buxton, 2002) proposed to calculate rCBF by measuring the initial slope of the bolus signal coming in. However, there is no reason to restrict ourselves to only using the initial slope. We can use equally well the signal intensity (ΔR2*) of the early concentration curve to obtain the rCBF value. Eqs (1) and (2) above already showed that the signal intensity was good enough for rCBF evaluation. Eqs (3) and (4) below provided an explanation on why the slope value works equally well as signal intensity:
Eq (4) shows that we can use the slope or the signal intensity to obtain the same rCBF values. To our understanding, signal intensity has not been proposed previously to be used for ET to evaluate rCBF. In fact, as long as the assumption of the early contrast agent concentration being proportional to rCBF holds, almost any number of arbitrary ways of manipulating the early rising bolus concentration time course would yield the same rCBF result. One can equally use the value of the second derivative, the 3rd derivative, and so on of the rising bolus concentration time course. While all these different methods yield theoretically the same flow results, one method could prove superior to another under specific conditions when noise is present. Eq (4), like Eq (2), is valid only if the AIF’s of the two pixels under comparison are the same, or if a global AIF is assumed as we have stated earlier for this manuscript. The question on the use of global AIF or local AIF is brought up in the Discussion section as the presence of bolus delay and dispersion between the artery and the tissue is an on-going concern of the research and clinical community.
If random noise is the only concern for our data, the signal intensity approach is in principle having better signal to noise than the slope approach. Due to propagation of error, a derivative evaluated by the simple difference between two time points would have lower signal to noise than the signal intensity alone. On the other hand, when physiological noise of a coherent character is present, the slope result which, being essentially a subtracted difference between subsequent data points, could be less sensitive to the coherent physiological noise oscillation common to the neighboring data points.
In theory, the signal intensity approach needs only one data point while the slope approach needs the difference of at least two neighborhood data points. However, a single data point or the difference of two data points suffers too much from signal fluctuation when the bolus time course is noisy (not just around TOA) and could return inaccurate values for rCBF. In practice, for the signal intensity approach, we used the running average of four time points and called it the signal averaging approach throughout the rest of this manuscript. For the slope approach, we derived a slope using Matlab (MathWorks Inc., MA) based on the linear least squares fit of the same four time points. Demonstration of consistency between rCBF results obtained by two different methods -- the slope method and the signal averaging method – will aid our confidence that the bolus data meet the demand of ET’s basic assumption. With simulation data, flow ratio results from the slope and the signal averaging methods can be checked for their accuracy against the true flow ratio. With experimental data when the true flow ratio is unknown, the flow results obtained by these two different methods can be checked for their consistency with each other.
It should be mentioned that the initial slope method bears resemblance to the maximum gradient-based method (Kimura and Kusahara, 2009; Koenig et al., 1998; Miles and Griffiths, 2003) used in Computed Tomography (CT) and MR perfusion analysis. We leave a more detailed comparison between the maximum gradient-based method and ET to the Discussion section.
For noiseless simulation data, accurate flow ratio calculation is in principle obtainable immediately at TOA. But a very useful property of ET’s basic assumption predicts that flow ratio does not need to be calculated at TOA alone. Indeed, for noiseless data, identical flow ratio values are expected to be obtainable over the whole range going from TOA to TOA+τ, with τ being the time window before the contrast agent starts washing out. For noisy simulation data, useful flow result is not expected to be obtainable over exactly the same range from TOA to TOA+τ. The reason is that inaccuracy of flow ratio is to be expected in the neighborhood TOA even if true TOA is known. However, we do expect the calculated flow ratio at some distance away from TOA (at the higher CNR portion of the bolus time course) to eventually settle into a steady state time course of fluctuating around a mean flow ratio value. This time segment of steady state flow ratio is named FRSS for convenience. Basically, we want to identify the range of FRSS over which ET is expected to yield a steady state mean flow ratio for noisy data. Given the interference of noise at TOA, the beginning of FRSS for noisy data is hypothesized to be TOA+toffset, with toffset being a time delay allowing a more accurate evaluation of rCBF to take place at a better CNR portion of the rising contrast concentration time course. The end of FRSS is to be investigated but it should mark operationally the time when the bolus concentration time course reaches beyond τ.
The existence of FRSS is expected for the simulation data if the chosen model R(t) allows the early time points to meet ET’s basic assumption. For experimental data where R(t) can only be surmised, a successful identification of FRSS will booster the claim that real world data can be compatible with ET’s basic assumption.
With experimental MRI data, we checked if rCBF results calculated by ET results fell within the published range of rCBF values collected by PET and MRI. Examples of published values can be found at Leender et al. and Ostergaard et al. (Leenders et al., 1990; Ostergaard et al., 1996a).
Our workflow included analyzing simulation data to examine the ideal behavior of ET and investigating experimental MRI data to see how well they matched the expected ideal ET behavior. In collecting experimental data, we imaged an anesthetized macaque in a clinical 3 Tesla MR scanner. While human clinical data would have been ideal, the untested ET technology and the lack of adequate brain coverage at short TR acquisition were some of the reasons why human patient studies are still on hold. Meanwhile, non-human primates were suitable subjects to provide an important model of the human brain. After presenting results to validate ET, we proposed an ET menu to account for the effect of noise and identify a useful range of early time points for rCBF calculation.
In our study, two categories of perfusion data were collected: data simulated with Monte Carlo method and MRI data acquired on monkey.
Tissue signal curves were simulated at TR=300ms with an AIF modeled by a gamma variate (Calamante et al., 2000) using techniques previously described by Wu et al. (Wu et al., 2003b). We performed simulations with and without noise added, generating bolus concentration time courses which reflect multiple flow conditions. There were seven levels (named level 1 to level 7) of flow value varying between 10–70 ml/100g/min in 10 ml/100g/min increments. Each flow level was given transit delays ranging from −5 to 5 s with respect to a given AIF. To avoid the need for the interpolation of data with discrete TR of 300ms, we selected to study only four steps of transit delays (−5 s, −2 s, 1 s, 4 s) in 3 seconds increments. The same blood volume was used for each flow level giving 7 τ values ranging from 3.4 sec to 23.8 sec. We investigated if the flow ratio calculated by ET would correctly match the true simulated flow ratio under a multitude of conditions. For ET, we selected a box-shaped or rectangle function model (Buxton, 2002) of the residue function R(t)) (Buxton, 2002; Ostergaard et al., 1996b). With a rectangle function R(t), no contrast agent leaves the tissue before τ. A rectangle function R(t) is an excellent model to be used for the clarification of the ideal behavior of ET even if the R(t) for real world data is a mixture of R(t) of many shapes (Buxton, 2002). Noiseless data were useful to separate the effect of noise from the ideal ET behavior. Signal-to-noise ratio (SNR) of the noisy signal curves were set to 20 modeled after the low SNR clinical imaging environment. The noisy simulation data set had not been modeled for physiological noise or recalibrated based the result of monkey’s imaging data which has an SNR of around 100 obtained by dividing the mean MR signal of the brain image by the standard deviation of the noise signal of a region of interest (ROI) outside the brain.
All monkey imaging experiments were performed according to NIH and the institutional animal care guidelines. We imaged the same healthy macaque (Macacca mulatta, 8 kg weight) in two different MRI sessions on a 3 Tesla Siemens Trio scanner (Siemens Medical System, Erlangen, Germany) using a customized surface coil. The monkey was anesthetized with ketamine/xylazin, held in place by a MR compatible head frame. In the first imaging session, a routine Gd-DTPA bolus imaging was performed. The monkey was imaged by gradient echo (GE) echo planar imaging (EPI) with TR=300ms and TE=31ms to acquire six slices with 600 time points lasting for 3 min. Gd-DTPA (0.2mmol/kg) was injected via leg vein at around the 250th image at 75 s and then flushed with 12ml of saline at an injection rate of 2ml/sec by a power injector. In the second imaging session, the animal received Gd-DTPA twice, about 30 minutes apart. The first Gd-DTPA injection was imaged by GE EPI, while the second injection was imaged by Simultaneous Echo Refocusing (SER) EPI. SER EPI delivered a larger brain coverage than GE EPI at the same TR of 300ms by offering a simultaneous acquisition of two slices with each EPI readout, obtaining a total of 10 imaging slices for each TR. Spatial resolution of both GE EPI and SER EPI was 3mm × 3mm × 3mm.
The same TOA search algorithm is used for both simulation and experimental MRI data. Since a vast number of different approaches can potentially be considered for TOA estimation, no attempt has been made here to identify an optimized algorithm. We used an algorithm which searches and defines TOA with the intersection of two lines. The first line was fitted over the baseline and the second line fitted over the early rising part of the bolus. Fitting was done using Matlab. Fig 2a gives a pictorial representation of how the algorithm was applied to the experimental MRI data. It worked the same way with the simulation data.
For the baseline linear fit, we collected a long baseline of more than 70 seconds and used around 50 seconds of data at the latter part of the baseline to make a linear least squares line fit. Based on the shape of the MR signal time course taken from a small region next to a major artery, we had an approximate idea when the contrast agent would first appear. So data points used for the baseline fit could include time points reasonably close to the time of the first appearance of the bolus signal. The second line fit for the rising bolus was taken with the consideration of the roughly biphasic appearance of the rising portion of most Gd-DTPA bolus concentration time courses. The biphasic appearance suggests a large second derivative at the transition between the roughly two segments of the rising bolus curve. We therefore identified the time location to reach the maximum second derivative of the rising concentration time course (Fig 2a). Then we approximated the early rising portion of the ΔR2* curve with a line fit using the concentration time point at the time location of the maximum second derivative and three preceding concentration data points. The location of the maximum second derivative of the rising bolus curve could easily be located by automated calculation. Pixel by pixel, the intersection of these two fitted lines (Fig 2a) was taken to be the TOA. Calculated time of arrival (TOA) map (Fig 2b) for MRI data confirmed that the contrast agent arrival time was different for each pixel. Arrival time was earliest for pixels in the neighborhood of arteries and in the gray matter.
The TOA derived from our algorithm was compared with true TOA in simulation data to assess the usefulness of our TOA searching algorithm and the possible effect of any miscalculated TOA.
Signal averaging and slope procedures written in Matlab were used for both simulation and experimental data. The signal averaging approach took the running average of a group of four early time points of the rising bolus concentration time course. The slope approach took the slope of the linear least squares fit of the same four early time points. Any group of four data points can be identified by the time location of earliest time point of the group (e.g. 3 s following TOA). We always compared the signal averaging and the slope results to see if they yielded similar rCBF values. Transit delays have always been corrected for early time points used for slope or signal averaging.
If the correct transit delay is found and properly accounted for each bolus time course, it is intuitively obvious that flow ratio results calculated by ET would be independent of transit delays. Nonetheless, given the knowledge of the four given steps of transit delays built into the simulation data, we double checked the claim of independence from transit delays by measuring the signal averaging and the slope values of the bolus time courses.
In order to find FRSS, we first made a time course of the ratio of two bolus time curves, already corrected for transit delays, from pixels of different flow levels. The two bolus time courses could take the form of contrast concentration signal intensity or its derivative (it would be average signal and smoothed derivative in practice). From the ratio of the two bolus time courses we then proceeded to identify FRSS, namely the portion of the ratio time course that maintains a steady state flow ratio.
There is a number of important parameters of FRSS calculated with respect to true TOA when it comes to noisy simulation data. Such parameters included the time TOA+toffset at the beginning of FRSS, the duration of FRss and the mean value of the data points of FRSS. The duration of FRSS was named τss for both the noiseless and noisy data. For noiseless data, τss would be equivalent to τ and the end of FRSS would be at TOA+τss = TOA+τ. For noisy simulation data, the end of FRSS was operationally marked by the time when the flow ratio no longer maintained a steady state status. The end of FRSS should indicate the practical upper time limit useful for rCBF evaluation. The calculated mean value of the data points of FRSS was basically the ET result of calculated flow ratio, to be used for the comparison with the true flow ratio.
It should be pointed out that in the correction for transit delays, true TOA was used for the simulation data analysis and calculated TOA was used for the experimental data analysis. The reason that true TOA instead of calculated TOA was used for simulation data was to find out the important information of the effect of noise on the various ET properties when the true TOA was known. It is important not to confound the problem of miscalculated TOA, which was one of many consequences of the noise effect, with other major consequences of the noise effect. The effect of miscalculated TOA was handled separately in the data analysis section above under “TOA search” with TOA derived from our chosen algorithm being compared with true TOA.
For the noisy simulation data, twenty-one true flow ratios could be obtained out of combination of the seven given flow levels (level 1–level 7, short for 10–70 ml/100g/min in 10 ml/100g/min increments). Three calculated flow ratios (level 7/level 1, level 3/level 1, level 7/level 5), presented in Table 1, were considered sufficiently representative for the demonstration of how the presence of noise could affect the value of toffset, the value of τss, and the calculated mean value of the data points of FRSS.
For experimental data, flow ratio values were taken from the regions of interest (ROI) of caudate gray and white matter (Fig 3) for the purpose of identifying FRSS.
For the experimental data, bolus data were preprocessed in the manner different from the simulation data. We took advantage of our short TR=300ms and applied a low pass filter of 0.16Hz to our experimental data to reduce respiratory noise. Care was taken to reduce respiratory noise without affecting unduly the higher frequency components of the bolus concentration time course. As any smoothing could modify the noise distribution, a more rigorous investigation of the smoothing process is warranted in the future.
Even after the respiratory noise filtering, there was left-over physiological noise which was time-coherent and slow varying (Biswal et al., 1995). Such coherent physiological oscillation could be observed clearly at the baseline of the bolus time course of Fig 2a. Just like random noise in simulation data, slow varying physiological noise would affect the accuracy of flow calculation even if the true TOA were known. The oscillating effect of the slow varying physiological noise at TOA could artificially reduce or increase the MR signal value of the arrival point of the bolus signal, resulting in an inflation or deflation of the true flow ratio near TOA. Fortunately, unlike random noise which is always present, cyclic and time coherent physiological noise could theoretically be managed and reduced. While there is virtually an industry on physiological noise reduction techniques for fMRI literature alone, it is beyond our scope to examine the large variety of methods here. We proposed instead a very simple step which is adding a compensating constant value to the bolus concentration time courses in order to manage slow-varying noise effect. Adding a constant served to shift the oscillating baseline above the zero value, greatly reducing the chance of the oscillating signal near TOA obtaining a false near-zero value which in turn could drive a calculated flow ratio to artificially and abnormally high value. The compensating constant would be small because the physiological fluctuation was usually in the order of a few percentage of the baseline MR signal. The size of the compensating constant, with an opposite sign to what it was used to compensate for, would be somewhere between zero and the maximum amplitude of the baseline physiological oscillation immediately preceding TOA. Since an added constant value becomes a smaller and smaller part of the bolus signal amplitude the further away the rising bolus is from the time of TOA, the flow ratio value is expected to be modified at the beginning of the bolus arrival but not at the later time, thereby providing correction to any distortion (inflation or deflation) of the flow ratio value only near TOA. The compensating mechanism was considered an ad-hoc manipulation because the compensating signal amplitude remains a somewhat subjective process. The compensating constant modified the effect of physiological noise but not that of random noise which, not being time coherent, was little affected by the addition of a constant. We expected that the signal averaging results, but not the slope results, would be affected by the ad-hoc manipulation because derivatives depended on the difference of signals and would not be affected by an added constant.
We made rCBF maps for experimental data, using both signal averaging and the slope approaches to check if calculated ET results fell within the range of literature values of rCBF. To stay relatively clear of the influence of noise, all the representative MRI maps in this manuscript were made with data taken from the time location of 3 s following TOA. Explanation and justification of the choice of 3 s come from results to be presented in the Results section.
For quick reference purposes, we re-analyzed our data with the deconvolution method. Perfusion calculation by deconvolution depends critically on factors like signal to noise of the data, correction to transit delay, adjustment of the free parameter of SVD threshold, image acquisition rate, etc. (Calamante et al., 2000; Carpenter et al., 2006; Ostergaard et al., 1996a; Ostergaard et al., 1996b; van Osch et al., 2001; van Osch et al., 2003; Wu et al., 2003a; Wu et al., 2003b). Hence, full comparisons between ET and deconvolution require extensive modeling and optimization of deconvolution and data parameters based specifically on the detailed characteristics of monkey data and will be reported elsewhere.
To get a better 2-dimensional view of the hypothesized match between the rCBF results made by the signal averaging and the slope approaches, we made a ratio map of the two rCBF maps generated by these two different ET calculating methods. The idea was that a ratio map should be devoid of perfusion patterns if the rCBF maps made by the two different calculating methods were similar.
Lastly, we made rCBF maps of 10 slices of SER EPI data using both signal averaging and slope approaches and then examined the gray-white flow contrast in comparison with those from the corresponding gray-white regions of interest of the regular 6-slice GE EPI data.
Topic I-IV raised in the Introduction section outlined the various conditions and limitations that ET needed to work under. When those conditions and limitations were taken into account, the flow ratio results calculated by ET not only matched the true flow ratio well for simulation data, but were also consistent with flow ratio values reported in the literature for experimental data. These matching flow ratio results increased the confidence on the compatibility of the real world bolus data to ET’s basic assumption.
Since there were three different data sets (noiseless simulation data, noisy simulation data and experimental data) and each data set was studied to obtain somewhat different type of information for ET, we offer a general description of results shared by all three data sets first and then details of the results special for each data set.
The useful range of early time points for rCBF calculation (Topic III) was identified to be FRSS, a time segment of steady state flow ratio for both simulation (Fig 1, Table 1) and experimental data (Fig 4). For noiseless data, no time delay was found as expected between TOA and the beginning time of FRSS (Fig 1a) and the duration of FRSS (τss) was demonstrated (Fig 1a) to match the given τ value of 3.4 s for a flow level of 70ml/100g/min. Noisy data, both simulated and experimental, showed that (Topic I) the accuracy of flow ratio calculation was compromised near TOA (Table 1, Fig 1b, Fig 4) and a time delay of toffset was required to be added to TOA before FRSS, the steady state time segment, could be reached. The size of toffset was shown to depend on CNR of the flow ratio near TOA (Table 1).
The upper time limit of the useful range of early time points was shown to be TOA+ toffset +τss (Table 1, Figs 1 and and4).4). After TOA+ toffset +τss, the flow ratio value started to decline slowly (Figs 1 and and4)4) for the signal averaging method but dropped rapidly for the slope method as the derivative went negative after the peak of the bolus. The values of τss of noisy and noiseless simulation data (Table 1, Fig 1) demonstrated a surprisingly close correspondence despite the difference in toffset.
Details more unique to each individual data set are presented below.
Results from noiseless simulation data clarified and confirmed the ideal behavior of ET for rCBF evaluation.
With the four steps of transit delays (in 3 second increments) given in simulation data, identical calculated flow ratio values were obtained by ET for each of the four steps demonstrating that ET results were independent of transit delays. For example, any calculated flow ratio at one particular transit delay could be reproduced at the other three steps of transit delay.
Given our TOA searching algorithm of two intersected lines, the same calculated TOA value was returned by ET for all seven flow levels of simulation data. But the calculated TOA value was delayed from the true TOA value by one time point (300ms). Why was our calculated TOA different from the true TOA, even for noiseless data? Our algorithm of two intersected lines did not reproduce the true TOA because the gamma variate signal of the simulated data was not a linear function so a linear line fit was unable to fit the rising bolus time course perfectly to locate the exact position of the true TOA. Despite this consideration, what is important to recognize is that the calculated TOA was different from the true TOA by the same amount (of one time point) for every flow level. Since the calculated TOA held an identical time relationship with the true TOA for all different flow levels, the consequence was that the same ET results would be obtained regardless whether calculated TOA or true TOA was used.
Results of calculated flow ratio from noisy simulation data were reasonably consistent with the true flow ratio, starting at a time delay of toffset after TOA. Fig 1b offers a visual representation of the example of a single calculated flow ratio corresponding to the true flow ratio of level 7/level 1 (level 1–7 is equivalent to 10–70ml/100g/min). Three calculated flow ratio results (corresponding to true flow ratios of level 7/level 1, level 3/level 1 and level 7/level 5) from the total of 21 possible calculated flow ratios were presented in Table 1 which reported the results of the following parameters: toffset, the mean value of the data points of FRSS, and τss. The numerical values of FRSS, toffset, and τss were all presented in Table 1 so they would in general not be repeated in the text except when clarification of the text was necessary. The calculated mean values of the FRSS (e.g. 6.6±0.14, 2.9±0.1, 1.4±0.02 for average signal) were consistent with the values of the three chosen true flow ratios (7, 3, 1.4). The calculated flow ratio (FRSS) of 6.6 did appear to be somewhat lower than the true flow ratio of 7. However, the original published deconvolution result reported in (Wu et al., 2003b) also showed that the mean flow result of the highest flow level calculated by the deconvolution method was 95% of the true flow level of 70ml/100g/min when a rectangle function R(t) was used, consistent with our mean FRSS value of 6.6 which was also 95% of the true flow ratio of 7.
It is clear from Table 1 that the CNR of the calculated flow ratio had a major effect on the size of toffset and the accuracy of FRSS in comparison with the true flow ratio. The worse CNR of the flow ratio, the larger toffset became (e.g. toffset of level 7/level 1 was larger than toffset of level 7/level 5). The better CNR of the flow ratio, the closer was the mean FRSS value to the true flow ratio (e.g. FRSS values of 6.6±0.14, 2.9±0.1 and 1.4±0.02 in comparison with the true flow ratios of 7, 3 and 1.4).
The duration of FRSS, which is represented by τss, was reasonably close to the true τ value of each of the three chosen flow ratios. It is interesting for the value of τss to be so close to that of τ. The implication was that for noisy data, the time duration of FRSS was not simply reduced by the size of toffset. Sample results from Table 1 are illustrative. Even though toffset was 2.4 s (row 1 and 2 of Table 1), the slope results gave τss = 3 s vs. τ = 3.4 s and τss = 7.2 s vs. τ = 7.9 s and the signal averaging results gave τss = 3 s vs. τ = 3.4 s and τss = 8.1 s vs. τ = 7.9 s. This intriguing observation was further addressed in the Discussion section.
Similar to noiseless data, flow ratio values obtained by ET were independent of the four given steps of transit delays.
With respect to our TOA searching algorithm, the same mean calculated TOA value, delayed from the true TOA value by four time points, was returned by ET for 6 out of 7 flow levels. Only the lowest flow level (10ml/100g/min) with the smallest CNR near TOA showed an extra delay in calculated TOA (delayed from true TOA by five time points instead of four). Since calculated TOA values of 6 out of 7 flow levels (20–70ml/100g/min) held the same fixed time relationship with true TOA, the same ET results would be obtained for these flow levels regardless whether calculated TOA or true TOA was used. The miscalculated and delayed TOA value for the lowest flow level would lead to an overestimate the bolus concentration value of the lowest flow level and therefore an underestimation of flow ratio.
With experimental data, the average results of the gray-white flow ratio from the ROI’s of caudate gray and white matter (Figs 3–5) taken from two different GE EPI imaging sessions were 3.2±0.2 for signal averaging and 3.2±0.3 for slope. Those values fell within the range of between 2.5 and 4 reported in the literature (Leenders et al., 1990; Ostergaard et al., 1996a). A certain amount of reproducibility could be inferred from the comparable ET results from these two separate imaging sessions even though the sample size was small. For reference purpose, the deconvolution results of flow ratio from the same gray-white ROI’s was 3.1±0.2.
The rCBF outcome generated by the slope and the signal averaging methods was comparable and further confirmed by the presentation of the ratio map (Fig. 6) of the slope result and the signal averaging result. Little perfusion contrast was observed in the ratio map except for pixel outliers which were either vascular pixels with enormous flow or pixels that showed partial volume mixture of tissues with vessels (van Osch et al., 2001). To go beyond a visual impression of the ratio map, an estimate of the degree of flatness of gray-white contrast in the ratio map could be assessed by the relative spread of the values of all the pixels in the brain volume of the ratio map. The relative spread of the values of all the pixels was expressed as Δy/y where Δy was the standard deviation and y was the mean of all the pixel values. A ratio map looking perfectly flat would show Δy/y = 0. We obtained a gray-white contrast with Δy/y = 0.07 for this ratio map as opposed to the gray-white contrast with Δy/y = 0.75 for the original rCBF map, suggesting that the ratio map of flow results calculated by the two different methods appeared relatively flat. It was also unlikely for Δy/y to merely reflect the random noise of the ratio map. If it were just random noise, Δy/y for the ratio map would always be larger than Δy/y for each one of the original maps.
By following the features of the time course of the flow ratio of experimental data (Fig 4a,b), it could be observed that the noise effect on flow ratio results showed an interesting difference between real world data and simulation data at the neighborhood of TOA. Fig 4a and 4b presented similar flow ratio results from two MR imaging sessions for reproducibility check-up. Unlike the observation of a leisurely rising flow ratio from the baseline of the noisy simulation data, a huge overshoot of flow ratio value near the beginning of TOA was demonstrated only for the signal averaging result of the experimental data. The huge overshoot in flow ratio value lead to a lengthy toffset before FRSS was reached. While toffset was large for the signal averaging result, it was quite small for the slope result. The wide gap in toffset values was a consequence of the much larger flow ratio overshoot at TOA for the signal averaging result. The big difference in flow ratio results between the signal averaging and slope approaches started to narrow and get close to each other at 1.8 s following TOA (Fig 4a,b) and remained close together forming the steady state of FRSS until 5.1 s when the slope signal began to drop off. Charting the behavior FRSS provided us the rationale to select data from a time location at 3 s following TOA to make all MRI maps.
Such flow ratio results at the beginning of TOA confirm our presentation in the Introduction section that the signal averaging result could be sensitive to coherent physiological noise oscillation while the slope result is not.
The ad-hoc manipulation was proposed with the aim to ameliorate the deleterious effect of physiological noise. When the ad-hoc manipulation was applied, the overshoot of the flow ratio at TOA of the average signal data was greatly reduced (Fig 4a,b) and the magnitudes of flow ratio and toffset near TOA were now almost the same for the signal averaging and the slope approaches. Fig 4a and 4b also confirm that the ad-hoc manipulation affects the later part of the flow ratio results much less than the early part. Adding a constant value to the bolus concentration time course did not change the slope value at all because the derivative of a constant is zero. So the ad-hoc manipulation did not do anything to the flow ratio obtained by the slope approach.
Finally, the gray-white flow ratio made by ET for the 10-slice SER EPI data (Fig 7) at the corresponding caudate and white matter ROI’s (slice 4) was 3.1 (signal averaging), and 3.2 (slope), consistent with ET results from 6-slice GE EPI data. Even though only the average signal map for rCBF was displayed in Fig 7 to reduce repetitious display of similar figures, the calculated flow ratio results confirmed once again the comparability of the signal averaging and the slope methods. Reproducibility of rCBF results calculated by ET was enhanced by the reasonable match between flow ratio results from the single imaging session of SER EPI and those results already presented for the two imaging sessions of GE EPI.
With our results showing that calculated flow ratios by ET were consistent with perfusion values obtained from simulation and experimental data, a reasonable recipe for rCBF calculation by ET can be offered as follows: 1) Apply a TOA searching algorithm to calculate TOA for each pixel. 2) Use the value of TOA to make the correction for transit delay for the bolus time course. 3) Choose the signal intensity (or signal average) or the time derivative (or smoothed slope) of the rising bolus concentration time course for rCBF evaluation. 4) Take a ratio of concentration time courses or of the derivative time courses selected from chosen gray and white matter ROI’s. The two bolus time courses could take the form of contrast concentration signal intensity or its derivative (average signal or smoothed derivative in practice). 5) Find the useful range of early time points for rCBF calculation by identifying the steady state flow ratio time segment (FRSS) from the ratio of the two bolus time courses. The range of FRSS give an idea on which early time points are far enough away from the influence the noisy neighborhood of TOA but still fall within the time window of τ. 6) Select time points from FRSS to make rCBF maps pixel by pixel.
In short, this recipe identifies and then manipulates the portion of early time points useful for rCBF calculation.
In order to gain a better understanding of some of the fine points and implications of the recipe for ET, we proceed to discuss below the question on global AIF vs. local AIF’s as well as questions related to the four topics raised in the Introduction section.
If there are local AIF’s due to delay and dispersion, as it had been reported in research literature (Calamante et al., 2000; Lorenz et al., 2006a; Lorenz et al., 2006b), there would be errors in rCBF evaluation by ET, in the same manner that errors would be expected in the deconvolution method using a single global AIF to estimate rCBF. Despite the many models provided in the deconvolution literature (Calamante, 2005; Calamante et al., 2000; Lorenz et al., 2006a; Lorenz et al., 2006b; Willats et al., 2006, 2008) to estimate the errors caused by the problem of dispersion and the presence of local AIF’s, the problem of quantitatively and definitively identifying local AIF’s at each pixel or region of the brain remain a technical challenge. Even when a particular local AIF model was selected for perfusion analysis, there were disagreements over whether a local AIF approach is necessarily superior to a global AIF approach. For example, a 2009 paper by Christensen et al (Christensen et al., 2009) showed that currently used local AIF methods did not outperform conventional whole-brain AIF approaches in a number of perfusion parameters which best predict tissue infarction with a population of more than 90 stroke patients. In a presentation by Stankiewicz et al at the most recent 2010 ISMRM conference meeting at Stockholm, the reproducibility between two baseline MRI scans separated by 48 hours for 13 tumor patients appeared better for global AIF maps than local AIF CBF maps. Hence, it is reasonable for us to develop our understanding of ET using the global AIF assumption and not be distracted by going deeply into this advanced and complicated topic of local AIF which will be left to future research.
While our TOA search algorithm of two intersected fitted lines did not reproduce the true TOA in noiseless simulation data, it remains useful because the calculated TOA holds a fixed time relationship with the true TOA for all different flow levels. Besides, the shape of the gamma variate function of the simulation data is not necessarily a perfect model for experimental data, so the difference obtained between calculated and true TOA might be less of a problem for real world data. It does suggest that the particular shape of the bolus time course needs to be taken into account on the choice of the TOA searching algorithm. For noisy simulation data, the calculated TOA was off with the lowest flow level of 10mg/100g/min, suggesting that the test for any future TOA searching algorithms will be scored on their ability to handle very low flow. Another lesson from our algorithm was the usefulness of anchoring the search for TOA with data points at some distance away (location of the maximum second derivative in our case) from the low CNR region of the immediately neighborhood of TOA. Continued research to better identify TOA will be a major direction for ET.
Signal averaging and slope methods yielded comparable rCBF results over all the different data sets we have studied, offering evidence that MR contrast-enhanced data was consistent with ET’s basic assumption. Other than their similar rCBF outputs, we have also looked into a number of differences between these two approaches. For example, beyond the end of FRSS, the value of signal intensity drops far more slowly than the value of slope which quickly drops below zero as is required by the bolus time course starting to turn direction. So signal intensity in theory allows some time beyond τ for an estimation of an approximate (even when ET’s basic assumption begins to be violated) flow ratio. We stated in the Introduction section that signal intensity had inherently better SNR than slope if we only considered random noise. The SNR report of different simulated flow ratios at Table 1 supported that statement. On the other hand, we also showed that the signal intensity results were a lot more sensitive to the oscillation of physiological noise near TOA than the slope results. If we prefer not to try some ad-hoc method to manage the effect of physiological noise, using the slope approach to calculate rCBF might sometimes be a better option for data analysis near TOA than using the signal averaging approach.
There is considerable interest to compare the maximum gradient-based method (Kimura and Kusahara, 2009; Koenig et al., 1998; Miles and Griffiths, 2003) and ET. If the time taken to get to the maximum gradient is always within the time window of τ, then the maximum gradient-based method would indeed be theoretically compatible with ET. The maximum gradient-based method is independent of TOA and is located at the better CNR portion of the bolus time course. Despite these advantages, there is however no a priori reason why the time taken to reach the maximum gradient always occurs within the time window of τ. After all, the time taken to reach the maximum gradient is dependent on factors such as the bolus volume and the rate of injection and the patient’s cardiac output (Miles and Griffiths, 2003). Contrast agent washout violating ET’s basic assumption of no contrast agent departure had been reported for the maximum gradient condition (Miles and Griffiths, 2003), resulting in significant drop of estimated perfusion values. On the other hand, the data utilized by the initial slope method is more likely to fall within the time window of τ. It is also possible to construct a simulation data set where the time taken to reach the maximum gradient is larger than τ while the data points selected for the initial slope are still within the τ window. The simulation data set we obtained from Wu et al (Wu et al., 2003b) included a bolus concentration time course with the time to reach the maximum gradient being larger than τ. Under that simulated condition, the flow ratio obtained by the maximum gradient-based method turned out to be lower than the true simulated flow ratio.
With our interest centering on gaining an understanding on some detail behavior of early time points for rCBF evaluation, including the time points immediately after the arrival of the contrast agent, the choice of focusing on two approaches - the initial slope method and the signal averaging method – is more appropriate for our purpose. Given the somewhat higher uncertainty than the initial slope method of not meeting ET’s basic assumption, the maximum gradient-based method deserves more exploration in the future but will be mentioned here mainly for reference purpose.
The excellent match (for noiseless simulation data) and the good match (for noisy simulation data) between calculated flow ratio by ET and true flow ratio showed that ET could be useful for rCBF calculation.
The length of the FRSS, namely τss, of noisy data was comparable to the length of τ and appeared to be not too affected by the amount of toffset added to TOA. The reason for the relative closeness between τss and τ in noisy data remains to be clarified. One possibility was that beyond τ, the loss of perfusion accuracy by the higher flow level (numerator of the flow ratio) was compensated by the gain in CNR of lower flow level (denominator of the flow ratio), accounting for the good correspondence between τss and τ. In other words, we hypothesize that while the numerator of the flow ratio starts violating the ET assumption leading to a worse estimation of the true flow ratio, the continual improvement of CNR of the denominator due to the rising bolus time course enhances significantly the accuracy of flow ratio estimation, providing some kind of compensating effect.
Just like noisy simulation data, the duration of FRSS, namely τss, for experimental data, was presumably related to τ of the chosen gray matter ROI.
In real world practice, we would like to get a single toffset and a single τss to represent the useful range of FRSS for all tissue types and not just between a pair of tissue regions. Thus, we would want to identify a toffset from an ROI with the lowest CNR and that would be white matter, and a τss taken from an ROI with the smallest τ, and that would be gray matter. That is what we did to obtain flow ratio data for Fig 4 to answer the questions on FRSS.
The proposed example of the somewhat arbitrary ad-hoc manipulation was used mainly to highlight the possibility and the importance for more research to manage physiological noise.
The match between calculated gray-white flow contrast for ET and published values of perfusion in the literature reinforced the usefulness of ET for perfusion measurement.
While we defer a more quantitative evaluation of advantages vs. disadvantage between ET and the deconvolution approach to future research, we are listing below a few of the hypotheses to highlight some of the possible reasons to consider ET.
The calculation of rCBF by ET and by the deconvolution method is fundamentally employing the same parameter for rCBF evaluation, as it would be expected. ET works in the regime of residue function R(t) = 1. The deconvolution method works with the height of the maximum amplitude of the deconvolved residue function, which in the ideal case is the same as looking for R(t)=1.
Hence, the early time points carry more weight for rCBF calculation than the later time points of the bolus curve. While more investigation and verification is needed, it is reasonable to hypothesize that the latter part of the full bolus time course used for the deconvolution process, for example the bolus recirculation components, may even introduce more noise for rCBF calculation. One can try to exclude the recirculation components by an arbitrary cut-off. But an arbitrary cut-off of the data could lead to more uncertainty in rCBF estimation by deconvolution. In addition, the deconvolution process itself may introduce noise which is not a question for the ET method.
There is interest in investigating whether the application of ET in leaky tissue could lead to better rCBF estimation than the deconvolution method. With leakage, two effects can impact the estimate of the concentration of the contrast agent and subsequently of rCBF. First, there is a continuous loss of concentration due to leakage which may affect ET and deconvolution in a similar way. There is a second effect that could affect the deconvolution method more. The additional T1 contrast from the leakage of contrast agent, mixed with the change of T2* contrast due to loss of concentration, introduces more unknowns and uncertainty to the total MR signal for the estimation of the concentration of the agent. We hypothesize that since the relative weight of the T1 contrast due to leakage is larger at the latter part of the bolus time course than at the early part, it will make a bigger impact on the deconvolution method which utilizes the full time course.
While the short TR of 300ms was valuable for the initial investigation of ET’s properties and the SER type of MR sequences helped to resolve the brain coverage problem of short TR, it would still be helpful to be able to apply longer TR in practice. The feasibility of obtaining reasonable clinical ET results with longer TR remains to be explored. On the other hand, there might be another unanticipated advantage of using shorter TR. A routine clinical TR of 1.5 s cannot, by Nyquist consideration, be used to properly estimate a τ value shorter than 3 s. Hence a TR of 1.5 s favors τ results longer than 3 s and could introduce unintentionally a bias into the reported data base of τ values in the MRI literature.
Future research would explore imaging in a combined PET/MRI scanner, applying SER EPI with parallel imaging (Reese et al., 2009) to reduce image distortion problems (Reese et al., 2002), testing other single-shot multi-slice imaging techniques (Harvey and Mansfield, 1990; Harvey and Mansfield, 1996; Larkman et al., 2001; Mansfield and Howseman, 1989; Song et al., 1994) and studying animals with compromised hemodynamics or extravasated contrast agent under conditions such as stroke or tumors.
Simulation data were used successfully to examine the ideal behavior of ET while monkey’s MRI results offered encouraging support to the utility of ET for rCBF calculation. We have shown that ET made rCBF maps of reasonable gray-white flow contrast which matched perfusion values in the literature. Our results supported the compatibility of experimental data with ET’s basic assumption. A better coverage of the brain for clinical consideration was obtainable by applying SER EPI. A recipe to run ET was presented, with attention paid to the noise problem around the time of arrival (TOA) of the contrast agent.
We would like to acknowledge the support of the Athinoula A. Martinos Center for Biomedical Imaging. This study was supported in part by the National Institutes of Health (grant R01NS59775).
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.