|Home | About | Journals | Submit | Contact Us | Français|
Time of arrival (TOA) of a bolus of contrast agent to the tissue voxel is a reference time point critical for the Early Time Points Perfusion Imaging Method (ET) to make relative cerebral blood flow (rCBF) maps. Due to the low contrast to noise (CNR) condition at TOA, other useful reference time points known as relative time of arrival data points (rTOA) are investigated. Candidate rTOA's include the time to reach the maximum derivative, the maximum second derivative, and the maximum fractional derivative. Each rTOA retains the same relative time distance from TOA for all tissue flow levels provided that ET's basic assumption is met, namely, no contrast agent has a chance to leave the tissue before the time of rTOA. The ET's framework insures that rCBF estimates by different orders of the derivative are theoretically equivalent to each other and monkey perfusion imaging results supported the theory. In rCBF estimation, maximum values of higher order fractional derivatives may be used to replace the maximum derivative which runs a higher risk of violating ET's assumption. Using the maximum values of the derivative of orders ranging from 1 to 1.5 to 2, estimated rCBF results were found to demonstrate a gray-white matter ratio of approximately 3, a number consistent with flow ratio reported in the literature.
Dynamic imaging of a bolus of contrast agent passing through the brain's vascular system has been used to study blood flow in a variety of approaches (Axel, 1980; Herfkens et al., 1982; Miles, 1991; Ostergaard et al., 1996a; Ostergaard et al., 1996b; Rempp et al., 1994). In dynamic susceptibility contrast-enhanced MRI (DSC-MRI) of Gd-DTPA bolus injection, relative cerebral blood flow (rCBF) can be measured using the strategy called “Early Time Points” method (ET) (Buxton, 2002; Klocke et al., 2001; Kwong et al., 2011; Lee et al., 2004). The concept of ET, an idea close to the microsphere perfusion model, proposes 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 (Calamante et al., 2000; Ostergaard et al., 1996a). Given that basic assumption of before the contrast agent having a chance to clear from the tissue, the quantity of Gd-DTPA present in the tissue will be proportional to local CBF. Since the basic assumption will be cited repeatedly, it will be named the ET's basic assumption for the rest of the manuscript.
Some of the differences of ET in comparison with the deconvolution method (Ostergaard et al., 1996b) for perfusion analysis had been presented previously (Kwong et al., 2011) so we will just give a brief summary here. Correction of bolus transit delay is a built-in component of ET. Blood flow results can be obtained from the early time points before τ so ET will not be affected by noise from the later part of the bolus time course which includes recirculation. Whatever noise that may be introduced by the deconvolution process itself is not a question for ET. The evaluation of rCBF by ET for any single patient is independent of the arterial input function (AIF) if a global AIF is assumed, an assumption broadly useful in clinical perfusion imaging. For cross subject comparison of rCBF, a reference based method (Kimura and Kusahara, 2009) which sets a flow reference region of interest (ROI) on some referential tissue (e.g. white matter ROI) is one possible option for ET. When global AIF is not assumed, the local AIF question is technically challenging and beyond the scope of this study. The answer to cross subject comparison in ET and to some of the questions of local AIF is provided in a companion manuscript (Kwong and Chesler, 2011).
While ET is useful for measuring rCBF, it should be mentioned that absolute CBF can in principle be obtained by ET just like the deconvolution method of Ostergaard et al. if a proportionality constant, required for absolute CBF estimation, can be obtained independently and separately from a different technique like PET. Absolute blood flow might also be obtainable by ET based on the well studied microsphere model if arterial blood samples can be obtained from the appropriate arteries.
Having provided a summary of some of the features of ET, we will present now the focus of this manuscript: the identification of useful reference time points to allow ET to correct the bolus transit delay for the purpose of carrying out rCBF analysis.
The correction of the transit delay of the contrast agent bolus for each voxel is required for ET to make relative cerebral blood flow (rCBF) maps. Intuitively, the time of arrival (TOA) of the contrast agent at the voxel can be used to correct the transit delay of the bolus. But TOA brings its own concerns. While the search for accurate TOA is an ongoing research topic studied elsewhere (Kwong et al., 2011), TOA is recognized as being located at a low contrast-to-noise ratio (CNR) region of the bolus signal due to the small amount of contrast agent entering the voxel at TOA.
Fortunately, one does not need to know the time location of TOA of each voxel for ET to work. In this manuscript, we propose that relative TOA (rTOA) can also be used as a reference time point for ET. The requirement for a time location to be rTOA is that it retains the same time distance from TOA for every voxel. If rTOA can be found, it is equivalent to TOA as a reference time point for rCBF calculation by ET. Due to the low CNR location of TOA, there is a strong motivation for selecting rTOA at the high CNR region of the rising bolus concentration time curve. Our aim is to identify and investigate some of those candidate rTOA's. For convenience, the bolus concentration time course is named C(t). The signal intensity of C(t) is described in terms of ΔR2*, an MR susceptibility signal related unit which represents the contrast agent concentration.
A preferred candidate for rTOA would likely have the following properties: reasonably far away from the low CNR neighborhood of TOA, easy to evaluate, and retains the same time distance from TOA for every voxel without at the same time requiring the knowledge of the value of TOA. Those properties will be made clearer as our presentation proceeds. It turns out that there are many data points from the rising C(t) that are good candidates for rTOA.
Meeting ET's basic assumption is a requirement for rTOA, meaning that rTOA < τ. Since we do not usually have the knowledge of whether a particular early time point picked for rTOA meets ET's basic assumption in real world experimental data, we made the initial assumption in this manuscript that ET's basic assumption had been met by the rTOA candidates. Then we investigated the implication for those rTOA candidates when ET's basic assumption was violated.
Based on the requirements of reasonable CNR locations, simplicity in calculation and independence of knowledge of TOA, we will first propose a list of rTOA candidates with the justification to be given in later paragraphs. Relative to the unknown (and uncalculated) TOA, the first candidate is the time to reach the maximum derivative (TMD1) of the rising bolus time course. TMD1 is of course just the time position of the classical maximum gradient approach to evaluate rCBF (Kimura and Kusahara, 2009; Koenig et al., 1998; Miles and Griffiths, 2003). The second candidate is the time to reach the maximum second derivative (TMD2). It turns out there are other choices as well. Intuitively, if the time to reach the maximum fractional derivative (TMD-FR) exists between TMD1 and TMD2, we can consider a TMD-FR value between 1 and 2 as a third rTOA candidate. For convenience, the fractional derivative selected between 1 and 2 is an arbitrary 1.5 and the time to reach its maximum value is TMD1.5. There is a motivation for considering fractional derivative. TMD1 has good CNR because it is farther away from TOA but TMD1 runs a higher risk of being longer than τ for some high flow values, thereby violating the ET's basic assumption. TMD2 has minimal risk of being longer than τ but it is also at a relatively low CNR region of C(t). Noise amplification can also occur with the second derivative due to the propagation of error if derivatives are calculated by finite difference, namely, taking the difference of neighborhood data points (Appendix A). Since TMD-1.5 is chosen to be a time point between TMD1 and TMD2, it is expected be able to offer a good compromise. There is one more choice of rTOA to be considered. If good CNR is of paramount importance, it is also reasonable to examine lower order fractional derivative such as ½ with the time to reach the maximum derivative of the order ½ (TMD0.5) being the desirable rTOA. In choosing TMD0.5, the size of error generated by violation of ET's basic assumption, if there is any, needs to be considered to be tolerable for a particular research goal.
Beyond TMD0.5, the time to reach the maximum peak of C(t) will not meet ET's basic assumption as C(t) starts making a turn at its maximum peak. Had ET's basic assumption or the microsphere condition been always met, C(t) would keep on growing and then leveling off but would never turn back down to baseline. When C(t) of experimental data starts turning away from its peak signal towards baseline, it is an indication that the time of the peak of C(t) has gone beyond τ. That is why the approach used by earlier reports (Klocke et al., 2001; Lee et al., 2004) of first-pass magnetic resonance perfusion imaging which estimates rCBF using data including the peak of C(t) will not be explored here.
Attention needs to be paid to the properties of fractional derivatives (Kilbas et al., 2006; Oldham and Spanier, 1974; Samko et al., 1993). It is well known that the local property of the fractional derivative at a time point t involves boundary conditions and is not well defined like the local property of the integer derivative. However, TMD-FR can still be a candidate for rTOA as long as a local maximum of the fractional derivative can be found at TMD-FR. The important test is whether TMD-FR can retain the same time distance from TOA for every voxel. We will verify that rTOA property of TMD-FR with simulation data.
For convenience, we also set the following abbreviations: TMDX is the time to reach a maximum derivative of arbitrary order. MD1 is short for the maximum first derivative, MD2 for the maximum second derivative, and MD-FR for the maximum fractional derivative. MD1.5 and MD0.5 are specific for the maximum values of fractional derivative of 1.5 and of 0.5, respectively. MDX is the maximum value of a derivative of arbitrary order. D1 is short for the first derivative value and D1(t) is short for the first derivative time course. D2, D2(t), D1.5, D1.5(t), D0, D0.5(t) are similarly defined for the derivative orders of 2, 1.5, and 0.5. For derivative of an arbitrary order X, we also made the abbreviation DX for the derivative value and DX(t) for the time course. D-FR is short for fractional derivative.
Fig. 1 shows schematically that for TMDX, TMD0.5 > TMD1 > TMD1.5 > TMD2 for a simulated bolus concentration time course.
The reason why TMD1, TMD2 and TMD-FR are good candidates for rTOA is based on two claims. First claim: the relative time between TOA and each of these rTOA candidates is independent of bolus transit delay. Second claim: the relative time between TOA and each of the rTOA candidates is independent of flow of each voxel. For example, the relative time between TOA and TMD2 is the same for all voxels, independent of whether one voxel has a flow level of 20ml/100mg/min and another voxel has a flow level of 50ml/100mg/min. Of course, the second claim is true only if ET's assumption is met (e.g. TMD2 < τ).
The first claim is obvious as the relative time (from TOA) to reach a maximum parameter (e.g. MDX) of C(t) does not depend on the value of TOA. No knowledge of TOA is required. As a corollary, any time course parameters, e.g. C(t), D2(t), etc. obtained at rTOA are independent of the bolus transit delay even though they are sensitive to the accuracy of rTOA estimation. With regard to the second claim that the relative time distance between TOA and each of the rTOA candidates is independent of flow of each voxel, a lengthier explanation, given below in Appendix B, is offered. Appendix B shows how an understanding of the way rCBF is calculated by ET leads to an explanation that TMDX relative to TOA is independent of flow.
Based on those two claims, TMD2, TMD1.5, TMD1, TMD0.5 are considered appropriate rTOA candidates which can be used as proper reference points by ET for rCBF calculation.
Out of all the rTOA candidates, TMD1 is of additional interest because it matches the well known maximum gradient-based method used in computed Tomography (CT) and MR perfusion analysis. Since TMD1 is located at the superior contrast to noise location of the tissue time course one might question the need to study TDM2 and TMD-FR at all. The basic concern is about the possibility of TMD1 > τ and thereby not meeting ET's basic assumption. It is important to recognize that the time taken to reach the maximum gradient is dependent on the shape of the external driving function AIF which depends in turn on bolus volume, 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, resulting in drop of estimated perfusion values (Miles and Griffiths, 2003). On the other hand, TMD2 is much more likely to lie within the time window of τ because it is so close to TOA. But it may not be far enough from the influence of the low CNR of the TOA neighborhood. That is why we also explore the option of TMD-1.5 because the time to reach the maximum fractional derivative between the value of 1 and 2 should define an rTOA time point between TMD1 and TMD2.
In a flow ratio calculation, it should be mentioned that whether TMDX exceeds τ is determined by the τ of the higher flow level. For example, in a gray-white matter flow ratio, τ refers to τ of gray matter.
The presentation so far does not imply that there are always voxels where TMD1 is longer than τ. But until it can be thoroughly demonstrated that TMD1 is always shorter than τ under all kinds of imaging settings and a wide range of clinical conditions, it is useful to consider TMD2 and TMD1.5 as available backups.
Our manuscript is built on the premise that TMD2, TMD1.5, TMD1 and TMD0.5 are all currently useful and valid rTOA's to be considered for ET. While in principle we have almost an unlimited number of ways (Kwong et al., 2011) of evaluating rCBF in ET (e.g. one can choose C(t) at TMD1 or at TMD1.5), we will focus in this manuscript only on the study of the maximum values, namely the values of MD2, MD1.5, MD1 and MD0.5 for rCBF evaluation. For noiseless data, the same calculated rCBF is expected to be obtainable by ET no matter which type of approach we choose to use. For noisy real life experimental data, there are tradeoffs between using signal intensity of C(t) and the various values of MDX. There are practical advantages of using those maximum signals MD2, MD1.5, MD1 and MD0.5 at the respective rTOA's. The reason is that in noisy data, the inaccuracy in the estimation of the TMDX can derail the estimation of rCBF by signal intensity. The use of MDX at TMDX is considered more robust because the value of MDX remains relatively stable ranging from TMDX-1 to TMDX+1 while C(t) could vary significantly in the same time range. Error in estimation of TMDX could potentially compromise the accuracy of calculated rCBF more if C(t) instead of MDX is used for ET.
Given that the use of MD1 is a well established perfusion protocol, The result of MD1 is a good reference for the results of MD2 and MD-FR. The work flow of this manuscript consists of first studying a set of noiseless simulation perfusion data to verify the theoretical properties of TMD2, MD2, TMD1.5, MD1.5, TMD1, MD1, TMD0.5 and MD0.5. The simulated data is a subset of the simulated perfusion data created using the algorithm of Wu et al. (Wu et al., 2003a). We then investigate experimental monkey 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. Short TR image acquisition, which is the best way to collect many data points between TOA and τ to test and verify the ET model, is not yet an approved technology for clinical imaging,. That is the reason why we acquired monkey instead of human perfusion imaging data with short TR MRI. Proposals to manage the shortcomings of reduced brain coverage of short TR acquisition have been presented previously (Kwong et al., 2011) and won't be investigated here.
A number of factors needs to be taken into account in comparing the perfusion results of simulation data and real world experimental data. Calculated perfusion results of simulation data using various rTOA's would be expected to match the pre-given perfusion values. The inevitable noise of experimental data (Kwong et al., 2011) would introduce uncertainties into the results of rTOA's and MDX, degrading the accuracy of the estimated perfusion results. In addition, the assumption of a single transit time in the residue function of the simulation data is bound to be just an approximation for experimental data because partial volume of gray, white matter and vascular tissue in a single voxel can give rise to a mixed compartment of tissue transit times. Despite that approximation, previous work showed that rCBF calculated by ET for experimental data yielded reasonable results consistent with what were reported in the literature (Leenders et al., 1990; Ostergaard et al., 1996a). As long as consistent perfusion results can be obtained, the new rTOA and perfusion information retrieved from simulation data by ET can serve as an approximate and useful model for the study of real world experimental data.
In this manuscript, we will investigate whether flow can be estimated by maximum values of higher order fractional derivatives instead of using the maximum derivative which runs a larger risk of violating ET's assumption. We will also study whether experimental rCBF calculated using MDX is consistent with what was previously reported in the literature. While we verify such consistency, any discrepancies observed between the ideal ET model and experimental data will serve to suggest ways to improve the modeling of experimental data and improve the calculated rCBF in future research.
Tissue signal curves were simulated TR=300ms with an AIF modeled by a gamma variate function (Calamante et al., 2000) using techniques previously described by Wu et al (Wu et al., 2003b). We simulated noiseless bolus concentration time courses which reflect multiple flow conditions. We picked eight levels of flow value ranging from 10ml to 80ml/100g/min in 10 ml/100g/min increments. For convenience, the eight flow levels were named as f10 (10ml/100mg/min), f20, f30, f40, f50, f60, f70 and f80 respectively. The same blood volume was given for each flow level so there were eight τ values ranging from 2.98 sec for f80 to 23.8 sec for f10. To offer an example of the properties of ET under the ideal condition, we picked a rectangle model (Buxton, 2002) of the residue function R(t) 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). A rectangle function is defined as rect(t). With rect(t), no contrast agent leaves the tissue before τ. An R(t) made of rect(t) which assumes a single transit time can be used to examine and clarify 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).
Experiments were performed in a 3 Tesla Siemens Trio scanner (Siemens Medical System, Erlangen, Germany) using a customized surface coil. For routine Gd-DTPA bolus imaging, we applied gradient echo (GE) echo planar imaging (EPI) with TR=300ms and TE=31ms to acquire six slices with 600 time points (3 min) and the Gd-DTPA injection taking place at around the 250th image (75 s). Spatial resolution of GE EPI was 3mm × 3mm × 3mm.
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 the Siemens 3T system. The monkey was anesthetized with ketamine/xylazin, held in place by a MR compatible head frame. The monkey was imaged by GE-EPI while being injected once with Gd-DTPA (0.2mmol/kg) flushed with 12ml of saline with an injection rate of 2ml/sec by a power injector. The injection site was a leg vein.
The baseline SNR of the monkey's bolus signal, ranging from 80-100 for different gray matter regions, was higher than SNR=20 of the simulation data modeled after clinical perfusion imaging at 1.5T MRI.
To evaluate flow ratio with noiseless simulation data, we pick f10 (10ml/100g/min) as a reference flow level and take flow ratio using other flow levels. We calculated TMD2, MD2, TMD1.5, MD1.5, TMD1, MD1, TMD0.5, MD0.5 using routines from Matlab (MathWorks Inc., MA). The fractional derivative Matlab code, written originally by Farshad Merrikh bayat and modified by Mykhaylo Kholodov, utilizes the Grünwald-Letnikov definition (Oldham and Spanier, 1974). The first element of the data input is assumed to have the value zero. Caution is required using Matlab which has an intrinsic problem of delivering erroneous results of gamma function dependent coefficients used in the Grünwald-Letnikov definition if the data sample size is larger than approximately 200 data points. This Matlab limitation can be overcome by either supplying precomputed gamma function dependent coefficients using programs outside of Matlab or by simplifying gamma function dependent coefficients into elementary functions (Mykhaylo Kholodov, personal communication) using a similar approach described in (Das, 2008).
The flow ratio result, evaluated by the ratio between MDX of a particular flow level and that of the reference flow level (f10), is compared with the true flow ratio of simulation data to assess the usefulness of MDX for flow calculation. Flow ratio obtained by MDX is also calibrated against flow ratio results obtained by the early time points preceding τ (Kwong et al., 2011). Based on the calibration and comparison results, we can assess whether the use of MDX yields the true flow ratio results as well as whether TMDX meets ET's basic assumption. The X of MDX and TMDX ranges over 2, 1.5, 1, and 0.5.
The same analysis approach is also applied to noisy experimental monkey's data. For experimental data, we also calculated TMDX and MDX, with X going from the derivative order of 2 to 0.5, on a voxel-by-voxel basis. MD2, MD1, MD1.5 and MD0.5 were then used to make rCBF maps. For gray-white flow ratio calculation, we picked gray matter regions of interest (ROI, 3×3 voxels) from the left caudate and white matter ROI from deep white tissue on the same left hemisphere.
The estimation of MDX and TMDX faces the challenge of noisy experimental data as well as the noise introduced by the process of differentiation. The simplest way of obtaining the slope value by the subtraction of two noisy neighborhood points could introduce significant fluctuation and serious error to the slope result. That is why the time course data generally need to be smoothed or fitted with a model function. For the monkey's data, we smoothed the signal intensity (ΔR2*) of our bolus concentration time courses with a running average of 8 data points and then took D2(t), D1.5(t), D1(t) and D0.5(t) of the smoothed data set.
Experimental gray-white matter flow ratio time courses were made in comparison with the flow ratio time courses of simulation data. In a previous report (Kwong et al., 2011), gray-white matter flow ratio time courses were made using TOA as the reference time point. In this manuscript, a particular rTOA was selected as a reference time point to be used by all the gray-white matter flow ratio time courses to correct their respective time delay. Our choice of rTOA was TMD1.5. TMD1.5 was chosen as rTOA over TMD1 because there is a larger chance for TMD1 to exceed τ. TMD2 was not picked as rTOA because CNR is lower at TMD2. TMD1.5 is assumed to be a good compromise. Once rTOA was set, time-delay corrected gray/white matter flow ratio time courses can be generated using D2(t), D1.5(t), D1(t) and D0.5(t). Then we can examine and compare flow ratio obtained by MDX and by information from the early time points of the flow ratio time courses.
As expected, TMD2 < TMD1.5 < TMD1 <TMD0.5 in simulation results. Fig. 2 demonstrates this ordering given flow level f70 (70ml/100mg/min). Fig. 2 shows that at τ of f70, TMD2, TMD1.5 are smaller than τ while TMD1and TMD0.5 are larger than τ. Hence, TMD2 and TMD1.5 meet ET's basic assumption while TMD1 and TM0.5 fail it. Basically, Fig 2 shows what TMDX of C(t) is like under the condition of ‘same f, different D’ (same flow level, different derivative orders of DX(t)). The maximum values of TMD2 and TMD1.5 (MD2 and MD1.5) can then be taken to represent rCBF accurately.
Given the same known TOA for all the simulation time courses, our results demonstrate that relative to TOA, each one of TMDX is independent of flow levels as long as TMDX < τ. Table 1 below summarizes the relationship between TMDX, flow levels and τ for the simulation data.
Fig. 3 a,b,c,d present sample results from Table 1 using a graphical presentation of DX(t) (X=2, 1.5, 1, 0.5) of three flow level examples of f10, f30 and f70. In contrast with Fig. 2, Fig. 3 demonstrates what TMDX is like with ‘same D, different f’ (same derivative order of DX(t), different flow levels). It shows graphically that TMDX is independent of flow levels as long as TMDX is shorter than τ of the particular flow level. If TMDX > τ, the values of TMDX would be different for different flow levels. In Fig. 3a, the maximum values of D2(t), namely TMD2, are the same for f10, f30 and f70. TMD2 < τ's of all three flow levels. In Fig. 3b, the maximum values of D1.5(t), namely TMD1.5, are again the same for f10, f30 and f70. TMD1.5 < τ's of all three flow levels. In Fig. 3c, the maximum values of D1(t), namely TMD1, are the same forf10 and f30. TMD1 < τ's of f10 and f30. But TMD1 is of a different value for f70. TMD1 > τ of f70. In Fig. 3d, the maximum values of D0.5(t), namely TMD0.5, are the same for f10 and f30. TMD0.5 < τ's of f10 and f30. But TMD0.5 is of a different value for f70. TMD0.5 > τ of f70.
The flow ratio result obtained by MDX is independent of X, the order of the derivative, and is equal to the true flow ratio. Flow ratio was obtained by dividing the MDX of any flow level by the MDX of reference flow level of f10. With X ranging from 2 to 0.5, the use of MD2 gives the correct flow ratio for flow levels ranging from f10 to f80. The use of MD1.5 also gives the correct flow ratio up to f70. The use of MD1 gives the correct flow ratio up to f50. The use of MD0.5 gives the correct flow ratio up to f30. In the case of MD1, the calculated flow ratio of f70/f10 and of f80/f10 is lower by 5% and 11% respectively than the true value. Similar error percentage could be obtained for MD0.5 at the high end of the flow ratio. It should be noted that such error results are obtainable due to designed characteristics of the simulated τ value increasing at a steady pace as the flow level decreases. Such results present a model of the effect of TMDX and τ on flow ratio but should not be considered literally to be applied to experimental data.
Fig. 4 provides an example (f70/f10) of the error of the flow ratio calculated by ratio of MDX of f70 and f10 when TDM1 and TDM-0.5 exceed τ of f70. In addition to showing the ratio of MDX, Fig. 4 also shows the full time course of flow ratio (f70/f10) starting from the known TOA of the simulation data. The flow ratio is 7 and its time course is evaluated four separate times using the separate values of D2, D1.5, D1 and D0.5. Fig. 4 shows that MD2 and MD1.5 give the correct flow ratio of 7 at their respective TMDX while MD1 and MD0.5 do not. The reason is simply that TMD2 and TMD1.5 are shorter than τ of f70 while TMD1 and TMD0.5 exceed τ. As expected, exactly the same flow ratio of 7 is obtained before τ for early time points of any order of derivatives. The flow ratio time courses calculated using D2(t), D1.5(t), D1(t) and D0.5(t) all start deviating from the true flow ratio value of 7 at τ. D0.5(t) deviates from the true flow ratio much more slowly beyond τ than the higher order derivatives. In general, the lower the order of the derivative, the more time beyond τ would be taken to reach the same percentage of calculated flow ratio error.
From two different MRI sessions, we obtained the gray-white matter flow ratio from brain tissue regions of interest of similar slice locations.
Between a caudate gray ROI and a white matter ROI (Fig. 5) from two MRI sessions separated by days, the gray-white matter flow ratio result was 3.1 ± 0.2 using MD2, 3.2 ± 0.3 using MD1.5, 3.2 ± 0.3 using MD1, and 3.1 ± 0.3 using MD0.5. Despite the inevitable slight variation in slice imaging positions from two separate days, a certain reproducibility could be inferred from the comparable ET results from these two MRI sessions. The flow ratio values (approximately 3) calculated with MDX were all reasonable in comparison with literature values (Leenders et al., 1990; Ostergaard et al., 1996a). The same gray-white matter flow ratio results were also consistent with those reported previously of 3.1±0.2 using the deconvolution method (Kwong et al., 2011) and of 3.2±0.3 using ET which employed TOA as time reference and evaluated flow with the slope values (but not MD1) of early time points (Kwong et al., 2011). Even MD0.5, with TMD0.5 giving the largest chance of exceeding τ of gray matter, delivered a mean value not far off from the reported flow ratio.
In addition to pointing out the gray and white matter ROI's used for flow ratio calculation, Fig. 5 also provides an example of a perfusion map given by a MD1.5 map showing good gray white contrast.
Fig. 6 offers a closer look of the relationship between flow ratio for experimental data. TMDX and MDX are provided (Fig. 6) by time courses made from DX(t) (X being 2, 1, 1.5 and 0.5) of experimental data from one imaging session, analogous to the simulation flow ratio time courses of Fig. 4. Unlike Fig. 4 where all the simulated time courses share the same known TOA, the experimental flow ratio time courses of Fig. 6 of different DX(t) are aligned using TMD1.5 as the reference rTOA to facilitate the comparison of TMDX. In Fig. 6, the flow ratio value reported by MD0.5 in Fig. 6 is slightly lower than that of MD1, MD1.5 and MD2, but the lower MD0.5 value is not statistically significant when results of MD0.5 from two different imaging sessions are averaged together (results reported in previous paragraph).
As expected from the example of simulation data of Fig. 4, D0.5(t) of Fig. 6 has the slowest rate of deviation from the early steady-state portion of the flow ratio time courses, in comparison with that of D1(t), D1.5(t) and D2(t), potentially limiting the damage of any flow ratio error generated by the possibility of TMD0.5 exceeding τ. Based on visual inspection, departure from the early steady state portion of the flow ratio time courses started showing up between TMD1.5 and TMD1 for D2(t) and D1.5(t), suggesting that an experimental “τ” may be close to that departure location. The experimental “τ” is presumably related to the mixture of transit times of gray matter ROI. But this observation of experimental “τ” is qualitative rather than quantitative due to experimental noise. For example, given that data error of DX grows with the order of the derivative, namely noise of D2(t) > noise of D1(t) and so on, it is difficult to ascertain whether the point of sharp departure of the deviation of D2(t) is more a consequence of noise or of the experimental “τ”. Also, even if the time location of experimental “τ” appears to be very close to or somewhat smaller than TMD1, MD1 and even MD0.5 are consistent with MD2 and MD1.5 within error of the two imaging sessions.
We have assumed global AIF. Global AIF means only one AIF is assumed for all the voxels in the brain. Local AIF means different regions may have their own AIF's. While one would expect that a local AIF model is more realistic than a global AIF model, capable of leading to more accurate rCBF estimation, such a conclusion has not been borne out in recent research (Christensen et al., 2009). The lack of clear results at this moment could be related to the basic assumption of local AIF itself or to the current tools used to analyze local AIF models. More research is anticipated in the future. 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 voxel or region of the brain remain a challenging technical issue. Even when a particular local AIF model was selected for perfusion analysis, a report (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. Hence, we focused on global AIF in this manuscript and leave some of the difficult questions of local AIF to a companion manuscript (Kwong and Chesler, 2011).
One of the difficulties of getting adequate information on the relationship between TMD1 and τ in the traditional maximum gradient based method for perfusion analysis is the relatively long TR acquisition time of 1.5 sec (Tombach et al., 2003) in most routine clinical MR perfusion imaging, a compromise between brain imaging coverage and sampling rate. A sampling rate of 1.5 sec will not be able to give an accurate estimate of τ shorter than 3 sec by Nyquist consideration alone. In other words, TR=1.5 sec favors the report of τ longer than 3 seconds. In experiments, the long TR of 1.5 sec is not well suited for slope analysis. With only two to three data points to cover the range of experimental τ values, the chance of obtaining an accurate slope value is diminished considering the noisy fluctuation of data and the propagation of error in the difference operation. Accurate information of the relationship between TMD1 and τ would benefit in the future from a more comprehensive data base of τ obtained by short TR imaging acquisition.
The inadequate sampling rate of 1.5 sec is a general consideration independent of any perfusion analysis technique. In perfusion literature reports based on the deconvolution method, it had always been observed that too long a data sampling interval led to an underestimate of flow (Knutsson et al., 2010), which is an indirect way of saying that an inaccurate overestimation of τ would result. We chose a short TR=300ms for our noiseless simulation and experimental monkey studies. TR=300ms is selected to take into account the shortest τ imaginable to cleanly demonstrate the principle of rTOA. But TR=300ms is not intended to be a requirement for all imaging studies. If the τ of the tissue of interest happens to be 4 sec or higher, than a TR=2 sec would be adequate provided that other external conditions such as the bolus volume and the bolus injection rate are favorable.
How much we need to worry about TMD1 or TMD0.5 being longer than τ can be resolved only by larger sample size future studies which can provide a good evaluation of rCBF in the context of τ vs. TMD1 or TMD0.5 under many different experimental and clinical conditions. Given the preliminary results in this manuscript, TMD1.5 appears to be a relatively safe compromise between CNR and ET's basic assumption. One possible way to increase the chance of TMD1/TMD0.5 meeting ET's basic assumption is to inject an extremely tight contrast agent bolus. While there is an upper limit to how fast a Gd-DTPA bolus can be injected, it is conceivable to choose some brand of contrast agent of especially higher concentration to reduce the bolus length and therefore TMD1/TMD0.5.
In Fig. 6 of experimental data, we selected rTOA to be TMD1.5 and obtained time-delay corrected flow ratio time courses from the early time points of D1(t), D1.5(t) and D2(t). The estimated values of the early steady-state portion of the flow ratio time courses could vary a little had we chosen TMD2 or TMD0 instead as rTOA because in the ratio, the time course of the gray matter ROI could be shifted relative to the time course of the white matter ROI due to inaccuracy of TMDX in experimental data. On the other hand, the values of MD1, MD1,5 and MD2 should be independent of the rTOA chosen. That is why using the maximum values of the derivatives should be more robust for rCBF estimation than using parameters whose accuracy is more sensitive to the accurate estimation of TMDX.
The objective of this manuscript is to study the feasibility and applicability of various MDX for rCBF estimation. The experimental flow ratio results obtained from MD1 are close to that of MD2 and MD1.5, suggesting that either TMD1 is close to an experimental “τ” value or, if flow ratio errors do exist due to TMDX > τ, difference between MD1 and MD1.5 may be masked by experimental noise. The mean flow ratio obtained from MD0.5 appears to be smaller (though not statistically significant after averaging two imaging sessions) than that of MD1, MD1.5 and MD2, as what is expected, but such a smaller value may also be maskable by noise. With our small sample size, it is possible for MDX with X ranging from 2 to 0.5 to all yield “reasonable” rCBF results difficult to distinguish from each other. Which order of the derivative is ultimately best for rCBF estimation can only be answered by a large and sufficient data base that will calibrate τ, TMDX and flow ratio adequately in the future. But our manuscript opens the door to using fractional derivatives for rCBF estimation to meet unique demands of different imaging occasions and we are no longer limited to a single maximum gradient based method. In our use of the MDX approach, the consistency of gray-white matter flow ratio of experimental data with the flow ratio reported in the literature (Leenders et al., 1990; Ostergaard et al., 1996a) lends support to the research direction of utilizing rTOA and fractional derivative.
The noise introduced by the process of differentiation is an ongoing concern shared by our rTOA approach and the traditional maximum gradient based method for perfusion analysis. The problem of noise amplification of differentiation is presented in Appendix B. The difficulty of searching for the derivative in noisy data is a long standing challenge in engineering and many fields (Culum, 1971; Hanke and Scherzer, 2001; Melanson et al., 2010), and differentiation of noisy data remains an on-going topic for research (Chartrand, 2011). Noisy data produce derivatives with slope alternating direction dramatically at every sample. Normal approaches to the problem include filtering and curve fitting by least squares but none is considered foolproof. With high SNR monkey images, we had shown at the Methods and Results sections above that simple filtering algorithms like running average were sufficient to obtain reasonably equivalent flow ratios calculated by MD1, MD1.5 and MD2. Simple filtering in this case is more convenient than curve fitting by least squares because there is one less concern over the choice of the curve fitting model.
Is it feasible to obtain meaningful MDX values at clinical SNR of around 20 at 1.5T MRI? The simulation data (Wu et al., 2003a) we adopted for this manuscript has a built-in SNR option of 20 with a Rician distribution. For such low SNR data, curve fitting is the preferred choice for clinical perfusion imaging data in the literature (Kimura and Kusahara, 2009; Koenig et al., 1998; Miles, 1991). The likely reason is the current lack of automatic criterion to determine the required amount of filtering to increase the SNR of derivatives. For illustration purpose, we also chose to make a preliminary assessment of the noisy simulation data by testing the gamma variate curve fitting procedure routinely used in the maximum gradient perfusion literature (Kimura and Kusahara, 2009; Koenig et al., 1998; Miles, 1991).
It should be noted that the validity of the gamma variate function fit for experimental bolus time courses had been disputed (Chen et al., 2005) and the gamma variate fit results are known to be quite sensitive to the choice of initial values. Nonetheless, for noisy simulation data we know that the AIF is a gamma variate function. By using the gamma variate fit, we can isolate the amount of flow estimate error attributed to low SNR and to the use of the maximum values of the derivatives (MDX) from the uncertain errors of having thoroughly inappropriate curve fitting models for the noisy data.
We actually applied the gamma variate function fit to the derivative of the simulated bolus time courses. This is a better choice than applying the fit to the bolus data itself because the derivative is theoretically close to the shape of AIF (a consequence of the ET/microsphere model; Eq (B3), Appendix B).
The gamma variate function for fitting a time course D(t) is
where t0, K, α, β are the fitting parameters and D^ (t) indicates the fitted result.
Once the fitted curve was obtained for the derivative of the noisy simulation time course, the fractional as well as the second derivative values were then derived from the fitted curve. The curve fitting results presented below demonstrate what flow estimate values are feasible at clinical SNR but are not considered to be complete or optimal.
Table D-I and Table D-II below give sample results of the ratio of the fitted MDX and true MDX, a proxy for the ratio of fitted flow and true flow. We started out with 50ml/100g/min of the simulation data because we know that the ET basic assumption was not met for TMD1 at 60ml/100g/min and above for the simulation data.
Table D-I is shown below. The results show that the standard deviation σ of the ratio between the fitted value of MDX and the true value of MDX gets larger as the flow gets smaller.
At Table D-II below, we also obtained ROI information by taking averaged simulation time courses from ROI's of 3×3 pixels, just like what was commonly done in clinical analysis of gray/white matter ROI. The results show a reduced σ as well as a smaller discrepancy between estimated and true flow for ROI, as expected for the increased SNR.
Inaccuracy in perfusion estimation occurs with any analysis method. In the evaluation of perfusion results using various deconvolution approaches on simulated data of SNR=20 (Wu et al., 2003b), ratio of estimated flow and true flow returned a Mean±σ ranging from 0.95±0.09 to 1.51±0.21 for flow levels from 70ml/100g/min to 10ml/100g/min. The results were obtained by the use of a block-circulant matrix for deconvolution (oSVD) and a box-car (rectangle) R(t). As there is no standard on how much inaccuracy in perfusion results is acceptable, the tolerance for perfusion inaccuracy by any methods depends on each unique clinical imaging requirement.
Discussion results of Table D-I and D-II are used mainly for illustrative purpose as the sensitivity of fitted outcome to initial values remains a focus of future investigations. Optimization constitutes a future goal. The results presented at Table D-I and D-II provide an idea on what flow estimate values are feasible at clinical SNR if an appropriate curve fitting model can be found for the data. Interesting results independent of initial values are patterns of how estimated results are related to high and low flow and/or to the SNR of the baseline signal. We expect that the accuracy in perfusion values estimated by ET would continue to improve in anticipation of the development of better filtering and curve fitting approaches to obtain more accurate derivative values in future research.
While TOA is critical for the evaluation of rCBF by ET, the low CNR of TOA encourages the search for other rTOA's such as the time to reach the maximum values of conventional derivatives or fractional derivatives. We have identified successful rTOA examples including TMD2, TMD1.5, TMD1, TMD0.5. While bolus signal intensity C(t) as well as a variety of orders of derivatives can all be utilized at TMDX to estimate rCBF, it is convenient and robust to use MD2, MD1.5, MD1 and MD0.5. The main message of this manuscript is the investigation of higher order fractional derivatives which enables us to avoid if necessary estimating flow with the maximum derivative which runs a higher risk of violating ET's assumption. For experimental data, the signals of MD1 and MD0.5 have better CNR than MD1.5 and MD2 but TMD1 and TMD0.5 run a stronger risk of longer than τ. Overall, MD1.5 is a reasonable compromise from our study results. The flexibility of using fractional derivatives could provide a solution for the trade-off between CNR and ET requirement. A comprehensive data base obtained by short TR imaging acquisition for many different kinds of diseases and disorders is needed in the future to provide meaningful information on the relationship between TMDX, τ and the acceptable amount of error for flow ratio calculation.
Good contrast to noise reference time points are proposed to correct bolus time delay in ET. > The time to reach the maximum fractional derivative of order 1.5 is a good reference time point. > In rCBF estimation, the maximum values of fractional derivatives can be used to replace the maximum derivative which runs a higher risk of violating ET's basic assumption. > Good rCBF results by ET are demonstrated with simulation and experimental data.
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). We would like to acknowledge Dr. Tuan Cao-Huu for his computation advice and his years of generous support of our research effort.
Given location x, we expect that the error of D1(x) < error of D-FR(x) < error of D2(x). The first derivative can be considered to be a difference between two consecutive sampled data points. For random noise, analytical solutions for errors of the first and second derivative are well known and can be obtained analytically using the propagation of error of arithmetic subtraction. The question is how to obtain a solution for errors of fractional derivatives.
In propagation of error, given S being a difference of consecutive and uncorrelated data points A and B of a data set f (x) which could be the baseline portion of a bolus time course
where , and represent variance.
To examine the errors of the derivatives, the results can be better visualized if we present different orders of the derivative of the data set f (x) in the following expressions:
For the (first order) derivative,
For the 2nd derivative,
For any integer (n) order of the derivative,
The binomial coefficient is what determines the results of the propagation of error. For example, if we assume the standard deviation σ of the f (x) to be 1, we obtain the variance for the (first order) derivative,
The variance of the 2nd derivative,
Now we know how to evaluate the theoretical error σ for D1(x) and D2(x), the next question is how to derive the error of a fractional derivative D-FR(x).
For fractional derivative, the term is generalized by in the Grünwald-Letnikov-derivative implementation (mentioned in Methods) where q could be a fractional number as well as an integer. In this generalization, the factorial terms in the binomial formula are replaced by the gamma function, namely, a continuous curve of fractional factorial values.
The values of the fractional derivatives are well characterized by the sum of the first three terms of the fractional derivative expression, with the contribution diminishing rapidly from subsequent terms. As an illustration, the coefficients of first seven terms of the expression for D1.5(x) are 1, -1.5, 0.375, 0.0625, 0.0234, 0.0117, 0.0068. It is therefore sufficient to approximate the values of different orders of the derivative by just the first three terms of the expression. Given the presence of fractional factorials, the absolute amplitude of the coefficients of the first three terms of the expression for a fractional derivative (between the order of 1 and 2) would be larger than that of the first derivative but smaller than that of the second derivative. For example, the coefficients are (1, -1.5, 0.375) for D.1.5(x) whereas they are (1, -2, 1) for D2(x) and (1, -1) for D1(x). Given that variance is the sum of squares of the (generalized) binomial coefficients, it can be derived that Var(D1(x)) < Var(D-FR(x)) < Var(D2(x)).
In Table A-I below, we calculated the variance and standard deviation of different orders of the derivative ranging from D1 to D2 by summing first the squares of the coefficients of the first three terms of the fractional derivative expression, and then summing those of the first 100 terms for comparison.
The theoretical curve of the variance or standard deviation of a continuous range of orders of fractional derivative can easily be graphed as a continuously rising curve. But it will not be done here.
Simulation data could be generated for comparison. We generated 1000 randomized data points f (x) of normal distribution from Matlab using an initial standard deviation of 1. The variance and standard deviation of the derivatives with orders ranging from D1 to D2 are shown at the Table A-II below.
Noise amplification for different orders of derivative becomes prominent if differentiation is done by taking the difference of neighborhood data points (the finite difference method). Filtering or curve fitting the data helps reduce the noise problem. Filtering is designed to increase signal to noise. In curve fitting, one often sidesteps the dependence of noise on the orders of the derivative. Instead, noise shows up in the quality of fit. For example, if data can be fitted by a simple polynomial, different orders of the derivative can simply be derived analytically as other polynomials. The real world challenge is finding the proper curve fitting model for the data. Many different filtering algorithms (e.g. moving average, Savitzky-Golay filter (Savitzsky and Golay, 1964), total variation regularization (Chartrand, 2011)) are available. Many curve fitting models (e.g. gamma variate function) also need to be tested. The necessary evaluation and comparison of many of these filtering and curve fitting techniques require future research.
Perfusion information can be obtained from the contrast agent concentration time curve in the general expression given below:
where C(t) is the tissue contrast agent bolus concentration signal time curve, f 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). In ET, the residue function R(t) = 1 when ET's basic assumption holds, leading to Eq (2) below which confirms that perfusion result evaluated by ET is independent of AIF (global AIF is assumed).
From Eq (B1), the convolution of AIF with R(t)=1 can be derived analytically as
Since ET is based on the microsphere model, the evaluation of rCBF by ET is however not limited to the use of C(t). One can use as well the integral of C(t) or derivatives of arbitrary order of C(t) (Klocke et al., 2001; Kwong et al., 2011; Lee et al., 2004). For our purpose here, we will focus only on the use of the derivatives. Given R(t) =1 based on ET's basic assumption, the first derivative of the concentration time course can equally be used by ET to calculate rCBF because
By the same consideration, the second derivative, the third derivative, and in fact almost any number of arbitrary ways of manipulating the early contrast agent bolus concentration time course can equally be used by ET to calculate and obtain the same rCBF values. This is indeed true with noiseless simulation data (Kwong et al., 2011). With real life experimental data, there will be signal to noise ratio (SNR) tradeoff with each different manipulation of C(t). But we are just laying down the principle here without going into the details of how each unique manipulation may have its own niche to play. While the study of (Kwong et al., 2011) is more comprehensive than the long established maximum gradient based method (Kimura and Kusahara, 2009; Koenig et al., 1998; Miles and Griffiths, 2003) which focuses on just the single maximum value of , the maximum gradient based method can indeed be recast in the framework of a version of ET analysis.
Having established that the time courses of derivatives of all orders can be used by ET for rCBF calculation, we can now explain the claim that the time distance between TMDX and TOA is independent of flow. It can be derived directly from Eq (B4) which shows that that if R(t) = 1 and AIF is global, TMD1 relative to TOA is independent of f. The idea is straight forward. The time to reach the maximum value of is the same as the time to reach the maximum value of for any f. The same argument goes for TMD2 and TMD-FR.
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.