WO2010045478A2

WO2010045478A2 - System and method for differential perfusion analysis of tissues with multiple vascular supplies

Info

Publication number: WO2010045478A2
Application number: PCT/US2009/060875
Authority: WO
Inventors: Michael H. Rosenthal; Francine L. Jacobson; Hiroto Hatabu
Original assignee: Brigham And Women's Hospital, Inc.
Priority date: 2008-10-15
Filing date: 2009-10-15
Publication date: 2010-04-22
Also published as: WO2010045478A3

Abstract

A method for calculating the contribution of multiple vascular blood flow sources to tissues in a region of interest undergoing perfusion analysis establishes a model time-enhancement reference curve comprised of the weighted combination of reference time-enhancement curves calculated for each blood flow source. The weighting coefficients are calculated for each region of interest being analyzed using the model time-enhancement curve and a time-enhancement curve of the tissue at that location.

Description

SYSTEM AND METHOD FOR DIFFERENTIAL PERFUSION ANALYSIS OF TISSUES WITH MULTIPLE VASCULAR SUPPLIES

CROSS-REFERENCE TO RELATED APPLICATIONS

[0001] This application claims benefit of U.S. Provisional patent application Serial No. 61/105,745 filed October 15, 2008 and entitled " System and Method for Differential Perfusion Analysis of Tissues With Multiple Vascular Supplies."

BACKGROUND OF THE INVENTION

[0002] Perfusion imaging is defined as any imaging study that produces a series of images over time that depict blood perfusion in an imaged tissue. From these series of images time-enhancement curves may be generated and hemodynamic parameters such as blood flow, blood volume, and mean transit time may be computed. As pathologic conditions often alter these properties, perfusion analysis provides improved detection, diagnosis, and assessment of therapeutic response for such conditions.

[0003] Typically, perfusion imaging is performed by administering a contrast agent to a subject, acquiring a series of MR or CT images as the contrast agent perfuses into the tissues of interest, and then computing a set of hemodynamic parameters describing blood flow, blood volume, and mean transit time. Perfusion imaging may be performed using CT systems with an iodinated contrast agent, for example, adrenal wash-in, wash-out studies, or using MRI systems with or without an intravenous gadolinium contrast medium. These methods are currently in clinical use for functional evaluation of the heart and tumor evaluation in the brain, but they are not widely used for the remainder of the body.

[0004] For example, hemodynamically weighted MR perfusion images of cerebral blood flow (CBF) may be acquired and used in combination with diffusion-weighted (DWI) MR images of the brain to delineate regions of viable brain parenchyma that are at risk of further infarction following a stroke. The DWI MR image shows ischemic regions where brain cells have died due to the stroke, and the CBF image shows regions with reduced blood flow that indicate at risk tissue. The size of the "ischemic penumbra" surrounding ischemic tissues is a critical factor in evaluating treatment options.

[0005] As described, for example, by K.A. Kemp, et al. "Quantification of Regional Cerebral Blood Flow and Volume with Dynamic Susceptibility Contrast Enhanced MR Imaging" Radiology 1994; 193:637-641, it is possible to assess regional cerebral hemodynamics by analyzing MR signal intensity changes after the first pass of the paramagnetic contrast medium. While passing through the capillary network, a short bolus of contrast medium produces local magnetic field inhomogeneities that lead to a reduction in the transverse magnetization relaxation time T2* of the bulk tissue. This susceptibility effect can be recorded by a series of rapid T2*-weighted gradient- echo images that reveal how the MR signal changes during the first pass of the contrast medium. The resulting MR signal-intensity-versus-time curves can be converted into contrast-agent-concentration-time curves, that is, time-enhancement curves. By using the indicator dilution theory, two important hemodynamic parameters can be determined from these curves. First, the cerebral blood flow (CBF), known as tissue perfusion can be determined. Second, the cerebral blood volume (CBV) can be determined. However, if absolute quantification of the CBV and CBF measurements are to be achieved, the concentration of contrast agent in the arterial blood pool - the "arterial input function" (AIF) - must be known. [0006] Typical methods used to measure the AIF require an additional step in which the operator manually selects a region of interest (ROI), based on anatomic information, that depicts an artery. The concentration-time curve from all voxels included in the ROI are then used to calculate the AIF as described, for example, by L. Ostergaard et al. "High Resolution Measurement of Cerebral Blood Flow Using Intravascular Tracer Bolus Passages. Part 1: Mathematical Approach and Statistical Analysis", Magnetic Resonance In Medicine. 36:715-725 (1996) and B. R. Rosen et al. "Perfusion Imaging With NMR Contrast Agents," Magnetic Resonance In Medicine, 14, 249-265 (1990). More recently, automatic methods for determining AIF have been used such as the method disclosed in U.S. Patent No. 6,546,275. [0007] These traditional perfusion analysis methods assume a singular vascular compartment with an accompanying tissue compartment that allows contrast material to flow between compartments according to standard kinetic diffusion equations. However, this assumption excludes the adequate assessment of tissues with dual vascular inputs such as the lungs, which receive vascular contributions from the pulmonary and bronchial arteries, and the liver, which receives vascular contributions from the portal vein and the hepatic artery. In particular, these traditional methods are unable to differentiate vascular contributions from multiple sources and cannot properly assess non-physiologic conditions, for example, cancer, that may alter the relative balance of the vascular contributions from each source.

[0008] Previous efforts to analyze dual-input tissues have been based solely upon the assumption of exponential kinetic models to describe the flow of contrast within a region of interest and the morphological analysis of time-enhancement curves. The fact that most voxels acquired during a perfusion imaging study contain more than one type of tissue means that this assumption is not valid and prohibits effective analysis using such methods.

[0009] It would therefore be desirable to develop a method of perfusion analysis, that is not constrained to an exponential kinetic decay model and that enables the estimation of fractional blood volumes.

SUMMARY OF THE INVENTION

[0010] [Insert Summary]

[0011] The present invention uses a class of statistical methods to estimate the relative contributions of multi-input vascular supplies, for example, the fractional blood volumes contributed by the bronchial and pulmonary arteries in a region of interest (ROI) in the lung and those contributed by the portal vein and the hepatic artery in tissues in the liver.

[0012] Various other features of the present invention will be made apparent from the following detailed description and the drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

[0013] Fig. 1 is a schematic diagram of a magnetic resonance imaging system used to perform a preferred embodiment of the invention; and

[0014] Fig. 2 is a flowchart setting forth the steps used to practice the present invention with a medical imaging system such as that in Fig. 1.

DETAILED DESCRIPTION OF THE INVENTION

[0015] As indicated above, perfusion imaging is commonly performed with either an x-ray CT system or an MRI system. In either case, a series of images are acquired of the tissues of interest which depict enhancement due to the inflow of arterial blood. In the preferred embodiment now to be described, these perfusion images are acquired with an MRI system. The processing of the images in accordance with the present invention may be performed on the MRI system, or the images may be off-loaded to a separate work station for processing.

[0001] Referring to Fig. 1, the MRI system includes a workstation 10 having a display 12 and a keyboard 14. The workstation 10 includes a processor 16 that is a commercially available programmable machine running a commercially available operating system. The workstation 10 provides the operator interface that enables scan prescriptions to be entered into the MRI system. The workstation 10 is coupled to four servers including a pulse sequence server 18, a data acquisition server 20, a data processing server 22, and a data store server 23. The workstation 10 and each server 18, 20, 22 and 23 are connected to communicate with each other. [0002] The pulse sequence server 18 functions in response to instructions downloaded from the workstation 10 to operate a gradient system 24 and an RF system 26. Gradient waveforms necessary to perform the prescribed scan are produced and applied to the gradient system 24 that excites gradient coils in an assembly 28 to produce the magnetic field gradients G_x, G_y and G_z used for position encoding MR signals. The gradient coil assembly 28 forms part of a magnet assembly 30 that includes a polarizing magnet 32 and a whole-body RF coil 34. [0003] RF excitation waveforms are applied to the RF coil 34 by the RF system 26 to perform the prescribed magnetic resonance pulse sequence. Responsive MR signals detected by the RF coil 34 or a separate local coil (not shown in Fig. 1) are received by the RF system 26, amplified, demodulated, filtered, and digitized under direction of commands produced by the pulse sequence server 18. The RF system 26 includes an RF transmitter for producing a wide variety of RF pulses used in MR pulse sequences. The RF transmitter is responsive to the scan prescription and direction from the pulse sequence server 18 to produce RF pulses of the desired frequency, phase and pulse amplitude waveform. The generated RF pulses may be applied to the whole body RF coil 34 or to one or more local coils or coil arrays (not shown in Fig. 1).

[0004] The RF system 26 also includes one or more RF receiver channels.

Each RF receiver channel includes an RF amplifier that amplifies the MR signal received by the coil to which it is connected and a detector that detects and digitizes the I and Q quadrature components of the received MR signal. The magnitude of the received MR signal may thus be determined at any sampled point by the square root of the sum of the squares of the I and Q components:

M = V^l2 +Q² > and the phase of the received MR signal may also be determined: φ = tan^"1 Q/r.

[0005] The pulse sequence server 18 also optionally receives patient data from a physiological acquisition controller 36. The controller 36 receives signals from a number of different sensors connected to the patient, such as ECG signals from electrodes or respiratory signals from a bellows. Such signals are typically used by the pulse sequence server 18 to synchronize, or "gate", the performance of the scan with the subject's respiration or heart beat.

[0006] The pulse sequence server 18 also connects to a scan room interface circuit 38 that receives signals from various sensors associated with the condition of the patient and the magnet system. It is also through the scan room interface circuit 38 that a patient positioning system 40 receives commands to move the patient to desired positions during the scan.

[0007] The digitized MR signal samples produced by the RF system 26 are received by the data acquisition server 20. The data acquisition server 20 operates in response to instructions downloaded from the workstation 10 to receive the realtime MR data and provide buffer storage such that no data is lost by data overrun. In some scans the data acquisition server 20 does little more than pass the acquired MR data to the data processor server 22. However, in scans that require information derived from acquired MR data to control the further performance of the scan, the data acquisition server 20 is programmed to produce such information and convey it to the pulse sequence server 18. For example, during prescans, MR data is acquired and used to calibrate the pulse sequence performed by the pulse sequence server 18. Also, navigator signals may be acquired during a scan and used to adjust RF or gradient system operating parameters or to control the view order in which k- space is sampled. And, the data acquisition server 20 may be employed to process MR signals used to detect the arrival of contrast agent in an MRA scan. In all these examples the data acquisition server 20 acquires MR data and processes it in realtime to produce information that is used to control the scan.

[0008] The data processing server 22 receives MR data from the data acquisition server 20 and processes it in accordance with instructions downloaded from the workstation 10. Such processing may include, for example, Fourier transformation of raw k-space MR data to produce two or three-dimensional images, the application of filters to a reconstructed image, the performance of a backprojection image reconstruction of acquired MR data; the calculation of functional MR images, the calculation of motion or flow images, and the like. [0009] Images reconstructed by the data processing server 22 are conveyed back to the workstation 10 where they are stored. Real-time images are stored in a data base memory cache (not shown) from which they may be output to operator display 12 or a display 42 that is located near the magnet assembly 30 for use by attending physicians. Batch mode images or selected real time images are stored in a host database on disc storage 44. When such images have been reconstructed and transferred to storage, the data processing server 22 notifies the data store server 23 on the workstation 10. The workstation 10 may be used by an operator to perform the perfusion analysis on the acquired images as described in more detail below with respect to Fig. 2.

[0016] Referring to Fig. 2, as indicated at process block 300, perfusion data from a tissue of interest having dual vascular supplies is acquired using a medical imaging system. For example, the MRI system of Fig. 1 can be used to obtain perfusion data from a subject that has been administered an intravenous contrast-enhancing agent. This is done by acquiring a series of images of the tissues of interest at frequent time intervals, for example, every 5 seconds. An observed time-enhancement curve for each subregion in the region of interest (ROI) is then determined from the signal magnitudes at each subregions over the image series. It is contemplated that the subregion may cover the entire or substantially of the ROI, may include a single voxel in the ROI, or may include a collection of voxels in the ROI. [0017] A T₂-weighted imaging pulse sequence is used to operate the MRI system so that the acquired MR signals are sensitive to changes in contrast agent concentration in the inflowing arterial blood. The characteristics of time- enhancement curves vary between tissues depending on the vascular content and flow of the tissues. For example, arteries demonstrate rapid contrast inflow and clearance, while connective tissue exhibits minimal levels of contrast enhancement and flow. It should be noted that a CT system or other imaging modality could also be used to acquire a similar series of images that could be analyzed using the methods of the present invention.

[0018] As indicated at process block 302, reference time-enhancement curves are produced by manually drawing ROIs over vascular system components that supply arterial blood to the tissues of interest. For example, the aorta, the pulmonary artery and vein, and the portal vein, might be identified in this manner within a time-series of perfusion images. It is contemplated that these reference time-enhancement curves can be automatically extracted by identifying subregions within time-perfusion imaging studies that cannot be decomposed into the sum of component time- enhancement curves and provide a satisfactory basis for describing the composition of subregions or voxel(s) having multiple component time-enhancement curves. A time-enhanced reference curve is thus produced for each source of blood flowing into the tissues of interest.

[0019] Referring still to Fig. 1, as indicated at process block 304 a model is then created to define a model time-enhancement curve for each region of interest in perfusion image data, S_mOd(t), as a statistical composition of the reference time- enhancement curves, S_r1(t) to Sm(t). Like subregions, ROIs can be defined as individual voxels or as arbitrary collections of voxels. When multiple voxels are included within a region of interest, then S_mocι(t) may, for example, be defined as the mean signal of all of the voxels within the region of interest at each time t. It is contemplated that a variety of summary methods can be used including median values, robust classifiers such as excluding statistically outlying voxels from the mean calculation, and the like. The model, shown below in equation 1 , includes "mixture terms," ai to a_n, that weight the contributions of the calculated reference time-enhancement curves and error terms, εt, for each time point t.

[0020] Smod(t)=ai * Sri(t) + a₂ * Sr2(t) + ... a_n * S_m(t) + £_t (1 ).

[0021] The structure of this model allows the calculated reference time- enhancement curves to define the time-enhancement model directly, whereas previous attempts constrained the observed time-enhancement curves to an exponential kinetic model. Moreover, a variety of estimator vectors, A, exist for the mixture terms, ai to a_N, and these estimators A can be determined using a number of statistical analysis techniques, indicated generally at 306, including least-squares and expectation maximization techniques.

[0022] In the preferred embodiment, an estimator vector A is determined using an iterative method outlined in process blocks 308-314. The observed reference time- enhancement curves are loaded into the model S_mOd(t) and the model is evaluated at process block 308 using initial values for the mixture terms ai to a_n to calculate error terms for each voxel in the ROI as the difference between the observed time- enhancement curve for the voxel S_Ob(t) and the estimated mixture of reference time- enhancement curves, S-ι(t) to S_n(t), weighted by the mixture terms, ai to a_n, as shown in equation 2.

[0023] B₁ =

S a₁ - S₁ Λ(I) ' (2). i=1

[0024] It is noted that appropriate initial values of the mixture terms may include setting the terms to zero.

[0025] As indicated at process block 310, the error terms are analyzed. If, at decision block 314, the error exceeds a threshold, then as indicated at process block 314, the model evaluation process is refined, generally by altering the mixture terms of the estimator vector A using an optimization method, such as the trust-region- reflective optimization variant of the nonlinear least squares method as implemented, for example, in the Matlab software package. Matlab is a registered trademark of MathWorks, Inc. of Natick, MA. The mixture terms may be constrained to be greater than or equal to zero and less than or equal to 1 at each iteration. [0026] Still referring to Fig. 1, the model is then reevaluated at process block 310 using the altered mixture terms. This iterative process is repeated until the values of the weighted sum of reference time-enhancement curves and the observed time- enhancement curve converge to an acceptable degree. That is, until the change in error terms between iterations decreases below the threshold defined at decision block 314. For example, a threshold value of 0.01 has proven to be effective in some implementations. In some configurations, no absolute error threshold is used; however, such a threshold could also be applied in some configurations. Then, the estimator vector A is stored at process block 316 and used to determine the contributions of vascular inputs in subsequent calculations as indicated at process block 318.

[0027] The estimator vector A can be calculated using any of a large set of optimization techniques. For example, one estimator for the mixture terms is the least squares estimator, A=Ia₁ a_n], such that the square of the difference between the time-enhancement curve of the tissue region, which is obtained from the acquired perfusion data, and the reference time-enhancement curves is minimized. This estimator, shown below in equation 3, can be obtained using, for example, any descent method of statistical analysis. One such method is the trust- region-reflective variation of non-linear least squares estimation. Mixture terms are again constrained to a range of, for example, [0, 1]. [0028] A = arg(min(S_u(t) - ai * Si(t) - a₂ * S₂(t) - ... a_n * S_n(t))²) (3)

[0029] Another method for calculating the estimator vector A for the mixture is the maximum likelihood estimator, which incorporates the variance of the signal from the family of voxels within each reference and tissue region of interest. A probability function for the composite unknown signal S_u(t) that is normally distributed with parameters N(mu_u, var_u) is obtained by assuming that Gaussian white noise is added to the baseline signal at each reference and observed voxel and that these variances are independent and not necessarily equal. As used herein, m_u is the weighted sum of background signals from reference regions-of-interest (ROIs) plus the unknown ROI and var_u is the sum of the variances of the unknown ROI voxels plus the weighted sums of variances of the references ROIs. Finding the zeros of the partial first derivatives of the natural log of this probability function produces the maximum likelihood estimate of the estimator of the mixture terms, A=[ai ... a_n]. [0030] Expectation maximization methods can be used to estimate the true fractional volumes of pulmonary arterial versus bronchial arterial vasculature with standard errors of, for example, less than 3% of parenchymal volume for 100 pixel ROIs with signal-to-noise-ratios (SNRs) as low as 0.20 at sampling frequencies of 0.1 Hz or with SNRs greater than or equal to one for any ROI of 25 pixels or greater area. Standard errors of 1% or less have been achieved with SNRs greater than or equal to one and ROIs greater than 81 pixels.

[0031] The present invention is not limited to the study of dual vascular input tissues like the lung and liver, as its methods may be applied to the assessment of abnormalities of tissues with a singular vascular input. It is further contemplated that the present invention has applications in benign processes such as differentiating infection from atelectasis. That is, the present invention has applications in compressed lung, and providing functional physiologic data about vascular flows, for example, cardiac output and relative lung perfusion. It should be appreciated by one of ordinary skill of the art that a large number of statistical methods and estimators utilize the methods of the present invention to determine the fractional blood components of tissues having multiple vascular sources.

[0032] The described invention allows the estimation of fractional blood volumes without constraint to an exponential kinetic decay model. That is, the invention views the time-enhancement curves as a composition of reference signals, instead of basing the time-enhancement curves on exponential kinetic models. [0033] The present invention has been described in accordance with the embodiments shown, and one of ordinary skill in the art will readily recognize that there could be variations to the embodiments, and any variations would be within the spirit and scope of the present invention. Accordingly, many modifications may be made by one of ordinary skill in the art without departing from the spirit and scope of the appended claims.

Claims

1. A method for determining the relative contributions of vascular blood flow from multiple sources to tissues within a region of interest (ROI) undergoing perfusion analysis, the steps comprising: a) acquiring a time series of images that depict the inflow of blood into the tissues; b) producing a time-enhanced curve for subregions in the ROI undergoing perfusing analysis using the acquired images; c) producing a reference time-enhanced curve for each source of vascular blood flow using the acquired images; d) producing a model time-enhanced curve for a selected subregion of the ROI undergoing perfusing analysis as a statistical combination of the reference time- enhancement curves produced in step c), wherein the contribution of each time- enhancement curve to the model time-enhancement curve is weighted accordingly to the contribution of corresponding blood flow source to blood in the selected subregion of the ROI; and e) calculating the weighting coefficients that indicate relative contributions from the multiple sources of blood flow using the model time-enhancement curve and a time-enhancement curve for the selected ROI.

2. The method as recited in claim 1 in which step e) is an iterative process in which the weighting coefficients are altered until the difference between the model time-enhancement curve and the time-enhancement curve for the selected subregion of the ROI changes less between iterations than a selected threshold.

3. The method as recited in claim 1 in which steps d) and e) are repeated for each subregion of the ROI undergoing perfusion analysis.

4. The method as recited in claim 1 in which the selected subregion is a voxel in the ROI.

5. The method as recited in claim 1 in which the selected subregion is a collection of voxels in the ROI.

6. The method of claim 1 wherein the subregion of he ROI includes substantially all of the ROI.

7. The method as recited in claim 1 in which the vascular blood flow is arterial blood flow.