El suiguiente texto ha sido generado con un teclado azerty, por ende se han omitido tanto la tilde como la "ñ". Presento excusas por los inconvenientes que dicho maltrato a la lengua española, haya podido generar.


Este es el trabajo que desarrolle durante el ultimo ano. La primera razon por la cual la publico en esta página es porque quería estar seguro de que no se perdiera la información encaso de una falla con mi PC. En segundo lugar, quería compartir lo que he hecho con la gente, especialmente con mis amigos, para que sepan un poco en que he andado metido. El titulo puede llegar a ser un poco pesado o poco claro, sin embargo, espero que la introducción les ayude en algo. Si no, me cuentan y con mucho gusto les cuento de que se trata en palabras sencillas.

Esta escrito en ingles pues aspiro homologarlo como mi trabajo de grado en la Javeriana y supuse que no seria tan facil si lo hubiera hecho en frances. Les recomiendo leer la introduccion unicamente pues no inclui las ecuaciones y al parecer se pone un poco denso...

Hasta pronto.

Atte:

JPS


Deconvolution Robustness for Local Perfusion
Parameters Estimation in Magnetic Resonance Imaging (MRI)
Juan Pablo SANTAMARIA
June 2004

INTRODUCTION

Cerebral and cardiovascular diseases caused by the obstruction of blood vessels, are one of the major causes of death in our modern society. The medical community has determined this class of disorders as an important public health threat.

Magnetic resonance imaging (MRI) technology has been widely used in clinical practice as a non-invasive technique to visualize the inside of the human body and to detect health-threatening situations in patients. In many cases, it is required to image the patient flowing blood, with the purpose to detect anomalies in the circulatory system. Perfusion MR imaging is a suitable method for this task, in which a contrast agent is introduced into the blood stream and then, it is digitally imaged as it travels across the tissue. MR perfusion imaging is especially a good choice for diagnosis and follow-up studies of several injuries in the arterial system.

Consequently, the study field for new concepts and applications within MRI medicine is today a vast and dynamic area, as it was proved last year by Paul C. Lauterbur and Sir Peter Mansfield , 2003 physiology/medicine Nobel prize winners, for their discoveries concerning magnetic resonance imaging.



Figure 1. A typical RM brain image sequence

When thrombosis is detected, it is also often desired to evaluate the regional cerebral blood flow (rCBF). This perfusion parameter quantifies the blood flow at the capillary level where the delivery of nutrients takes place. Its value is an indicator of the tissue viability, which is useful in taking appropriate treatment decisions. However, in clinical practice there is still a need for robust diagnostic tools for determining perfusion parameters in environments affected by random electrical noise.

The purpose of this study was to determine the robustness of two deconvolution techniques in rCBF estimation: The singular value decomposition (SVD) which is commonly used in the literature and the ARMA model that was developed at the CREATIS laboratory.

Each algorithm was tested on 150 data sets in two distinct cerebral blood volumes, and we defined error measurements in order to compare the performance of the different methods. Computer simulations here were not intended to substitute practical data but were only applied as a first approach before analyzing the clinical samples.

2. THEORETICAL BACKGROUND

2.1. Biological system modeling

2.1.1. Dynamic systems

The physiological behavior of some organs can be modeled using the dynamic systems theory. The heart, veins, arteries and capillaries may be treated as a biological system, mathematically representable, where variables and parameters interact regularly. This representation provides a powerful tool to estimate variables when some parameters are unknown. Another advantage of the dynamic systems theory is that we can handle various -apparently different- physical parameters (e.g., temperature, flow, concentration and voltage) without modifying the equation’s essence. In recent years, the study field for new concepts such as the biochemical circuit theory analysis is today a fundamental field of interest for many laboratories.

The physical elements treated in this medical context, are the perfusion parameters, i.e., regional cerebral blood flow rCBF, regional cerebral blood volume rCBV, mean transit time MTT. However, the same parameters and relationships can be deduced and applied in a different area such as the myocardial circulatory system. In the published literature, the notation is quite similar except for that the letter M (myocardial) replaces the letter C (cerebral), in the perfusion parameter abbreviations.

The human circulatory system is divided into two main parts with the heart acting as a double pump. This organ pumps of a special kind of fluid, which is the blood. The across variable in this case, the blood flow, leaves the left side of the pump (heart) and travels through arteries which gradually divide into capillaries. In the capillaries, the nutrients are released to the body cells. The blood then travels in veins back to the right side of the pump, and the whole process begins again.

One of the most common ways in which this kind of process is analytically described is the black-box model. The biodynamic system has, therefore, the following components:

 Parameters: the rCBF and rCBV.
 Dependent variables: the arterial, tissular and venous concentration, denoted Ca, Ct and Cv respectively.
 Independent variable: time.

Figure 2.
The blood circulation as a dynamic system

2.1.2. The Stewart-Hamilton physiological model

The Stewart-Hamilton model applied to MRI acquisition techniques can estimate perfusion parameters such as the regional blood flow rBF, mean transit time MTT or regional blood volume rBV, describing the tissue condition.

The concentration of tracer within a tissular volume, at a given time t, during the passage of a bolus injection of a contrast agent, is given by:



Where denotes the tissue density (i.e., the tissue mass per unit volume), Cartery is the arterial concentration and R(t) is the residue function and represents the tissular retention of the contrast agent. If this biodynamic system has a single fluid storage element hence, the residue function can be treated as the response of a first-order circuit. The general solution for this differential equation is of the form:



In this particular case we have: is the mean time taken by the blood to pass trough the system, from the arterial entrance to the vascular exit and k=1.

Let

Therefore



Assuming the tissular density being close to the water density, i.e., ,
then:


Once the regional blood flow is estimated, the unknown perfusion parameters, can be computed from the central volume theorem (Stewart, 1984; Meier and Zierler, 1954):



where rCBV and MTT are the regional cerebral blood and mean transit time, respectively and



However, the concentration is measured from the MR digital images, therefore, the continuous time form of the equations is not considered. The Stewart-Hamilton model becomes a discrete convolution:



Moreover, the discrete convolution can be viewed as a matricial product:



Where N is the number of samples taken in the RM sequence.

This product can be compactly denoted as:



Where is unknown, then:


Hence, the deconvolution process becomes a linear algebra problem, specifically, a matrix inversion problem. However, this process, when working with noisy signals, could sometimes lead to determinants close to zero, the matrix A is said to be “ill-conditioned”, and consequently, unexpected results may be obtained.

Moreover, the concentration of the contrast agent is not directly measured in clinical practice; instead, MR images are used to determine signal intensity (SI) variations over time, which are related with the change in the spin-spin relaxation, also known as the transverse relaxation . The SI functions are obtained by selecting a pixel or a region of interest (ROI) in the MR gray-scale digital images.

Previous studies have shown the existence of a linear relation between the measured concentration C(t) and . On the other hand, the transverse relaxation and the signal intensity have an exponential dependence, therefore, we write:



Where TE is the sequence echo time, S(t) is the signal intensity over time and S0 is the baseline MR intensity.

Consequently, by measuring SItissue and SIaif we can estimate the regional blood flow using a deconvolution technique and evaluating the expression r(t) at zero.

2.2. Deconvolution methods

2.2.1.1. General Singular Value decomposition (SVD)

The singular value decomposition is a widely used technique to solve ill-conditioned problems with several applications, (e.g., in image compression, watermarking, image filtering). The general SVD statement can be expressed as:

Every real matrix A can be decomposed into a product of three matrices of the form:



Where, U and V are orthogonal matrices and with . Therefore, S is a diagonal matrix whose elements are the singular values of the original matrix.

Consequently, the inverse can be expressed as:



Where W is also a diagonal matrix of elements

Depending on the signal-to-noise ratio of the signal intensity function, a tolerance threshold PSVD is set. The values of wi corresponding to values where si is less than PSVD are set to zero. Typically the threshold is given as a percentage of the greater singular value . Generally the threshold value varies between 20 and 30% and the general principle is the higher the noise, the higher the PSVD.


2.2.1.2. Singular Value decomposition adaptive threshold (pSVD)

Liu et al. (2) showed that the PSVD selection had a significant influence in the shape of the residue function and an apparently inaccuracy in the rCBF estimation. An oscillation index O, was proposed by Østergaard et al. (1) to measure the distortion in the residue function as follows:


Where f is the scaled estimated residue function, fmax is the maximum amplitude of f, and L is the number of sample points.

This oscillation index is the discrete form of the convolution product between the residue function and the second derivative of a Gaussian distribution. The use of this digital filter, which is the same as calculating the sum of the absolute value of the second derivative of f(t), determines the change in the curvature over the time.
Figure 3.
The Gaussian distribution (a) and its second derivative (b)


2.2.2. ARMA

The estimation of the residue function can be treated as a spectral parametric analysis problem. It is known that when a linear time-invariant (LTI) system is excited by the Dirac's delta function, the output is the transfer function of the system.
Figure 4.
Discrete time system - ARMA model representation

The discrete time filtering theory states that the impulse response of a linear system h(n) can be modeled using one of the following frequency responses:
The “zeros model” or moving average (MA):



The “poles model” or autoregresive (AR):



The “zeros and poles” model (ARMA)





Which can be written as:

The aim of the ARMA technique modeling is to find the transfer function H(z), whose impulse response h(k) approximates the scaled tissue response r(k), such that the sum of the squared error, denoted e, is minimum:

Considering the system illustrated in figure 2, and using the ARMA model, the output y(n) can be computed as



This mathematical statement generates the following set of simultaneous equations:



Which can be compactly expressed, in matrix notation as

,

where Y represents the output matrix, A is the input-output matrix and is the ARMA parameter matrix (coefficients ak and bk).

Therefore, the general ARMA-model matricial representation can be written as



To solve the linear equations system , an estimation square error is defined as


where

The minimization of the error by the least mean square process gives the coefficients ai and bi, and consequently the searched solution.



where

To minimize this equation the SE must be differentiated with respect to and equated to zero.

Solving the equation we have


Where represents the unknown ARMA parameters.

3. METHODS

3.1. Monte–Carlo simulations
In order to assess the performance of the different deconvolution techniques, realistic MR signals were generated. The statistical data samples of the Monte Carlo simulation were the estimated blood flows, which were calculated from noisy concentration functions. These simulated signals were considered as a suitable approximation of what is typically measured in patients.

Range of simulated perfusion parameters
Representative rCBF values were selected from the scientific literature. The rCBF range varying from 5 to 35 ml/100g/min in 5 ml/100g/min increments was chosen for a volume rCBV = 2%, corresponding to healthy white matter. Afterwards, the rCBF was used to calculate the MTT characteristic values, applying the central volume theorem.

The same process was used for rCBV=4%, corresponding to healthy gray matter except for that the range varied from 10 to 70 ml/100g/min, in 10 ml/100g/min increments. These values are summarized in the following table:


MTT (s) 24 12 8 6 4,8 4 3,43
rCBF(ml/100g/min) RBV=2% 5 10 15 20 25 30 35
rCBF(ml/100g/min) RBV=4% 10 20 30 40 50 60 70

Table 1.
Simulated perfusion parameters. Mean transit time, regional blood flow and regional blood volume.

Noise simulation procedure
In order to simulate each of the rCBF in both blood volumes values (i.e., rCBV= 2% and rCBV= 4%), the following procedure was applied:

First, the simulated residue functions were generated as

for each of the 14 different rCBF values.

The arterial input function (AIF) was simulated as a gamma-variate function


where u(t) is the step function (due to causality in the modeled system).

The concentration tissue response (Ct) was then obtained by convolving the AIF with each residue function. SI functions were calculated assuming an exponential relationship between C(t) and the pixel gray-level:


where S0 is the baseline MR image intensity, k is a proportionality factor and TE is the echo time.

Additive white Gaussian noise was then used to impose the SI function samples with SNRtissue = 6 or 18 dB and SNRAIF = 15 dB, measured from our set of patient data. The low SNR in tissue corresponds to one pixel and higher SNR to a typical region of interest (ROI). It is very important to notice, that the noise can not be directly added to the concentration functions since the existence of a non linear relationship between C(t) and the simulated signal intensities. To simplify the notation, let the function be the noisy version of :

Once the noisy SI were generated, arterial and tissue concentrations were calculated using the inverse of the equation above:




From the moment the concentration signals were calculated, the scaled residue function was computed, applying the Stewart Hamilton model and one of the deconvolution methods:





For the ARMA deconvolution algorithm, the first poles-model order and a second order for the zeros-model was used. On the other hand, a fixed 30% threshold was first chosen for the SVD method.

Finally, the blood flow was calculated as



Every deconvolution method was tested in the 7 different blood flows values and the number of replicates for each one, denoted N, was of 150.

Deconvolution performance criteria

The mean of the simulated samples, for each of the estimated rCBF, was calculated and subsequently compared with the true values (see table 1).

To evaluate the performance of the deconvolution techniques, the percentage error (PE) and the standard deviation was calculated for different rCBF and blood volumes. The standard deviation for a discrete variable made up of N observations is the positive square root of the variance and is defined as:


where is the arithmetic mean of the set of samples and N=150.

The mean percentage error, denoted MPE, was defined as follows:


where N is the number of simulated flows and,



However, the deconvolution performance was not only assessed in terms of its error and standard deviation. Sensibility to SNR shifts or to sample rate was also considered in the analysis.

3.1.1. ARMA vs. SVD

The first part of the simulation was the comparison of the ARMA and the SVD deconvolution in the cerebral perfusion context. Similar studies have been done in the myocardial environment (Neyran et al. 2002), so different blood flow ranges are found in this project. Moreover, the sample rate was included as an additional issue in the analysis and it was tested in different noisy environments.

Two distinct tissular signal-to-noise ratios were used, in order to simulate an approximated realistic random noise when using single pixel or ROI selection. The different simulation environments are summarized in the following scheme:



Figure 5.
ARMA vs. SVD perfusion simulation strategy.

3.1.2. SVD vs. Adaptive threshold SVD

The adaptive threshold SVD deconvolution (aSVD) performance was tested and compared with the non-adaptive method for a sampling rate of 1 Hz in two different noisy environments, as shown in the following scheme:

The percentage threshold, was first set at 30 % and then a first oscillation index was calculated. The aSVD algorithm consisted in systematically change the PSVD in the range [20% , 40%] in 4 % increments in order to recalculate the respectively oscillation index of the residue functions. Afterwards, all the oscillation indexes were compared and the minimum was selected.

3.2. Perfusion MRI in stroke patients

Patients and image acquisition:
Nineteen patients (age range, 32-94) with symptoms of acute hemispheric stroke were retrospectively included in the present study, observing all respective medical and ethical regulations.

MRI studies were obtained within 6 hours from symptom onset using a 1.5-T Magnetom Vision whole body MR imager (Siemens, Erlangen, Germany). Perfusion-weighted MRI was performed with a T2*-weighted gradient-echo echo-planar imaging sequence, using the bolus passage of contrast agent (repetition time: 2000 ms; echo time: 60ms; 7 slices; slice thickness: 5 mm; interslice gap: 0.5 mm; field of view: 240 mm; matrix 128 x 128 voxels; 30 measurements obtained at intervals of 2s).

The contrast injection (15 ml of Gd-DTPA) was performed after the third scan, using a power injector at a rate of 5 ml/s via access through anantecubital vein; the bolus of contrast medium was followed by a 15-ml bolus of saline solution at the same injection rate.

Data post-processing
Perfusion data were transferred onto an independent work-station and analyzed using homemade programs written in Matlab language (Math Works, Natick, MA). To obtain the SI functions, tissular and arterial ROIs were manually selected and converted to concentration signals applying the logarithmic relationship. The cerebral blood flow estimation was then calculated using the ARMA and SVD deconvolution techniques.

4. RESULTS AND DISCUSSION

4.1. Monte-Carlo simulations

4.1.1. ARMA vs. SVD

Sampling Rate: 1 Hz (1 sample per second)

Non-noisy data
As expected, in the case of non-noisy data for a PSVD threshold of 30%, the SVD estimated flow was underestimated. The mean percentage error was –34%. In contrast, using the ARMA deconvolution, the true and identified cerebral blood flows were identical for both blood volume values hence the percentage error was 0%.

Figure 6.
CBF estimation without random noise using ARMA and the SVD deconvolution.
The selected threshold percentage was 30% for the SVD technique.

The estimation in a non-noisy environment was used to show the influence of the Psvd percentage threshold selection in blood flow estimation. The underestimation in the SVD method was due to the elimination of non-noisy information contained in the singular values matrix.

Those elements are not supposed to be removed without the presence of random noise. Obviously, if a small (close to zero) Psvd percentage threshold had been chosen, both results would have been the same. The variation effect of this parameter is shown in figure 5. Note that the appropriate Psvd for a high SNR tends to zero.

Noisy data in ROI simulation

For SNRtis = 18 dB and SNRaif = 15 dB the identified rCBF was underestimated in the ARMA simulation, with a mean percentage error of –6.9% and a mean standard deviation of 5.8. In SVD, there was also an underestimation, however, the absolute value of the mean percentage error was greater than the ARMA one. On the other hand, the standard deviation of the SVD was less than the ARMA one. This phenomenon was valid for both blood volumes as is illustrated in figure 6.

Noisy data in pixel simulation

When the tissular SNR was decreased to 6 dB without changing the arterial signal-to-noise ratio (SNRaif = 15 dB) the calculated data was, in this case, overestimated using the ARMA model, with a mean percentage error of +86%. On the contrary, the identified CBF with the SVD deconvolution technique remained underestimated with a mean percentage error of -34.5 %.
Figure 7.
Influence of Psvd threshold selection in flow estimation without random noise in SVD deconvolution method.




Figure 8.
Standard deviation comparison for ARMA and SVD deconvolution techniques, for R=1Hz, SNRtis= 18 dB and SNRaif = 15 dB.


Figure 9.
This figure illustrates the SNR shift sensitivity in flow estimation using both deconvolution methods for a CBV=2% and a sampling rate of 1Hz.

Sampling Rate: Rs= 0.5 Hz (1 sample every 2 seconds)

Non-noisy data

When the sample rate was halved, there was no significant change in the rCBF estimation. The performance without the presence of random noise for Rs= 0.5 Hz, was similar for both techniques. That means a zero mean percentage error (MPE) using ARMA and a variable MPE in the SVD deconvolution, depending on the selected threshold:



Noisy data in ROI simulation

Sensitivity to sample rate shift:

The estimated blood values that were first underestimated in the ROI SNR, for Rs=1 Hz in both deconvolution methods changed to overestimation, when the sample rate was halved to Rs=0.5 Hz.

Specifically, the effect in the ARMA performance for a tissular SNR=18dB and an arterial SNR=15dB, for this sample rate shift, was similar to the simulation for Rs=1 Hz when the SNR was decreased. That means that the ARMA model passed from underestimation to overestimation when the sample rate was halved. For Rs=0.5 Hz the mean percentage estimation error was of 70%. This variation represents a mean relative increment of 0.4 in the initially identified flow with a sample rate of 1Hz. The ARMA standard variation increased in average 50%.
Figure 10.
Comparison of SNR shift sensitivity in blood flow estimation for both deconvolution methods. CBV=2%, SNRtis= 18 dB and SNRaif = 15 dB.

The SVD deconvolution passed from underestimation to overestimation as well, when the sample rate was changed to 0.5 Hz. This change represented a mean increment of 0.5 into the estimated flows for Rs=1 Hz. The standard variation for the SVD method increased in average 60 % in both blood volumes after the sample rate shift.

Noisy data in pixel simulation
Sensitivity to SNR shift:
When the random noise was increased (SNRtis was changed from ROI to pixel level) the ARMA deconvolution passed from an MPE overestimation of 68% to 135%. Its mean standard deviation increased as well from 7.5 to 15.6 in other words there was a standard deviation increment of 108%.

On the other hand, SVD was less sensitive to the SNR shift, when the tissular signal-to-noise ratio was decreased to SNRtis=6dB as seen in the following figure. The MPE passed from a mean percentage error of 32% to a MPE=35% with this SNR change. The mean SVD standard deviation passed from 3.7 to 5.








Figure 11.
Comparison of tissular SNR shift sensitivity in blood flow estimation for both deconvolution methods. Rs = 0.5 Hz and SNRaif = 15 dB.



4.1.2. SVD and aSVD deconvolution

For PSVD = 30 %, the performance of the SVD method was not enhanced in terms of its mean percentage error, when the adaptive threshold deconvolution was used. The estimated flow was still underestimated in both blood volumes in the same proportions.

For pixel and ROI simulations, neither the MPE nor the standard deviation was improved for both cerebral blood volumes. The reason for this lack of improvement was because the selected PSVD was appropriate for the S/N ratio and therefore, the optimal oscillation index was already chosen in the SVD deconvolution.

Figure 12.
Example of oscillation index calculation of four different residue functions from the adaptive SVD algorithm.

Under these noisy conditions, the main difference found between the SVD and aSVD deconvolution was in terms of its execution time. It is important to keep in mind that the aSVD algorithm needed to test 5 different matrix thresholds, recalculating at each time by deconvolution, the estimated residue function. Therefore aSVD was much more expensive in terms of execution time than SVD.

4.2. Perfusion MRI in stroke patients


Figure 13.
Arterial and tissular concentrations calculated from SI functions for patient 6.



Due to the lack of a reference method for the assessment of perfusion, it was not possible to determine if there was under or overestimation for each case. However, this part of the experiment, was useful to determine the orders of magnitude of the rCBF in clinical practice and to assess the relationship of the ARMA and SVD deconvolution.

A set of perfusion signals from one of the patients is shown in figure 11 and its corresponding identified residue functions using both deconvolution methods are presented in the following figure.

Figure 14.
Estimated residue function using the ARMA and SVD deconvolution model for patient 6.

Note that the residue function falls monotonically to zero in ARMA while the SVD function presents fluctuations over the time. The blood flow range in the 18 patients varied from 84 to 370 ml/min/100g in ARMA and from 43 to 300 ml/min/100g using SVD.
The clinically acquired human MRI data were coherent with the simulations since the estimated blood flow applying ARMA was in average 1.21 times greater than the SVD method. All the identified rCBF are summarized in the following figure.



Figure 15.
Comparison of identified cerebral blood flow in 19 stroke patients using ARMA and SVD deconvolution.




5. CONCLUSIONS AND FURTHER WORK

5.1. CONCLUSIONS

The performances of the ARMA and SVD deconvolution method have been compared in ROI and pixel SNR, using two different sample rates. Each method had its advantages and weaknesses depending on the application environment. This is why it is not desirable to search the perfect or the best deconvolution method, without knowing the context. Instead the main question should reformulated as:

Under which circumstances is it better to apply one algorithm instead of another?
Or, Which is the most appropriate algorithm for some specific situation?

The ARMA deconvolution was generally closer to the true blood flow value, especially when the random noise corresponded to a ROI selection and for Rs=1Hz. The standard deviation in ARMA was generally greater that the SVD one. However, when the sample rate was halved, the increment of the SVD standard deviation was greater than ARMA. On the other hand, SVD showed to be less sensitive to SNR shifts. Consequently, SVD is a suitable technique when the sample rate is halved and if ROI or pixel selection is simultaneously being used.

The adaptive threshold deconvolution can be a useful method when the order of magnitude of the signal-to-noise ratio is unknown. Otherwise, this model is not advisable, since its execution time was much greater than the SVD deconvolution.

The simulations suggested that there is no difference in blood flow estimation in the two simulated volumes. In other words, the behavior of both deconvolution methods remained was similar or analogous when dealing with healthy white matter or healthy gray matter.

5.2. FURTHER WORK

The results from this study showed the influence of random noise in blood flow estimation and it was seen how strategies such as the ROI selection could optimize the rCBF identification. Further research could include regularization solutions to solve these kind of discrete ill-posed problems. Different regularization strategies would be compared in order to determine the most appropriate technique for this MR perfusion-imaging context.

It would be interesting to quantify the influence of image registration in the correct blood flow estimation.

The integration of both deconvolution methods into the software developed by the Cardiotools Creatis team is highly conceivable.

In this study, the ARMA deconvolution was applied using the first and second order for the poles and zeros models respectively. However, could the rCBF estimation be improved if the moving average (zeros model) order was modified?


REFERENCES


Wu O, Ostergaard L, Weisskoff RM, Benner T, Rosen BR, Sorensen AG. Tracer arrival timing-insensitive technique for estimating flow in MR perfusion-weighted imaging using singular value decomposition with a block-circulant deconvolution matrix. Magn Reson Med. 2003 Jul;50(1):164-74.

Neyran B, Janier M, Casali C, Revel D, Canet E. Mapping myocardial perfusion with an intravascular MR contrast agent : Robustness of deconvolution methods at various blood flows. Magn. Reson Med, 2002 ; vol 48, pp. 261-268.

Gobbel GT, Fike JR. A deconvolution method for evaluating indicator dilution curves.
Phys med Biol 1994;39:1833-1854.

Smith M, Lu H, Frayne R. Signal-to-noise effects in quantitative cerebral perfusion using dynamic susceptibility contrast agents. Magn. Reson Med 49: 122-128 (2003)

Wiart M, Imagerie par résonance magnetique (IRM) de la perfusion cérébrale. Modelisation de la cinétique d’un produit de contraste pour la quantification de la perfusion. Thèse de doctorat, Lyon : Université Claude Bernard Lyon 1; 2000, pp.266

Wiart M, Rognin N, Berthezene Y, Nighoghossian N, Froment JC, Baskurt A. Perfusion-based segmentation of the human brain using similarity mapping. Magn. Reson Med, 2001 ; vol 45, n°2, pp. 261-268.

Carme S, Modelisation de la perfusion quantitative en imagerie par resonance magnetique (IRM) cardiaque in-vivo. Thèse de doctorat, Lyon : Ecole doctorale : Electronique, Electrotechnique Automatique ;2004.

SAMPLE  SELECTION
How  large  should  a  sample  be  from  any  given

population? This question takes us into the mathematics
of  probability.  Do  not  worry,  you  will  not  have  to
understand  the  statistics  behind  national  surveys

produced by the likes of George Gallup. But perhaps a
quote  from  Mr.  Gallup  might  help  put  the  question  of
sample  size  in  perspective.  “Both  experience  and

statistical theory point to the conclusion that no major
poll  in  the  history  of  this  country  ever  went  wrong
because  too  few  persons  were  reached.”

Gallup conducted a number of experiments on the
effects  of  sample  size.  In  1936,  he  used  30,000  ballots
to ask the question: “Would you like to see the National
Recovery  Act  (NRA)  revived?”  The  first  500  ballots

showed  a  “no”  vote  of  54.9  percent.  The  complete
sample of 30,000 ballots returned a “no” vote of 55.5
percent.  In  other  words,  the  addition  of  29,500  ballots

to the first 500 ballots only made a difference of 0.6
percent.  (See  fig.  9-2.)
Through  the  mathematics  of  probability,  we  know
there is a real but unknown distribution of all possible

answers to a question. If we then know that our sample
is random (meaning that every person in our audience
is just as likely to get a survey as any other person) and
that our techniques are capable of obtaining a reliable
response (without bias) from each person, we will be
able to tell how representative the responses are.
For  the  purposes  of  this  chapter,  the  Sociology

Department at the University of West Florida supplied
a quick and easy formula often used in social science
research. It is shown in figure 9-3.
N> (t*SD/p)^2
where p is the accepted margin of error



http://www.tpub.com/content/photography/14129/css/14129_270.htm

METHODS

Monte Carlo simulations

To obtain realistic MR signals, representing the statistical data samples in the Monte Carlo simulation, in order to test the deconvolution techniques, the concentration functions were generated as follows:

Representative rCBF values were selected from the scientific literature. The rCBF range varying from 5 to 35 ml/100g/min in 5 ml/100g/min increments was chosen for a rCBV = 2%. Afterwards, the rCBF was used to calculate the MTT characteristic values, applying the central volume theorem. The same process was used for rCBV=4% except for that the range varied from 10 to 70 ml/100g/min, in 10 ml/100g/min increments. The simulated residue functions were generated as for each of the different rCBF values. The arterial input function (AIF) was simulated as a gamma-variate function


where u(t) is the step function (due to causality).
The concentration tissue response (Ct) was obtained by convolving the AIF with each residue function. SI functions were calculated assuming an exponential relationship between C(t) and the pixel gray-level:


where S0 is the baseline MR image intensity, k is a proportionality factor and TE is the echo time.

Additive white Gaussian noise was then used to blur the SI function samples with SNRtissue=6 or 18 dB and SNRAIF=15 dB, considered as typical values found in clinical practice. It is very important to notice, that the noise can not be added to the concentration functions since the existence of a non linear relationship between C(t) and the simulated signal intensities.

NUMBER OF SAMPLES
STATISTICAL ANALYSIS

ARMA VS SVD
SVDs VS SVD psvd
PATIENT DATA



===========


The number of simulated flow values

The mean percentage error was defined as the mean of the percentage difference between true and identified cerebral blood flow.

CONTENTS

1. Introduction

2. Theoretical background
2.1. Biological system modeling
2.1.1. Dynamic systems
2.1.2. The Stewart-Hamilton physiological model
2.2.Deconvolution methods
2.2.1. SVD
2.2.1.1. SVD
2.2.1.2. Adaptive threshold SVD
2.2.2. ARMA

3. Methods
3.1. Monte–Carlo simulations
3.1.1. ARMA vs. SVD
3.1.2. SVD vs. Adaptive threshold SVD

3.2. Clinically acquired human MRI

4. Results
4.1. Monte-Carlo simulations
4.2 Clinically acquired human MRI

5. Conclusions

6. References
INTRODUCTION

Cerebral and cardiovascular diseases caused by the obstruction of blood vessels, are one of the major causes of death in our modern society. The medical community has determined this class of disorders as an important public health threat.

Magnetic resonance imaging (MRI) technology has been widely used in clinical practice as a non-invasive technique to visualize the inside of the human body and to detect health-threatening situations in patients. In many cases, it is required to image the patient flowing blood, with the purpose to detect anomalies in the circulatory system. Perfusion MR imaging is a suitable method for this task, in which a contrast agent is introduced into the blood stream and then, it is digitally imaged as it travels across the tissue. MR perfusion imaging, is especially a good choice for diagnosis and follow-up studies of several injuries in the arterial system.

Consequently, the study field for new concepts and applications within MRI medicine is today a vast and dynamic area, as it was proved last year by Paul C. Lauterbur and Sir Peter Mansfield , in their discoveries concerning magnetic resonance imaging.
Figure 1.
A typical RM brain image sequence

When thrombosis is detected, it is also often desired to evaluate the regional cerebral blood flow (rCBF). This perfusion parameter quantifies the delivery of nutrients to the tissue and can determine whether the affected tissue is viable or necrosed, in order to take appropriate treatment decisions. However, in clinical practice there is still a need for robust diagnostic tools for determining perfusion parameters in environments affected by random electrical noise.

The purpose of this study was to determine the robustness of the SVD deconvolution technique and the ARMA method in rCBF estimation. Each algorithm was tested on N data sets, and we defined error measurements in order to compare the performance of the different methods. Computer simulations here were not intended to substitute practical data but are only applied as a first approach before analyzing the clinical samples.

4. RESULTS

4.1. Deconvolution performance criteria

MTT (s) 24 12 8 6 4,8 4 3,43
RBF(ml/100g/min) RBV=2% 5 10 15 20 25 30 35
RBF(ml/100g/min) RBV=4% 10 20 30 40 50 60 70
Table 1.
Simulated perfusion parameters. Mean transit time, regional blood flow and percentage regional blood volume.

The mean of the simulated samples, for each of the estimated rCBF, was calculated and subsequently compared with the true values, showed in the above table. To evaluate the performance of the deconvolution techniques, the percentage error (PE) and the standard deviation was calculated for the different rCBF and blood volumes. The mean percentage error, denoted MPE, was defined as follows:



where N is the number of simulated flows and

4.2. ARMA and SVD deconvolution

Sampling Rate: 1 Hz (1 sample per second)
In the case of non-noisy data, as expected, for a threshold (PSVD) of 30%, the SVD estimated flow was underestimated, the mean percentage error was of –34%. In contrast, using the ARMA deconvolution the true and identified cerebral blood flows were identical for both blood volume values, hence, the percentage error was of 0%.

Figure 4.
Comparison of performances without random noise of the ARMA and SVD deconvolution technique in estimating blood flow with a non-optimal selected threshold percentage

Non noisy-simulation signals were here used to show the influence of the Psvd percentage threshold selection. The underestimation in the SVD method was naturally due to the elimination of non-noisy information contained in the singular values matrix. Those elements are not supposed to be removed without the presence of random noise. Obviously, if a small (close to zero) Psvd percentage threshold was chosen, both results would have been the same. The influence of this parameter is shown in figure 5. Note that the appropriate Psvd for a high SNR tends to zero.

Figure 5.
Influence of the selection of the Psvd threshold percentage in estimating blood flows in a non-noisy environment.

For SNRtis = 18 dB and SNRaif = 15 dB the identified rCBF was underestimated with in the ARMA simulation, with a mean percentage error of –6.9% and mean standard deviation of 5.8. In SVD, there was also an underestimation; nevertheless, the absolute value of the mean percentage error was greater than the ARMA one.
Figure 6.
Standard deviation comparison for ARMA and SVD deconvolution techniques.
On the other hand, the standard deviation of the SVD was always considerably less than the ARMA one. This phenomenon was valid for both blood volumes as it is illustrated in figure 6.

When the tissular SNR was decreased to 6 dB without changing the arterial signal-to-noise ratio (SNRaif = 15 dB) the flow was in this case, remarkably overestimated using the ARMA model, with a mean percentage error of +86%. The identified CBF with the SVD deconvolution technique, on the contrary, remained underestimated with a mean percentage error of -34.5 %.
SNRtis=6dB





Sampling Rate: 0.5 HZ
(1 sample every 2 seconds)


The identified rCBF was not affected by undersampling in the ARMA deconvolution, hence, the error remained the same (i.e., PE=0%). The SVD method underestimated, again, the rCBF.
Sampling Rate: 0.5 HZ
(1 sample every 2 seconds)

SNRaif = 15 dB
SNRtissue=18 dB


After undersampling the rCBF is underestimated in both ARMA and SVD deconvolution. The MPE was less in ARMA.




Sampling Rate: 0.5 HZ
(1 sample every 2 seconds)
SNR=6 dB





When the SNRtissue was changed into 6 dB the ARMA estimation overestimated in this case the flow on the other hand, SVD remained underestimating.


Time widow
Non-noisy SI
Noisy SI






The perfusion parameters obtained with the ARMA deconvolution technique were more accurate than SVD in a non-noisy environment.

In average, The SVD deconvolution was more expensive than ARMA, in terms of execution time.


4.3. SVD and pSVD deconvolution

The standard deviation did not vary significantly with the volume changes


However it is important to keep in mind that the execution time in the adaptive deconvolution was always greater than SVD. That is quite natural, considering that the pSVD algorithm tested 4 different matrix thresholds, recalculating at each time by deconvolution, the residue function. Nevertheless, one advantage of the pSVD is its standard deviation



===
Was less sensitive
Demonstrates the variability of

2. Theoretical background

2.1. Biological system modeling

2.1.1. Dynamic systems

The physiological behavior of some organs can be modeled using the dynamic systems theory. The heart, veins, arteries and capillaries may be treated as a biological system, mathematically representable, where variables and parameters interact regularly. This representation provides a powerful tool to estimate variables when some parameters are unknown. Another advantage of the dynamic systems theory is that we can handle various -apparently different- physical parameters (e.g., temperature, flow, concentration and voltage) without modifying the equation’s essence. In recent years, the study field for new concepts such as the biochemical circuit theory analysis is today a fundamental field of interest for many laboratories.

The physical elements treated in this medical context, are the perfusion parameters, i.e., regional cerebral blood flow rCBF, regional cerebral blood volume rCBV, mean transit time MTT. However, the same parameters and relationships can be deduced and applied in a different area such as the myocardial circulatory system. In the published literature, the notation is quite similar except for that the letter M (myocardial) replaces the letter C (cerebral), in the perfusion parameter abbreviations.

The human circulatory system is divided into two main parts with the heart acting as a double pump. This organ pumps of a special kind of fluid, which is the blood. The across variable in this case, the blood flow, leaves the left side of the pump (heart) and travels through arteries which gradually divide into capillaries. In the capillaries, the nutrients are released to the body cells. The blood then travels in veins back to the right side of the pump, and the whole process begins again.

One of the most common ways in which, this kind of process is analytically described is the black-box model. The biodynamic system has, therefore, the following components:

 Parameters: the rCBF and rCBV.
 Dependent variables: the arterial, tissular and venous concentration denoted Ca, Ct and Cv respectively.
 Independent variable: the time.

Figure 2. The blood circulation as a dynamic system

2.1.2. The Stewart-Hamilton physiological model

The Stewart-Hamilton model applied to MRI acquisition techniques can estimate perfusion parameters such as the regional blood flow rBF, mean transit time MTT or regional blood volume rBV, describing the tissue condition.

The concentration of tracer within a tissular volume, at a given time t, during the passage of a bolus injection of a contrast agent, is given by:



Where denotes the tissue density (i.e., the tissue mass per unit volume), Cartery is the arterial concentration and R(t) is the residue function and represents the tissular retention of the contrast agent. If this biodynamic system has a single fluid storage element hence, the residue function can be treated as the response of a first-order circuit. The general solution for this differential equation is of the form:



In this particular case we have: is the mean time taken by the blood to traverse the system, from the arterial entrance to the vascular exit and k=1.

Let

Therefore



Assuming the tissular density being close to the water density, i.e., ,
then:


Once the regional blood flow is estimated, the unknowns perfusion parameters, can be computed from the central volume theorem (Stewart, 1984; Meier and Zierler, 1954):



where rCBV and MTT are the regional cerebral blood and mean transit time, respectively and



However, the concentration is measured from the MR digital images, therefore, the continuous time form of the equations is not considered. The Stewart-Hamilton model becomes a discrete convolution:



Moreover, the discrete convolution can be viewed as a matricial product:



Where N is the number of samples taken in the RM sequence.

This product can be compactly denoted as:



Where is unknown, then:


Hence, the deconvolution process becomes a linear algebra problem, specifically, a matrix inversion problem. However, this process, when working with noisy signals, could sometimes lead to determinants close to zero, the matrix A is said to be “ill-conditioned”, and consequently, unexpected results may be obtained.

Nevertheless, the concentration of the contrast agent is not directly measured in clinical practice; instead, MR images are used to determine signal intensity variations over time, which are related with the change in the spin-spin relaxation, also known as the transverse relaxation . The SI functions are obtained by selecting a pixel or a region of interest (ROI) in the MR gray-scale digital images.

Previous studies have shown the existence of a linear relation between the measured concentration C(t) and . On the other hand, the transverse relaxation and the signal intensity have an exponential dependence, therefore, we write:



Where TE is the sequence echo time, S(t) is the signal intensity over time and S0 is the baseline MR intensity.

Consequently, by measuring SItissue and SIaif we can estimate the regional blood flow using a deconvolution technique and evaluating the expression r(t) at zero.

2.2. Deconvolution methods

2.2.1.1. General Singular Value decomposition (SVD)

The singular value decomposition is a widely used technique to solve ill-conditioned problems with several applications, (e.g., in image compression, watermarking, image filtering). The general SVD statement can be expressed as:

Every real matrix A can be decomposed into a product of three matrices of the form:



Where, U and V are orthogonal matrices and with . Therefore, S is a diagonal matrix whose elements are the singular values of the original matrix.

Consequently, the inverse can be expressed as:



Where W is also a diagonal matrix of elements

Depending on the signal-to-noise ratio of the signal intensity function, a tolerance threshold PSVD is set. The values of wi corresponding to values where si is less than PSVD are set to zero. Typically the threshold is given as a percentage of the greater singular value . Generally the threshold value varies between 20 and 30% and the general principle is the higher the noise, the higher the PSVD.


2.2.1.2. Singular Value decomposition adaptive threshold (pSVD)

Liu et al. (2) showed that the PSVD selection had a significant influence in the shape of the residue function and an apparently inaccuracy in the rCBF estimation. An oscillation index O, was proposed by Østergaard et al. (1) to measure the distortion in the residue function as follows:


Where f is the scaled estimated residue function, fmax is the maximum amplitude of f, and L is the number of sample points.

This oscillation index is the discrete form of the convolution product between the residue function and the second derivative of a Gaussian distribution. The use of this digital filter, which is the same as calculating the sum of the absolute value of the second derivative of f(t), determines the change in the curvature over the time.
Figure 2. The Gaussian distribution (a) and its second derivative (b)


2.2.2. ARMA

The estimation of the residue function can be treated as a spectral parametric analysis problem. It is known that when a linear time-invariant (LTI) system is excited by the Dirac's delta function, the output is the transfer function of the system.

Figure 2. Discrete time system - ARMA model representation

The discrete time filtering theory states that the impulse response of a linear system h(n) can be modeled using one of the following frequency responses:
The “zeros model” or moving average (MA):



The “poles model” or autoregresive (AR):



The “zeros and poles” model (ARMA)




Which can be written as:

The aim of the ARMA technique modeling is to find the transfer function H(z), whose impulse response h(k) approximates the scaled tissue response r(k), such that the sum of the squared error, denoted e, is minimum:

Considering the system illustrated in figure 2, and using the ARMA model, the output y(n) can be computed as



This mathematical statement generates the following set of simultaneous equations:



Which can be compactly expressed, in matrix notation as

,

where Y represents the output matrix, A is the input-output matrix and is the ARMA parameter matrix (coefficients ak and bk).

Therefore, the general ARMA-model matricial representation can be written as



To solve the linear equations system , an estimation square error is defined as


where

The minimization of the error by the least mean square process gives the coefficients ai and bi, and consequently the searched solution.



where

To minimize this equation the SE must be differentiated with respect to and equated to zero.

Solving the equation we have


Where represents the unknown ARMA parameters.