WO2005104936A1

WO2005104936A1 - Method and system of obtaining improved data in perfusion measurements

Info

Publication number: WO2005104936A1
Application number: PCT/AU2004/000821
Authority: WO
Inventors: Qing Yang
Original assignee: Apollo Medical Imaging Technology Pty Ltd
Priority date: 2004-04-30
Filing date: 2004-06-23
Publication date: 2005-11-10
Also published as: EP1635703B1; EP1635703A1; US8855985B2; US20110118615A1; EP1635703A4; US20060083687A1; US8285490B2

Abstract

A method of deriving blood perfusion indices for a region of interest (ROI) of a subject, comprising the steps of administering a contrast agent to the subject during a dynamic imaging scan, converting signal intensity data from raw images of the scan into contrast agent concentration data, deriving parameters from the contrast agent concentration data using at least one transport function that accounts for delay and dispersion of the contrast agent, and calculating the blood perfusion indices from the derived parameters.

Description

Method and System of Obtaining Improved Data in Perfusion Measurements

Field of the Invention This invention relates to a method and system of obtaining improved data in blood perfusion measurements, and more particularly to a method and system of deriving blood perfusion indices for a region of interest of a subject.

Background to the Invention The process of measuring blood flow within a body of a subject non-invasively is useful in diagnosing and treating the subject. This is particularly the case where a part of a subject or patient, such as a tissue or organ, suffers from ischaemia due for example to a stroke. Determining perfusion indices including the blood flow through such a tissue or organ can provide important information to a physician in order to determine an appropriate treatment regime for the patient. A number of systems pertaining to blood flow information have been disclosed. In general, the systems involve a contrast agent delivered as an intravascular bolus during a dynamic imaging session such as computerized tomography (CT), nuclear medicine (NM) or magnetic resonance imaging (MRI). The temporal profile of the image intensity in a pixel or region of interest (ROI) reflects the characteristics of the contrast agent hence the blood passing through the vasculature. The typical method of obtaining quantitative perfusion indices involves several steps including: (a) convert the signal intensity profile to the contrast concentration profile depending on the type of imaging modality; (b) measure the arterial input function (AIF) from a feeding vessel to the tissue of interest; (c) measure the tissue profile; (d) extract the tissue impulse residue function (IRF) from the AIF and tissue profile using deconvolution; (e) calculate quantitative perfusion indices including blood flow (BF), blood volume (BN) and mean transit time (MTT) using the IRF. Furthermore, the tissue IRF contains information about the flow heterogeneity associated with dispersion of blood transit time through capillaries, which is an important factor determining the efficacy of oxygen delivery to tissue. However, in the case of major vessel disease, such as acute stroke or carotid artery stenosis, the measured AIF is often associated with a delay and dispersion before it reaches the tissue of interest, and causing overestimation of the MTT and underestimation of the BF. In United States Patent No. 5,190,744 a contrast agent is injected into a patient for the purpose of detecting blood flow abnormalities. This disclosure describes in some detail the different types of agents that can be used and the administration of those agents into the patient. However results of the perfusion process may not be evaluated until some time after the initial injection of the contrast agent and is thus not a real time process. Furthermore there is no disclosure of obtaining quantitative data relating to blood flow and blood volume which can assist a physician to make a relatively quick and accurate diagnosis and decide on what steps can be taken to treat the patient. More particularly this document does not account for any delay or dispersion of the contrast agent in an initial bolus injection. In US Patent No. 6,542,769 there is disclosed an imaging system associated with MRI whereby a bolus containing optical and MRI contrast agents is administered to a patient in order to determine perfusion indices. It uses an optical contrast agent which is injected into the patient and is used to define the arterial input function. The optical contrast is injected so as to overcome the problem of the signal intensity of the vasculature not being proportional to the amount of contrast agent with MRI. A disadvantage of measuring the signal change in arteries using MRI is that it does not provide a true indication of the contrast volume as MRI depends upon electromagnetic fields that are altered due to the contrast agent. By using an optical contrast agent the invention disclosed in this document tries to overcome these disadvantages. Again there is no taking into the account the delay and dispersion associated with the bolus progressing through the artery selected and through the tissue or organ in the region of interest. The present invention seeks to substantially overcome, or at least ameliorate, any one or more of the abovementioned disadvantages. Summary of the Invention According to a first aspect of the invention there is provided a method of deriving blood perfusion indices for a region of interest (ROI) of a subject, the method comprising the steps of: administering a contrast agent to the subject during a dynamic imaging scan: converting signal intensity data from raw images of the scan into contrast agent concentration data; deriving parameters from the contrast agent concentration data using at least one transport function that accounts for delay and dispersion of the contrast agent; and calculating the blood perfusion indices from the derived parameters. The transport function may represent a probability distribution function of transit times of the contrast agent through the subject. More particularly, the method may further comprise using a first model to represent an arterial transport function h_a(t) through a vessel leading to the ROI, and using a second model to represent a tissue transport function h_s(t) through the ROI. The transport function preferably accounts for the delay and dispersion of the contrast agent simultaneously. The method may further comprise selecting an arterial input function AIF_a(t) in the vessel, preferably an artery, leading to the ROI by searching pixels taken of the contrast agent concentration data. The method may further comprise measuring the contrast agent concentration C(t) remaining in the ROI. The method may further comprise representing h_a(t) using a gamma-variate function (GNF) in the first model such that:

A_λ M = 0 (t < t_λ )

where A - σ_x "TCL+α,) , T ( ) ≡ x ^{a ~} e^{~ x}dx is the Gamma function, ti is J 0 the time taken for the contrast agent to move from the initial measurement of AIF_a(t) to a vessel, preferably an artery, at the entry to the ROI, o and α₁ are related to the mean transit time and dispersion of h_a(t). The method may further comprise estimating h_a(t) after deriving values for parameters ti and σι and setting a_λ =0 through the equation:

The method may further comprise determining an estimate for the arterial input function AIF_t(t) of the vessel at the entry to the ROI through the equation: AIF_t(t) = AIF_a (t) ® h_a(t) ≡ {' AIF _a (τ)h_a (t - τ)dτ JO The method may further comprise determining an estimate of blood flow F and an estimate of the tissue IRF R_e(t) from the deconvolution of: C(t) = (F_t/ky) ΛIF_t(t) ®R_e(t)

where kH l-H_a)/(l-H_t) is a correction constant taking into account different values of arterial hematocrit H_a and tissue hematocrit H since the contrast agent remains in the extracellular fraction of blood (plasma). The hematocrit is the volume fraction of cells in the blood, which has a typical value of H_a « 0.45 for large vessels such as the artery and a value of H_t « 0.25 for small vessels such capillaries in tissue. The method may further comprise determining an estimate for the tissue transport function h_e(t) from the estimated R_e(t) using: h_e(t) = - ^- R_e(t) at

The method may further comprise determining a rise time and a mean transit time of h_e(t) in order to determine values of parameters α₂ and σ₂ by assuming t₂=0_j or determining a peak height and a mean transit time of h_e(t) in order to determine values of parameters σ₂ and t₂ by assuming ₂=0, where α_2, σ₂ and t are parameters related to the mean transit time and dispersion of h_e(t). The method may further comprise representing a simulated transport function h_s(t) using a GNF in the second model such that:

(t - t₂ y* e -(. -.₂ ) / σ₂ (t ≥ t₂ ) K (t) = 0 (t < t₂ )

where A₂ = σ^ ²T( + ₂) , t₂, σ₂ and α are parameters related to the mean transit time and dispersion of h_s(t) through the ROI. The method may further comprise estimating h_s(t) using the derived values for parameters α and σ₂ by setting t₂=0 through the equation: h_s(t) = —t'^he^{' ,σ}* (t > 0) or using the derived values for parameters σ₂ and t₂ by setting ₂=0 through the equation: l_ _e-( -t₂ )/_{σ 2} (^ > £₂ ) σ₂ A, (0 = 0 (t < t₂ )

The method may further comprise determining a simulated tissue IRF R_s(t) by: R_s (t) = l - ^h_s (τ)dτ The method may further comprise determining a simulated contrast agent concentration C_s(t) as: (t)= (F_t/kA) AIF_t(t) ®R_s(t) The method may further comprise fitting the simulated C_s(t) to C(t) using a least squares method according to: S = ∑(C(t) - C„(t))² . The method may further comprise the step of optimising the parameters F_t, tι_; σ_1; oti, σ₂, α₂ and t₂by minimizing S iteratively. The method may further reduce the number of adjustable parameters by fixing

leading to five adjustable parameters. The method may further reduce the number of adjustable parameters by fixing a relative dispersion,

of h_a(t) resulting in σi dependent on ti, hence leading to four adjustable parameters. The method may further comprise calculating quantitative blood perfusion indices from the optimized parameters of F_t, t_1; σ_1; α_l5 σ , α₂ and t₂. The perfusion indices may include any one or more of blood flow, blood volume, mean transit time, arterial delay time, arterial dispersion time or relative arterial dispersion, tissue dispersion time or relative tissue dispersion. Preferably the ROI is a tissue. The ROI may be a pixel or a plurality of pixels in a tissue. The scan may be any one of CT, MRI or NM. In the brain, many cerebral arteries are small subjecting to partial voluming. The method may further comprise determining a venous input function (VIF_a(t)) from a draining vein to estimate an AIF_a(t) where a selected artery has partial voluming, the vein being larger than the artery. The method may further comprise the step of determining the profile of a venous input function (NIF_a(t)) from a large draining vein. The AIF_a(t) may then be scaled up to have the same first-pass bolus peak area as the NIF_a(t) to minimize partial voluming (PN) effect from the AIF_a(t). The first-pass AIF_a(t) and NIF_a(t) profiles can be obtained by fitting the profiles to gamma-variate function (GNF) profiles respectively to remove contrast recirculation effects. The method may further comprise determining a simulated tissue IRF R_s(t) in the case that the contrast agent does not always remain in the vascular system, such as in a tumour in the subject in order to determine blood perfusion indices and permeability indices using: R, (0 = 1 - fa (τ)dτ + Ee^'1* J λ, (τ)e^kτ dτ r i (t - t₂ y* e -( t-t₂ ) l σ₂ (t ≥ t₂ ) A where ^ (0 = 1 0 (t < t₂ )

E is the extraction fraction of the tracer in the blood stream that leaks out of the vessel into tissue, and the tracer clearance rate constant k=E*F_t /N_e is a rate constant at which the leaked contrast agent diffuses back into the blood stream and leaves the tissue, N_e is volume fraction of the extravascular and extracellular space (EES). The permeability surface area product PS can be determined by RS = -F_t ln(l - E) . The method may further comprise the step of repeating the entire process (except for selecting the AIF and/or NIF) on a pixel-by-pixel basis to produce quantitative maps of the perfusion indices for further analysis and presentation. According to a second aspect of the invention there is provided computer program means for deriving blood perfusion indices for a region of interest (ROI) of a subject by directing a processor to: retrieve raw image data from a dynamic imaging scan of the subject after a contrast agent is administered to the subject; convert signal intensity data included in the retrieved raw image data into contrast agent concentration data; derive parameters from the contrast agent concentration data using at least one transport function that accounts for delay and dispersion of the contrast agent; and calculate the blood perfusion indices from the derived parameters. The computer program means may further direct the processor to select an arterial input function AIF_a(t) in the vessel, preferably an artery, leading to the ROI by searching pixels taken of the contrast agent concentration data. An optimal AIF_a(t) may be selected on the basis of early arrival and high and narrow peak for the arterial input function. The program means may further direct the processor to measure the contrast agent concentration C(t) remaining in the ROI. The program means may further direct the processor to estimate the arterial transport function through a vessel leading to the ROI, h_a(f), using a GNF in a first model such that:

and thereafter estimating h_a(t) after deriving values for parameters ti and σi and setting a_γ =0 through the equation:

MO o (* < ⁾

The program means may further direct the processor to determine an estimate for the arterial input function AIF_t(t) of the vessel at the entry to the ROI through the equation: AIF_t(t) = AIF _a (t) 0 h_a(t) where ® is the convolution operator.

The program means may further direct the processor to determine an estimate of blood flow F_t and an estimate of the tissue IRF R_e(t) from the deconvolution of: C(t) = (F_t/km AIF_t(t) ® R_e(t)

where k_H=(l-H_a) (l-H_t) is a correction constant taking into account different values of arterial hematocrit H_a and tissue hematocrit H since the contrast agent remains in the extracellular fraction of blood (plasma). The hematocrit is the volume fraction of cells in the blood, which has a typical value of H_a « 0.45 for large vessels such as the artery and a value of H_t « 0.25 for small vessels such capillaries in tissue. The program means may further direct the processor to determine an estimate for the tissue transport function h_e(t) from the estimate R_e(t) using: h_e(t) = - ~ R_e(t) dt

The program means may further direct the processor to determine a rise time and a mean transit time of h_e(t) in order to determine values for parameters α₂ and σ₂by assuming t =0_; or to determine a peak height and a mean transit time of h_e(t) in order to determine values for parameters σ₂ and t₂ by assuming α =0, relating to mean transit time and dispersion of h_e(t). The program means may further direct the processor to estimate a simulated transport function h_s(t) using a GNF in a second model such that:

~(t - t₂y -{t- ) l σ₂ (t ≥ t₂ ) A₂ Λ. (0 0 (t < t₂ )

where A₂ = σ^¹ "¹ T(l + a₂) , t₂, σ₂ and α are parameters related to the mean transit time and dispersion of h_s(t) through the ROI. Thereafter the program means may direct the processor to estimate h_s(t) using the derived values for α₂ and σ₂ by setting t₂=0 through the equation:

or using the derived values for σ₂ and t₂ by setting ₂=0 through the equation:

The program means may further direct the processor to determine a simulated tissue IRF R_s(t) by:

The program means may further direct the processor to determine a simulated contrast agent concentration C_s(t) as: C_s(t) = (F_t / AIF_t(t) ®R_s(t) The program means may further direct the processor to fit the simulated C_s(t) to

C(t) using a least squares method according to: S = ∑(C(t) - C.(t))²

The program means may further direct the processor to optimize the values F_t, t_1; σ_1; αi, σ₂, α and t₂by minimizing S iteratively. The program means may direct the processor to reduce the number of adjustable parameters by fixing

leading to five adjustable parameters. The program means may direct the processor to further reduce the number of adjustable parameters by fixing a relative dispersion,

of h_a(t) resulting in o_\ dependent on t_1; hence leading to four adjustable parameters. The program means may further direct the processor to calculate quantitative blood perfusion indices from the optimized values of parameters F_t, t_1; σ_i; αi, σ₂, α₂ and t₂. The perfusion indices may include any one or more of blood flow, blood volume, mean transit time, arterial delay time, arterial dispersion time or relative arterial dispersion, tissue dispersion time or relative tissue dispersion. Preferably the ROI is a tissue. The ROI may be a pixel or a plurality of pixels in a tissue. The scan may be any one of CT, MRI or NM. The program means may further direct the processor to determine a venous input function (NIFa(t)) from a draining vein to estimate an AIFa(t) where a selected artery has partial voluming, the vein being larger than the artery. The program means may further direct the processor to determine the profile of a venous input function (NIF_a(t)) from a large draining vein. The AIF_a(t) may then be scaled up to have the same first-pass bolus peak area as the NIF_a(t) to minimize partial voluming (PN) effect from the AIF_a(t). The first-pass AIF_a(t) and NIF_a(t) profiles can be obtained by fitting the profiles to gamma-variate function (GNF) profiles respectively to remove contrast recirculation effects. The program means may further direct the processor to determine a simulated tissue IRF R_s(t) in the case that the contrast agent does not always remain in the vascular system, such as in a tumour in a subject in order to determine blood perfusion and permeability indices using: R_s (0 = 1 - J_o\ (τ)dτ + Ee-^M \ _s (τ)e^kτdτ

(t - t₂ )^aa^> ₂ e -'t-t₂ ) l σ₂ (t ≥ t₂ ) A, where h_s (t) = 0 (t < t₂ )

E is the extraction fraction of the tracer in the blood stream that leaks out of the vessel into tissue, and the tracer clearance rate constant k=E*F_t /N_e is a rate constant at which the leaked contrast agent diffuses back into the blood stream and leaves the tissue, V_e is volume fraction of the extravascular and extracellular space (EES). The program means may further direct the processor to calculate the permeability surface area product PS byRS = -R in(l - E) . The program means may further direct the processor to repeat the entire process (except for selecting the AIF and/or NIF) on a pixel-by-pixel basis to produce quantitative maps of the perfusion indices for further analysis and presentation. According to a third aspect of the invention there is provided a system of deriving blood perfusion indices for a region of interest (ROI) of a subject, the system comprising: scanning means for providing a dynamic image scan of the subject during which a contrast agent is administered to the subject; processor means linked to the scanning means for retrieving raw image data from the scan; the processor means further: converting signal intensity data included in the retrieved raw image data into contrast agent concentration data; deriving parameters from the contrast agent concentration data using at least one transport function that accounts for delay and dispersion of the contrast agent; and calculating the blood perfusion indices from the derived parameters.

Brief Description of the Drawings Preferred embodiments of the invention will hereinafter be described, by way of example only, with reference to the drawings wherein: Figure 1 is a side view of the head of the subject indicating flow of the contrast agent through a region of interest, such as a tissue; Figure 2 is a block diagram showing a communications network including a number of scanners linked to a data storage system and a processing system; Figure 3 shows various graphs against time at different parts of the subject's head as the contrast agent traverses the region of interest and the input artery; Figure 4 shows a plot whereby an input arterial profile for a small artery exhibiting partial voluming is scaled up based on a vein profile; and Figure 5 is a flow diagram showing steps performed by a computer program to obtain values for the blood perfusion indices such as BF, BN and MTT.

Detailed Description of Preferred Embodiments The present invention is particularly applicable to CT, MRI and MΝ imaging systems. A bolus of contrasting agents is introduced via a needle into a patient at, for example, the arm of the patient. However the bolus can be input to any other part of the patient. A region of interest (ROI) may be a tissue 6 in a part of the patient's brain as shown in Fig. 1. Alternatively, the ROI may be a pixel or a plurality of pixels, where many pixels represent a calculated image to produce one or more perfusion maps. Blood circulating throughout the patient will contain the contrast agent and in particular may be delivered to the tissue 6 via artery 8 and the blood flowing through the tissue 6 is returned to the heart via vein 10. Raw data and/or images collected by a scan, such as from a CT scanner 20, MRI scanner 30 or ΝM scanner 35 are forwarded to a data storage system 40 in the form of a Picture Archiving Communications System (PACS) in Fig. 2. A computer program operating on a processor 50, in the form of a computer, is used to retrieve the various images or raw data from any one of the scanners 20, 30 or 35 or from the data storage system 40. The program then processes those images to provide an improved data set for a clinician to use, particularly in relation to perfusion indices including blood flow, blood volume, mean transit time, arterial delay time, arterial dispersion time or relative arterial dispersion, tissue dispersion time or relative tissue dispersion.. The computer program need not reside on computer 50, but may reside in a console computer linked to any one of the scanners 20, 30 or 35. Alternatively the program may reside in a workstation (stand-alone or in a system) or in the PACS 40. In order to select an appropriate arterial input function, AIF, various images (slices) from a scan are analysed to identify a major artery of interest. In CT the signal changes are directly proportional to the contrast agent concentration profile. However in MRI a mathematical conversion is used in order to convert the measured signal time- curve into contrast agent concentration profile. From the raw data retrieved, the program stored in the memory of system 50 automatically calculates the contrast concentration based on the measured signal intensities of the contrast agent. It then searches all pixels to find the optimal AIF (or NIF) based on the criteria of early arrival, high and narrow peak for AIF, and later arrival, high and broad peak with maximum peak area for NIF. The program displays the searched AIF and NIF pixels on the corresponding images and plots the AIF and NIF time-curves. A user may further select a particular pixel while dynamically viewing its profile against the selected AIF in order to confirm the best arterial input function. A better arterial pixel can be saved to replace or average with the saved AIF and then the user may "click" on further pixels in order to compare the further pixels with the updated AIF until the user is satisfied with the selected AIF. The selection of the best AIF may be done through different slices with the effort to minimize partial voluming (PN) where only a part of the artery is contained in the pixel. Similarly, the best NIF profile can be confirmed by the user. Referring to Fig. 5, the computer program at step 100 thus retrieves raw data and/or images from any one of the scanners 20, 30, 35 or PACS 40, including the signal intensities containing information of the contrast agent. At step 102 the program calculates the contrast agent concentration based on the signal intensities. It then plots the contrast agent concentration profile C(t) against time at step 104. Where the data is retrieved from an MRI scan, the signal intensities are converted mathematically to obtain C(t) at step 106. At step 108 the program searches pixels taken from the plots to find an optimal AIF (NIF) based on given criteria such as arrival times and peaks. At step 110 the program displays the searched pixels of the AIF (NIF) and plots these as a function of time at step 112. The best pixel(s) to date are then stored in memory means, such as located at computer 50, at step 114. A decision is made at step 116 to determine if the optimal pixel has been found and stored, which decision can be made by the user. If an optimal pixel has not been found, the program keeps reverting to step 118 to find a better pixel than the pixel stored, which is subsequently stored at step 114. When an optimal pixel has been found the process moves to step 120, to be described hereinafter. The amount of contrast agent passing through the tissue 6 may then be measured by the computer program, the contrast agent concentration being represented as C(t). Thus two known profiles are used by the computer program, one for the concentration of the contrast agent C(t) and the other for AIF_a(t), being the arterial input function in the vessel (artery) leading to the ROI, against time. By knowing these two particular profiles the tissue blood flow F and tissue impulse residue function (IRF), R(t), can be derived from a deconvolution of the equation

AIF_a(t)®R(t), where k_H=(l-H_a)/(l-H_t) is a correction constant taking into account different values of arterial hematocrit H_a and tissue hematocrit H_t since the contrast agent remains in the extracellular fraction of blood (plasma). The hematocrit is the volume fraction of cells in the blood, which has a typical value of H_a « 0.45 for large vessels such as the artery and a value of H_t ∞ 0.25 for small vessels such capillaries in tissue. In other words the concentration of the contrast agent is derived by a convolution of the arterial input function and the tissue IRF multiplied by the tissue blood flow. This is the case where there is no delay or dispersion so that the selected AIF_a(t) from a major artery is taken to be the same as the AIF (t) directly feeding the tissue. In practice however it is difficult to measure the arterial input function at the input to different tissues of interest. This is due to the practicalities that, arteries directly feeding the tissues are usually small in size and subject to a substantial partial voluming effect. In the case of major vessel disease such as acute stroke or carotid artery stenosis, the AIF selected from a major artery is often associated with a delay and dispersion before it reaches the abnormal tissue of interest. As an example reference is made to Fig. 3 where the arterial input function is measured in artery 8 resulting in the graph of Fig. 3(A). It can be seen from the graph that there is a time t_a taken from injection for the contrast agent to arrive at the point where the arterial input function is measured in artery 8. It results in a narrow 'pulse' having a large amplitude. Then in Fig. 3(C) there is shown the arterial input function if measured at the tissue 6 input artery designated by 60. It can be seen that the graph has dispersed somewhat or broadened, as well as involving a time delay ta in traversing the smaller artery 62 where a vessel disease such as stroke or stenosis may occur. Other normal small arteries supplying different tissues may have little delay and dispersion. Therefore it is practically useful to select a normal large vessel such as the internal carotid artery (ICA), middle cerebral artery (MCA), anterior cerebral artery (ACA) or posterior cerebral artery (PCA) through multiple slices in the head and neck of the patient. The transit of the contrast agent from the artery 8 through the artery or arteries 62 up to the entry point 60 of the tissue 6 is described by the vascular transport model (NTM) using the delay t and an estimated dispersion parameter to derive a simulated function AIF_t(t) at the tissue input artery 60. The next part of the transit of the concentration of the contrast agent is described by the tissue perfusion model (TPM) where the contrast agent traverses across the tissue 6 from an input 60 to an output 64. The measured contrast concentration profile C(t) represents the contrast agent remaining in the tissue 6 as represented by the curve shown in Fig. 3(E) and the tissue blood flow Ft and impulse residue function (IRF) R_e(t) can be estimated using a model-free deconvolution technique such as the singular value decomposition (SND) method. However, such deconvolution is sensitive to noise, which may produce some mathematical solutions of R_e(t) but without any physiological meaning. Further, the estimated F_t and R_e(t) may not be accurate due to uncertainties associated with unaccounted delay and dispersion effects. It is desirable to use a constrained deconvolution process using a model derived IRF R_s(t) with a typical shape as shown in Fig. 3(D). The estimated R_e(t) can be used to derive parameters for R_s(t). A simulated tissue contrast concentration curve derived using convolution as C_s(t)=(F_t/kn)AIF_t(t) ®R_s(t) can be fitted to the measured C(t) curve by optimizing the model parameters through an iterative least square method. The ga ma- variate function (GNF), represented by equation (1) below, has been generally used to describe the temporal profile of contrast during blood circulation through the vascular system.

In order to account for the delay and the dispersion through the artery 62 the computer program employs a first model of GNF to represent a vascular transport function as

Where A_γ = e ^xdx is

the Gamma function, ti is the time taken for the initial AIF_a(t) measured from artery 8 to arrive at artery 60 and σi and αi are related to the mean transit time and dispersion of h_a(t). If setting

equation (2) becomes,

-{t- ) l σ₁ (t ≥ ) σ. M0 = (3) 0 (* < *ι )

An example plot of h_a(t) versus time is shown in Fig. 3(B). The value of ti can be estimated by t_d as the delay between C(t) in Fig. 3(E) and AIF_a(t) in Fig. 3(A). Starting with an estimate

the mean transit time of h_a(t) is ti + σi and a relative arterial dispersion is defined as β_t = σi / (ti + σ^ ranging from 0 to 1. A relative dispersion value of βi = 12% is chosen based on previous measurements of dispersions typical for arteries (12%), vein (30%) and whole organs (40%). Thus an initial estimate of σi can be obtained as tiβi / (1-βι). Referring again to Fig. 5, at step 120 the computer program applies the GNF to represent h_a(t) in a first model. At step 122 an estimate of ti is made from the plots of C(t) and AIF_a(t). Then at step 124, the program estimates σi using ti and a relative dispersion value βi assuming αι=0. The process then moves to step 126. With the estimated ti, i and σi values, the estimated h_a(t) function in equation (2) can be calculated by the computer program at step 126 to describe the arterial transport by AIF_t(t) = AIF_a(t) ®h_a(t) (A)

where ALF_t(t) is the arterial input function at the input to the tissue designated by 60, AIF_a(t) is the initial AIF at artery 8 and ® is the convolution operator. The contrast concentration profile in the tissue of interest can be further determined by

C(t) = (F_t/kn) AIF_t(t) ®R_e(t) (5) where k_H=(l-H_a)/(l-H) is a correction constant taking into account different values of arterial hematocrit H_a and tissue hematocrit H_t because the contrast agent remains in the extracellular fraction of blood (plasma). The hematocrit is the volume fraction of cells in the blood, which has a typical value of H_a « 0.45 for large vessels such as the artery and a value of H « 0.25 for small vessels such capillaries in tissue. With measured C(t) and the model derived AlF_t(t), an estimate of F_t and R_e(t) can be obtained using a model-free deconvolution technique such as the singular value decomposition (SND) method. The deconvolution is very sensitive to noise, which may produce some mathematical solutions of R_e but without any physiological meaning. Further, the estimated F_t and R_e(t) may not be accurate due to uncertainties associated with the initial estimate of t_1; αi and σι values. It is desirable to use a constrained deconvolution process using a model derived IRF with a typical shape as shown in Fig. 3(D). Again referring to Fig. 5, the computer program stored in memory of the computer 50 directs the computer at step 128 to calculate an estimate for AIF_t(t) from the convolution of AIF_a(t) and h_a(t) in equation (4) and at step 130 to calculate an estimate for F_t and R_e(t) from equation (5). A more realistic (simulated) profile of the tissue IRF can be provided by the second model of GNF, which describes the tissue transport function as

Where A₂ = σ₂ ²r(l + a₂) , t₂ , σ₂ and α₂ are parameters related to the mean transit time and dispersion of h_s(t) through the tissue. If assuming t₂=0, equation (6) becomes h_s(t) = —t^a*e-^llσ2 (t≥ O) (7a) ÷2 Several typical characteristic parameters of h_s(t) are determined as

Peak rise time (RT) = σ₂ α₂ Mean transit time (MTT) = σ₂ (1+ α₂) (7b )

Alternately, if assuming α₂=0, equation (6) becomes

Several characteristic parameters of h_s(t) are determined as

Peak height (PH) = l/σ₂ Mean transit time (MTT) = t₂+σ₂ (8b)

The relationship between the tissue IRF R(t) and transport function h(t) is

R(t) = l - ['h(τ)dτ = h(t) = -^ - (9) ^o dt

Since h(t) is a probability density function, R(t) is a positive, decreasing function of time with a property of with R(0) = 1 as shown in Fig. 3(D). From the estimated R_e(t) profile deconvolved from equation (5), an estimated transport function h_e(t) can be derived as h_e(t)= - dR_e(t)/dt. The peak rise time and mean transit time of h_e(t) can then be calculated and used to estimate σ₂ and α₂ using equation (7b ) or to estimate σ₂ and t using equation (8b) respectively. Knowing the estimates of σ and α with t₂=0, or knowing σ₂ and t₂ with α₂=0, these are then input to equation (6) to determine a simulated transport function h_s(t). The simulated tissue IRF R_s(t) can then be determined from equation (9) as below:

Once R_s(t) is determined, then the simulated concentration curve can be determined as follows: C_s(t)=(F_t/kn) AIF_t(t)®R_s(t) yF tø ^AJF^R^t-τ τ (11) Using the computer program, the user selects the initial AIF and NIF, the program will automatically derive the AIF_t(t) input to the tissue 6 based on the first model and the convolution thereof. Secondly the program will estimate tissue blood flow F_t and IRF R_e(t) and derive parameter values used to build the simulated tissue IRF R_s(t) in the second model. The program further calculates a simulated contrast curve at the tissue of interest. The seven parameters Ft, t_l5 σ_l5 α_l5 σ₂, α₂ and t₂ are optimized through a least squares method in order to fit the simulated C_s(t) to the measured tissue curve C(t). A least squares fit can be represented by a minimization process of the quantity S defined in equation (12) below:

S = ∑(C(t) - C_s(t))² (12) With the optimized seven parameters F_t, t_1} σi, α_l5 σ₂, α₂ and t_2j several quantitative perfusion indices can be determined as

Blood Flow (BF) = F_t Mean Transit Time (MTT) = t₂ + σ₂ (1+ α₂) Blood Volume (BN) = BF*MTT (13) Arterial Delay Time (DT) = ti +σι(l+αι) Arterial Dispersion Time (ADT) = σi ■Xl + a_l Tissue Dispersion Time (TDT) = σ₂ -<Jl + c ₂ Relative Arterial Dispersion (RAD) = ADT / DT Relative Tissue Dispersion (RTD) = TDT/ MTT

These indices can be determined on a pixel-by-pixel basis to produce quantitative perfusion maps respectively for further analysis and interpretation. This provides more accurate information to a clinician so that the clinician can decide on appropriate therapy for the patient on retrieving the above results or data. Thus referring again to Fig. 5, at step 132 an estimate of the transport function across the ROI is calculated by the computer program using the equation h_e(t) = R_e(t). At step 134 the program derives t₂, σ₂ and α₂ (with either α₂=0 or dt t =0) from the h_e(t) curve using the equations (7b) or (8b). At step 136 h_s(t) is derived by the program knowing the values for t₂, σ and α using the second model. At step 138 R_s(t) is derived from equation (10) by the program. At step 140 C_s(t) is determined by the program using the estimates for R_s(t), AIF_t(t), kπ and F_t. At step 142 a least squares method is used by the program to fit C_s(t) to C(t) and to optimize the parameters F_t, t_ls σ_l5 α_ls σ₂, α₂ and t₂ by minimising S in equation (12) iteratively. Finally at step 144 the program calculates values for perfusion indices such as BF, MTT and BN etc using equation (13). An artery is usually selected in the process of obtaining an arterial input function, however in the brain it is not always easy to obtain a major artery. A smaller artery in the brain may be selected instead leading to partial voluming. To compensate for partial voluming, a vein that is much larger than the artery and is usually easy to identify may be used. The user and/or computer program searches for a large vein which should have minimal partial voluming effect. A smaller artery can be selected and scaled against a vein profile. Thus, a profile of a NIF from a large draining vein is determined. The AIF is then scaled up to have the same first-pass bolus peak area as the NIF to minimise the PN effect from the AIF. The first-pass AIF and NTF profiles can be obtained by fitting them to the GNF profiles respectively to remove contrast recirculation effects. The area under the vein profile should be the same as the arterial profile. However, this approach of using a NIF_a(t) to correct for partial volume effects of AIF_a(t) is not applicable outside the brain as the contrast agent does not always remain within the vascular system during transit through the body. Usually in the body a large artery without partial voluming can be found on the imaging slices. Thus in Fig. 4 the AIF profile 80 of the original artery selected is shown, which is much smaller than the expected profile due to partial voluming. Therefore a vein is selected and it has the NIF profile 84. Due to recirculation effects, each profile shows a local maximum 82 (on the AIF curve) and 86 (on the NIF curve). A GNF is fitted by the computer program to the NIF to obtain an estimate of the total area (BN) under the fitted NIF curve whilst eliminating the local maximum 86 and following contour 87. Then the GNF is applied by the computer program to the selected AIF to eliminate the local maximum 82 and extend the profile along contour 89. The program then uses this estimate to scale up the original AIF 80 to AIF 88 to obtain an estimate of the concentration of contrast agent from the scaled up AIF 88. This approach overcomes the problem when the NIF is represented partly with a missing end portion in the data due to a shorter scanning time in order to keep to a minimum the amount of time a patient has to spent being scanned by a particular scanner. It is to be noted that the use of a GNF in each of the first and second models is for relative ease of calculations. Should the full GNF be used in both the first and second models, there will be seven adjustable parameters that need to be optimized by the least square fitting process. The computer program may provide various options to allow the user to fix certain parameters such as αi=0 and t₂=0 (or α^O and α₂=0) throughout the least square fitting process, in which only five parameters F_t, σ_2; α₂₎ σi and ti (or F_t, σ _j t_2; σi and ti) would then need to be optimized. The computer program further allows the user to fix the relative arterial dispersion βi thus σi can be calculated dependent on ti. A fixed value of βi=12%o can be chosen based on previously published results. Alternately, if one can measure AIF_t(t) by identifying a small artery showing a delay relative to AIF_a(t), optimized σi and ti values can be determined by fitting the simulated AIF_t(t) from equations (3) and (4) to the measured AIF_t(t). Then a relative dispersion βi value can be determined and applied to all other pixels of the same subjects assuming a constant relative dispersion. Thus there will be only four parameters Ft, σ _; α₂ and ti (or F_t, σ_2; t₂ and ti) that need to be optimized for increased robustness of the fitting process. Furthermore, one may apply the above approach to various subjects with vascular abnormalities accompanied by delay and dispersion, such as in acute stroke and stenosis. Once a consistent relative dispersion value βi is determined from all the representing cases, the vascular transport function in equation (3) can be described by a single variable ti together with a constant βi. A two-step method can be implemented to account for delay and dispersion. At first, an initial IRF Ro(t) can be derived by deconvolution of AIF_a(t) from C(t) using the model-free SND method. The delay time ti can be determined by the maximum position of Ro(t), i.e. Ro_{m x} ≡ Ro(t=tι). The ti value determined this way is less sensitive to curve noise because the deconvolution involves all data points of the time curve. In the second step, the AIF_t(t) feeding the ROI can be derived from equation (3) with ti and the constant β_1; which determine σi. Then value of F_t and corrected IRF R_e(t) can be obtained by deconvolution of the model derived AIF_t(t) from C(t) using the SND method. Perfusion indices can be determined from the calculated R_e(t) curve as MTT= JO R_e(τ)dτ , BF=F_t and

BN=BF*MTT. This approach can be implemented via a computer program for fast processing of perfusion maps by accounting for delay and dispersion without a time- consuming least-square-fitting process. Alternatively, as the transport function h(t) is simply a probability distribution function of the transit times, it is possible to use other functions such as a modified Gaussian function in equation (14) below to substitute equation (1) hence to describe h_a(t) and h_s(t) respectively.

1 _e- (t-t₀ )² I 2 σ ² (t > 0) A A( = (14) 0 (t < 0)

= the error

function. Using the Gaussian function to substitute the first and second models in equations (3) and (8a) respectively, there are five parameters (F_t, σ_2; t₂₎

and ti) that need to be optimized through the fitting process. Furthermore, the two models are not limited in scope to use in major vessel disease associated with the head of a patient, such as acute stroke or carotid artery stenosis. The models can be used in any intra-vascular application and therefore can apply to different parts of a patient's body, such as the cortex of the kidneys, lungs or spleen. The models can be further extended to other cases where contrast may not totally remain intravascular but leak into the tissue, such as in a tumour. For a tissue ROI with a mean transit time of τ, the tissue IRF can be described by the adiabatic approximation to the tissue homogeneity model as

where the first term is the intravascular component and the second term is the leakage component. E is the extraction fraction of the tracer in the blood stream that leaks out of the vessel into tissue, and the tracer clearance rate constant k=E*F_t /N_e is a rate constant at which the leaked contrast agent diffuses back into the blood stream and leaves the tissue, N_e is the volume fraction of the extravascular and extracellular space (EES) in the tissue. Normally there is perfusion heterogeneity associated with a distribution of transit time τ of blood vessels in tissue. Such a distribution can be described by a probability density function h_s(τ) such that the average tissue IRF involving leakage becomes R_s (t) = l - ^h_s (r)dr + Ee-^kt^h_s (τ)e^kτdτ (16) where h_s(τ) can be described by the GNF model of equation (1) or by a Gaussian distribution function of equation (14). The above described method for intravascular perfusion can be extended for perfusion measurements in a tumour by substituting equation (10) with (16) for the simulated C_s(t) in equation (11). With two additional parameters, E and N_e (or k), the method described above can be used to derive parameters for measuring both blood perfusion and permeability related indices including F_t, E and N_e. The parameters E and

N_e have a value between zero and one. The program selects certain starting values of E and N_e such as E = 0.2 and N_e = 0.4, then further optimizes E and N_e together with all other adjustable parameters (e.g F_t, σ_2j αi_, σi, ti, α _; and t₂) using the least squares method of equation (12). The permeability surface area product can then be determined by RS = -F_t ln( 1 - E) , where PS = E*F_t when E « 1. It will be appreciated by persons skilled in the art that numerous variations and/or modifications may be made to the invention as shown in the specific embodiments without departing from the spirit or scope of the invention as broadly described. The present embodiments are, therefore, to be considered in all respects as illustrative and not restrictive.

Claims

CLAIMS:

1. A method of deriving blood perfusion indices for a region of interest (ROI) of a subject, the method comprising the steps of: administering a contrast agent to the subject during a dynamic imaging scan: converting signal intensity data from raw images of the scan into contrast agent concentration data; deriving parameters from the contrast agent concentration data using at least one transport function that accounts for delay and dispersion of the contrast agent; and calculating the blood perfusion indices from the derived parameters.

2. A method according to claim 1 wherein the transport function represents a probability distribution function of transit times of the contrast agent through the subject.

3. A method according to claim 2 further comprising the step of using a first model to represent an arterial transport function h_a(t) through a vessel leading to the ROI.

4. A method according to claim 3 further comprising using a second model to represent a tissue transport function h_s(t) through the ROI.

5. A method according to claim 4 further comprising the step of selecting an arterial input function AIF_a(t) in the vessel leading to the ROI by searching pixels taken of the contrast agent concentration data.

6. A method according to claim 5 further comprising the step of measuring the contrast agent concentration C(t) remaining in the ROI.

7. A method according to claim 6 further comprising the step of representing h_a(t) using a gamma-variate function (GNF) in the first model such that:

÷-e- 'i⁾ -(f-., ) / σ, (t ≥ ^A. ω o (t < t ) where A = σ_x " ^ + α.) , T (a ) ≡ x e^{~ x} dx is the Gamma function, ti is J o the time taken for the contrast agent to move from the initial measurement of AIF_a(t) to a vessel at the entry to the ROI, σi and αi are related to the mean transit time and dispersion of h_a(t).

8. A method according to claim 7 further comprising the step of estimating h_a(t) after deriving values for parameters ti and σi and setting a_x =0 using the equation: -(. -!, ) t cr, (t ≥ ) 0 = o (t < )

9. A method according to claim 8 further comprising the step of determining an estimate for the arterial input function AIF_t(t) of the vessel at the entry to the ROI using the equation: AIFfi) (τ)h_a (t - τ)dτ

where ® is the convolution operator.

10. A method according to claim 9 further comprising the step of determining an estimate of blood flow F_t and an estimate of the tissue IRF R_e(t) from the deconvolution of: C(t) = (Ft/tø AWtf) ®R_e(t)

where kH=(l-H_a)/(l-H_t) is a correction constant taking into account different values of arterial hematocrit H_a and tissue hematocrit H_t since the contrast agent remains in the extracellular fraction of blood (plasma).

11. A method according to claim 10 further comprising the step of determining an estimate for the tissue transport function h_e(t) from the estimated R_e(t) using the equation: h_e(t) = - Re(t) dt

12. A method according to claim 11 further comprising the step of determining a rise time and a mean transit time of h_e(t) in order to determine parameters α₂ and σ₂ by assuming t₂=0, where t₂, α₂ and σ₂ are parameters related to the mean transit time and dispersion of h_e(t).

13. A method according to claim 11 further comprising the step of determining a peak height and a mean transit time of h_e(t) in order to determine parameters σ₂ and t₂ by assuming α₂=0, where t, α₂ and σ₂ are parameters relating to mean transit time and dispersion of h_e(t).

14. A method according to claim 12 or claim 13 further comprising the step of representing a simulated transport function h_s(t) using a GNF in the second model such that:

(t-t₂y^>e -(.-.₂)/σ₂ (t≥ t₂)

*,(') = 0 (t< t₂)

where A₂ =

+ o₂), t₂, σ₂ and α₂ are parameters related to the mean transit time and dispersion of h_s(t) through the ROI.

15. A method according to claim 14 further comprising the step of estimating h_s(t) using the derived values for parameters α₂ and σ₂ by setting t₂=0 using the equation: h_s(t) = — t^ae ^tlσ (t>0)

16. A method according to claim 14 further comprising the step of estimating h_s(t) using the derived values for parameters σ₂ and t₂ by setting α₂=0 using the equation: __eAt- )/σ₂ (_{t ≥} t₂) σ ■ Kit) (t < t₂)

17. A method according to claim 15 or claim 16 further comprising the step of determining a simulated tissue IRF R_s(t) using the equation:

18. A method according to claim 17 further comprising the step of determining a simulated contrast agent concentration C_s(t) using the equation:

19. A method according to claim 18 further comprising the step of fitting the simulated C_s(t) to C(t) using a least squares method according to: S = ∑(C(t) - C,(t))²

20. A method according to claim 19 further comprising the step of optimising the parameters Ft, ti.σi, αi, σ₂, α₂ and t₂by minimizing S iteratively.

21. A method according to claim 20 further comprising the step of reducing the number of adjustable parameters by fixing

leading to five adjustable parameters.

22. A method according to claim 20 or claim 21 comprising the step of further reducing the number of adjustable parameters by fixing a relative dispersion,

of h_a(t) resulting in σi dependent on t_lι and therefore leading to four adjustable parameters.

23. A method according to claim 22 further comprising the step of calculating quantitative blood perfusion indices from the optimized parameters of Ft, ti, σi, \, σ₂, α₂ and t₂.

24. A method according to claim 23 wherein the perfusion indices include any one or more of blood flow, blood volume, mean transit time, arterial delay time, arterial dispersion time or relative arterial dispersion, tissue dispersion time or relative tissue dispersion.

25. A method according to claim 24 further comprising the step of repeating each previous step, apart from the step of selecting the AIF, on a pixel-by-pixel basis to produce quantitative maps of the perfusion indices for further analysis and presentation.

26. A method according to any one of claims 1 to 25 wherein the ROI is a tissue.

27. A method according to any one of claims 1 to 25 wherein the ROI is a pixel or a plurality of pixels in a tissue.

28. A method according to any one of claims 1 to 27 wherein the scan is any one of CT, MRI or NM.

29. A method according to any one of claims 3 to 24 wherein the vessel is an artery.

30. A method according to claim 29 further comprising determining a venous input function VIF_a(t) from a draining vein to estimate an AIF_a(t) where a selected artery has partial voluming, the vein being larger than the artery.

31. A method according to claim 30 further comprising the step of detemiining the profile of NIF_a(t) from the draining vein.

32. A method according to claim 31 further comprising the step of scaling AIF_a(t) to have the same first-pass bolus peak area as the NIF_a(t) to minimize partial voluming effect from the AIF_a(t).

33. A method according to claim 32 wherein the first-pass bolus peak areas of the AIF_a(t) and NIF_a(t) profiles are obtained by fitting the profiles to gamma-variate function (GNF) profiles respectively to remove contrast recirculation effects.

34. A method according to any one of claims 17 to 33 further comprising the step of determining a simulated tissue IRF R_s(t) in the case that the contrast agent does not always remain in the vascular system, such as in a tumour in the subject in order to determine blood perfusion indices and permeability indices using:

where (t) =

E is the extraction fraction of the tracer in the blood stream that leaks out of the vessel into tissue, and the tracer clearance rate constant k=E*F_t /N_e is a rate constant at which the leaked contrast agent diffuses back into the blood stream and leaves the tissue, V_e is volume fraction of the extravascular and extracellular space (EES).

35. A method according to claim 33 wherein a permeability surface area product PS is determined by PS = -F_tln(l - E) .

36. Computer program means for deriving blood perfusion indices for a region of interest (ROI) of a subject by directing a processor to carry out any of the method steps according to any one of claims 1 to 35 apart from the step of administering a contrast agent to the subject during a dynamic imaging scan.

37. Computer program means according to claim 36 further directing the processor to retrieve raw image data from the dynamic imaging scan of the subject after a contrast agent is administered to the subject.

38. A system of deriving blood perfusion indices for a region of interest (ROI) of a subject, the system comprising: scanning means for providing a dynamic image scan of the subject during which a contrast agent is administered to the subject; processor means linked to the scanning means for retrieving raw image data from the scan; the processor means further: converting signal intensity data included in the retrieved raw image data into contrast agent concentration data; deriving parameters from the contrast agent concentration data using at least one transport function that accounts for delay and dispersion of the contrast agent; and calculating the blood perfusion indices from the derived parameters.

39. A system according to claim 38 wherein the transport function represents a probability distribution function of transit times of the contrast agent through the subject.

40. A system according to claim 39 wherein a first model is used to represent an arterial transport function h_a(t) through a vessel leading to the ROI.

41. A system according to claim 40 wherein a second model is used to represent a tissue transport function h_s(t) through the ROI.

42. A system according to claim 41 wherein the processor means selects an arterial input function AIF_a(t) in the vessel leading to the ROI by searching pixels taken of the contrast agent concentration data.

43. A system according to claim 42 wherein the processor means measures the contrast agent concentration C(f) remaining in the ROI.