Extracting parametric images from dynamic contrast-enhanced MRI studies of the brain using factor analysis

Abstract

Factor analysis of dynamic studies (FADS) is a technique that allows structures with different temporal characteristics to be extracted from dynamic contrast enhanced studies without making any a priori assumptions about physiology. These dynamic structures may correspond to different tissue types or different organs or they may simply be a useful way of characterising the data. This paper describes a method of automatically extracting factor images and curves from contrast enhanced MRI studies of the brain. This method has been applied to 107 studies carried out on patients with acute stroke. The results show that FADS is able to extract factor curves correlated to arterial and venous signal intensity curves and that the corresponding factor images allow a distinction to be made between areas of the brain with normal and abnormal perfusion. The method is robust and can be applied routinely to dynamic studies of the brain. The constraints described are sufficiently general to be applicable to other dynamic MRI contrast enhanced studies where an increase in contrast concentration produces an increase in signal intensity.

Introduction

Stroke is a major cause of death and long term disability. The ability to measure regional brain blood flow in vivo is vital in order to understand the pathophysiology of stroke and to assess the effect of new therapies such as neuroprotective agents and thrombolytic drugs on the disease process. It is also useful to assess brain perfusion in patients with other cerebral vascular diseases such as carotid artery disease and in determining the vascularity of brain tumours.

Both quantitative and qualitative measurements of cerebral perfusion can be made by injecting a bolus of contrast agent intravenously and then imaging the changes in signal intensity as the bolus of contrast agent makes its first-pass through the brain. Many different imaging modalities have been used to carry out such studies in the past (Axel, 1980, Lassen, 1982, Hatazawa et al., 1997)^.and recent advances in MR imaging techniques, in particular the development of echo planar imaging (EPI) and rapid gradient echo (GE) sequences, have lead to considerable interest in the use of MRI to obtain non-invasive measurements of brain blood flow.

In dynamic contrast enhanced MRI studies the passage of a bolus of paramagnetic contrast agent such as gadolinium diethylene triamine pentaacetic acid (Gd-DTPA) is tracked through the brain. The contrast agent decreases both the T1 and T2 relaxation rate of tissue. An decrease in T1 typically causes an increase in signal intensity whereas an decrease in T2 causes a decrease in signal intensity – the net effect on the MR signal intensity of these opposing effects will depend on the type of imaging sequence used.

T2* weighted gradient echo images are the most commonly used to assess cerebral perfusion (Rempp et al., 1994, Sorensen et al., 1996). As the contrast agent passes through the capillary network of the brain it produces a reduction in signal intensity. This technique has two main advantages. The changes in signal intensity due the passage of a bolus of contrast agent through the capillaries are large, which reduces the effect of noise on the subsequent numerical analysis. It also allows the acquisition of multiple slices through the brain, allowing the extent of any abnormalities to be assessed. There are however several disadvantages of this technique. If the blood brain barrier is no longer intact then the contrast agent leaks into the extravascular space and T1 effects begin to compete with the T2* effects. This makes the T2* method unsuitable for assessing cerebral blood flow in the presence of tumours. The images may be distorted in regions where bone and air are in close proximity due to susceptibility artefacts and only MRI scanners equipped with fast gradient systems suitable for EPI sequences are able to provide adequate spatial and temporal resolution. The T2* signal intensity curves are generally transformed to give curves where the signal intensity is proportional to concentration using the expression ΔR2*=−ln[S_t/S₀]/TE where S₀ is the pre-contrast signal intensity and S_t is the post-contrast signal intensity at time t. TE is the echo time which will depend on the image sequence and ΔR2* is equal to 1/ΔT2*.

The T1 relaxation effect may also be used to measure the concentration of the contrast agent in the tissue (Dean et al., 1992, Hacklander et al., 1996a, Hacklander et al., 1996b, Moody et al., 2000). The main advantage of the T1 imaging techniques is that quantitative measurements of cerebral blood flow can be made from the image data. If very low doses of contrast agent are used then the increase in signal intensity produced by the relaxation effect is directly proportional to contrast concentration (Canet et al., 1995). The use of a small dose means that the volume of contrast agent injected is 10–20 times smaller than that required for the T2* technique allowing a bolus injection to be administered without the use of a power injector. Other investigations requiring the use of additional doses of contrast agent such as repeated first pass perfusion studies or MR angiography can also be carried out during the same session without exceeding the maximum permissible contrast dose. For these reasons we have chosen to use T1 weighted imaging to assess perfusion. The main disadvantage of the T1 weighted technique is that the changes in signal intensity due to the contrast agent are small, making the signal intensity curves more noisy.

The techniques used to extract information from the dynamic MRI studies are similar regardless of which imaging sequence has been used to acquire the data. A simple parametric image showing time to peak signal intensity is the most commonly used since it is quick and simple to calculate (Sorensen et al., 1996, Sunshine et al., 1999). The relative cerebral blood volume (rCBV) which can be approximated by calculating the area under the signal intensity curve is also frequently reported (Dean et al., 1992, Hacklander et al., 1996a, Hacklander et al., 1996b, Rordorf et al., 1998). From the indicator dilution theory (Axel, 1980) the relationship between the CBV, the mean transit time (MTT) and cerebral blood flow (CBF) is given by the expression MTT=CBV/CBF. A parametric image of MTT can be estimated by fitting a gamma variate to each individual signal intensity (SI) curve. The CBF image is then generated from the CBV and MTT images. This is possible with T2* data (Hagen et al., 1999) which has a high signal to noise ratio but is not practical with T1 weighted images which are more noisy.

An alternative approach to estimating the CBF is to use the gradient of the signal intensity curves. This method was originally developed for use with nuclear medicine data (Peters et al., 1987) and has been used successfully with T1 weighted MRI data (Moody et al., 2000). The technique makes the assumption that no output from the tissue occurs until the whole bolus has entered the region of interest therefore it is not suitable for use with T2* data since the bolus injection is much broader than is the case with the T1 weighted images.

The indicator dilution theory imposes a mathematical model on the data and this may not be valid in certain cases. The technique cannot be used where the blood brain barrier is damaged since the leakage of contrast agent into the interstitial space invalidates the assumptions made. The theory also makes assumptions about vascular topology (Weisskoff et al., 1993) and these may be incorrect in patients with complex vascular anomalies such as Moya Moya disease.

Attempts have been made to extract information from dynamic contrast enhanced studies without imposing any a priori model to the data and similarity mapping is one such approach. This involves defining a reference region over an area of interest, e.g. a tumour, and then correlating the pixel signal intensity curves with this reference curve to form a similarity map (Rogowska et al., 1995). This method has been applied to CT images of rabbits with focal cerebral ischemia and was able to identify small differences in the temporal dynamics around the infarct (Lo et al., 1996). The main disadvantage of this approach is that it requires the manual selection of a reference region.

In this paper we investigate the use of factor analysis to extract functional images. Factor analysis of dynamic studies (FADS) is a technique that allows structures with different temporal characteristics to be identified automatically without making any a priori assumptions about physiology. These dynamic structures may correspond to different tissue types or different organs or they may simply be a useful way of characterising the data. FADS was originally developed for use with nuclear medicine data (Barber, 1980, Di Paola et al., 1982) and has more recently been applied to MRI (Bonnerot et al., 1992, Zagdanski et al., 1994).

In the factor model it is assumed that the correlation between individual pixel concentration-time curves can be explained in terms of a small number of underlying factor curves. These factors will often describe some physical property of the system being observed. Alternatively, they may be theoretical constructs which have no physiological significance but which simplify the task of interpreting the data. The dynamic study can be represented by the equation

$$\text{d}\text{x, t}\text{=}\text{∑}\text{k=1}\text{M}\text{a}{}_{\mathrm{k}}\text{x}\text{f}{}_{\mathrm{k}}\text{t}\text{+e,}$$

where d(x, t) is the signal intensity in pixel x at time t, f_k is the kth factor curve, a_k is the kth factor image and e is the residual error. M should be equal to the number of kinetic compartments present in the data. If sufficient factors have been identified then the variance represented by the error term e is due to random noise only.

Since there are an infinite number of possible solutions to Eq. (1) it is necessary to apply constraints in order to obtain a unique solution. Principal components analysis (PCA) uses a statistical constraint to obtain a unique set of orthogonal factors with no a priori assumptions being made about the data. The dynamic study can be represented by the (T×N) matrix D, where T is the number of frames in the study and N is the number of pixels in each frame. Each row of the matrix represents an image in the sequence and each column represents a pixel signal intensity curve. The principal component (PC) curves are obtained by extracting the eigenvectors of the covariance matrix C (given by DD^T) in decreasing order of importance, i.e. λ_i>λ_i+1 where λ_i is the eigenvalue corresponding to the ith eigenvector $\text{u}{}_{\mathrm{i}}$. The first PC accounts for most of the information in the study, with subsequent PCs containing progressively less. We have assumed that there exists a subset of M PCs which account for all of the useful information while the remaining (T−M) PCs represent random noise. The dynamic study can therefore be represented by

$$\text{D}\text{=}\text{∑}\text{i=1}\text{M}\text{u}{}_{\mathrm{i}}\text{v}{}^{\mathrm{<mtext>T</mtext>}}{}_{\mathrm{i}}\text{λ}{}_{\mathrm{i}}\text{+}\text{∑}\text{i=M+1}\text{T}\text{u}{}_{\mathrm{i}}\text{v}{}^{\mathrm{<mtext>T</mtext>}}{}_{\mathrm{i}}\text{λ}{}_{\mathrm{i}}\text{,}$$

where $\text{u}{}_{\mathrm{i}}$ is the (T element) ith principal component curve and $\text{v}{}_{\mathrm{i}}$ is the corresponding (N element) vector of coefficients. By definition the PC curves are mutually orthogonal which means that they do not represent true physiological structures. The problem of identifying a set of physiologically meaningful oblique factors can be simplified by extracting the first M principal components in order to reduce the dimensionality of the data set, and then rotating these components under the control of certain constraints. Substituting for D from Eq. (2) in Eq. (1) and ignoring the residual matrix E gives

$$\text{FA}\text{=}\text{U*}\text{R*}\text{V*}\text{,}$$

where U* is the (T×M) orthonormal column matrix containing the first M PC curves, V* is the (M×N) row matrix of PC images, R* is the diagonal matrix of eigenvalues, F is the (T×M) column matrix containing the M factor curves and A is the (M×N) row matrix of coefficients or factor images. From Eq. (3) F and A can be represented by

$$\text{F}\text{=}\text{U*}\text{R*}\text{T}\text{,}$$

$$\text{A}\text{=}\text{T}{}^{\mathrm{-1}}\text{V}\text{,}$$

where T is an (M×M) rotation matrix and $\text{TT}{}^{\mathrm{-1}}\text{=}\text{I}$. The aim of any FADS technique is to find a matrix T which rotates the principal components in such a way that the final factor solution satisfies certain constraints. The most commonly used is the positivity constraint (Barber, 1980) which assumes that neither the factor images nor the factor curves should contain any negative values. This constraint has been applied to various types of nuclear medicine studies (Cavailloles et al., 1984, Cinotti et al., 1985) and more recently to contrast enhanced MRI studies (Bonnerot et al., 1992, Zagdanski et al., 1994) with some success. The positivity constraint is not sufficient on its own to produce a unique solution (Houston, 1984) and the use of additional constraints, both spatial (Nijran and Barber, 1985, Samal et al., 1987) and temporal (Nijran and Barber, 1987), has been investigated. The constraints we have used to extract factors from dynamic MRI studies are described in detail in Section 2.2.

The aim of this study was to establish a method of extracting factor images from dynamic contrast-enhanced MR images of the brain. In particular we wanted to determine whether FADS could be carried out routinely in the setting of acute stroke and whether it could provide useful information about cerebral haemodynamics.