Sage Journals: Discover world-class research

Abstract

The goal is to quantify the fraction of tissues that exhibit specific tracer binding in dynamic brain positron emission tomography (PET). It is achieved using a new method of dynamic image processing: clustering-initiated factor analysis (CIFA). Standard processing of such data relies on region of interest analysis and approximate models of the tracer kinetics and of tissue properties, which can degrade accuracy and reproducibility of the analysis. Clustering-initiated factor analysis allows accurate determination of the time–activity curves and spatial distributions for tissues that exhibit significant radiotracer concentration at any stage of the emission scan, including the arterial input function. We used this approach in the analysis of PET images obtained using ¹¹C-Pittsburgh Compound B in which specific binding reflects the presence of β-amyloid. The fraction of the specific binding tissues determined using our approach correlated with that computed using the Logan graphical analysis. We believe that CIFA can be an accurate and convenient tool for measuring specific binding tissue concentration and for analyzing tracer kinetics from dynamic images for a variety of PET tracers. As an illustration, we show that four-factor CIFA allows extraction of two blood curves and the corresponding distributions of arterial and venous blood from PET images even with a coarse temporal resolution.

Keywords

Alzheimer imaging Dynamic PET image-derived input function time series analysis

INTRODUCTION

The principal data processing tasks in dynamic brain positron emission tomography (PET) are determining the amount of radiotracer in the blood or arterial input function (IF) and tissue classification. Image-derived IF is important because arterial cannulation, the standard method for determining the IF, is an invasive and potentially risky procedure¹ that should be avoided for patient health and clinical efficiency considerations. Tissue classification groups voxels into tissue types based on the fact that the shape of the time–activity curve (TAC) for the same input function is determined by the interaction of the tissue and the radiotracer. In this manuscript, we describe an application of an extension of factor analysis of dynamic structures (FADS) technique to solve these two data processing tasks.

Despite a physical impossibility to differentiate between plasma concentration and metabolites or bound tracer, a variety of data-based techniques of determining the IF have been suggested in the last decades.^2–4 In these and other reports, the IF is determined from the reconstructed images with some amount of bias caused by partial volume and/or inaccurate determination of the region of interest (ROI) of the image where the blood activity is digitally sampled, usually in the carotid artery. A recent review² concludes that for most PET tracers and dynamic protocols, the IF cannot be determined accurately only from the dynamic images. Some of the limitations of image-derived IF techniques are not present in FADS, a powerful method used in both SPECT and PET for over a decade,^5–8 although still underused. The FADS circumvents the purely methodological constraint of voxel-sampling-based analysis by decoupling the temporal and the spatial components of the tracer concentration distribution. It has been demonstrated to measure the IF in cardiac imaging by estimating the concentration in both right and left ventricular blood pools.⁹ In brain imaging, the IF is more challenging because of the lack of large blood pools as seen in the heart chambers. Nevertheless, in this manuscript, we demonstrate the ability of FADS to recover the IF from dynamic brain PET data.

The quantitative task addressed in the paper is related to dynamic tissue separation. Most reported methods of quantifying tissues from dynamic images are based on kinetic modeling, e.g., the graphical analysis method¹⁰ or alternative approaches.¹¹ The FADS can serve as the tissue classification tool that does not imply the knowledge of the IF or a specific kinetic model. A state-of-the-art application of FADS methodology to brain imaging⁸ achieves successful separation of tissue types and TACs for two different radiotracers injected simultaneously. For a single tracer, a typical classification task is differentiating tissues that exhibit nonspecific binding and specific binding. Specific binding reflects binding of a PET ligand to a receptor, which is saturable and which can be blocked by another ligand with known activity at the receptor, whereas nonspecific binding is neither saturable nor blockable. A variety of methodologies have been used for this task.^10,12 These analyses generally rely on a reference tissue model,¹³ which assumes the knowledge of the tracer dynamics in a region of the brain that contains only the nonspecific binding tissue. When such a region is known, the dynamic data from the ROIs and from the reference region are analyzed using any of the standard techniques^10–12 to provide the quantitative estimate of the amount of specific binding tissue in the target ROIs. Although well-studied, the reference tissue model can yield results that are inaccurate or biased as the result of manual reference region selection; spatial inhomogeneity of the tissue concentrations; patient-to-patient or study-to-study variability of the tissue properties; and the time when the tracer concentrations in the target and reference regions reach a steady state. We seek to avoid some of these problems by using FADS, which does not rely on the named assumptions and on significant operator intervention.

One of the common tracers in dynamic brain imaging is ¹¹C-labeled Pittsburgh Compound B (¹¹C-PIB). It binds to β-amyloid, a protein commonly linked to Alzheimer's disease.^14–15 The physiological mechanisms and the diagnostic value of this link are outside the scope of the presented work; we are interested in the very fact that there is persistent need for quantitative evaluation of β-amyloid. All the images and analysis presented in this paper are based on the ¹¹C-PIB imaging data acquired in an Alzheimer research study,¹⁶ however, the presented methodology is tracer-independent and can be used for other tracers, imaging protocols, and diagnostic objectives.

One of the challenges of FADS in dynamic emission tomography is the initial choice of the factors.¹⁷ This paper presents our modification of the factor analysis-based approach, CIFA, with the novel way of initializing the factors. We select the starting values of the factors using a representative clustering algorithm developed by us. It operates on the dynamic voxels similar to the independent-component analysis¹⁸ and principal component analysis¹⁹ techniques. Our approach establishes a number of small clusters in the space of dynamic image voxels, so that the elements of one cluster are as distinctive from the elements of other clusters as possible, thus corresponding to different dynamics. Clustering helps to estimate the number of factors for the FADS calculation, and average intensities for different clusters are used to initialize the factors.

The paper is organized as follows. First, we explain the theory of CIFA components and the complete algorithm and describe the dynamic ¹¹C-PIB PET experiments to which we apply our algorithm. Then, the algorithm is applied to 13 human subject data sets and our estimates of β-amyloid in the brain are compared with the results obtained from the Logan graphical analysis. The results are followed by the discussion of the new method and by conclusions of our work.

MATERIALS AND METHODS

CIFA: Clustering-Initiated Factor Analysis

Clustering-initiated factor analysis is an extension of the FADS algorithm obtained by adding a specially designed clustering method for initial processing of dynamic data.

In emission imaging, the unknown radiotracer concentration is a distribution V(r), where r is the three-dimensional position vector. In a discrete representation, V(r) is an array of N voxel values V_n, n = 1,2,... N. In dynamic imaging, these variables also depend on the time t, becoming V(r, t) and V_n(t), respectively. If the time is discretized into time frames t_k, k =1,2,...K, the dynamically changing voxels are denoted as V_n(t_k) or V_nk. In dynamic PET, the array V_nk is typically computed from projections by reconstructing a sequence of K arrays V_n before tissue classification, kinetic modeling, or other analysis. The principal challenges of the subsequent analysis are tightly related to the imperfections in the dynamic voxel reconstructions: high noise content, large voxel size causing partial volume effects, and insufficient temporal resolution.

Factor analysis of dynamic structures allows partial mitigation of the challenges in dynamic imaging. In FADS, the tracer concentration dynamics in the nth voxel is approximated by a linear expansion in terms of several time-dependent factors f_j with elements $f_{j k} \equiv f_{j} (t_{k})$ . The expansion coefficients C_nj are stationary:

V_{n} (t_{k}) \equiv V_{n k} = \sum_{j = 1}^{J} C_{n j} f_{j k} .

(1) The number of factors J is significantly smaller than both temporal and spatial dimensions of the system:

J < < K, J < < N .

(2) Both the N × J array C and the K × J array f are the unknowns and have to be determined either from the pre-reconstructed images V_nk or directly from the dynamically acquired projections. In this paper, we choose the former approach and determine C and f using multi-dimensional optimization by minimizing the following cost function:

X^{2} = \sum_{n = 1}^{N} \sum_{k = 1}^{K} {(\sum_{j = 1}^{J} C_{n j} f_{j k} - V_{n k})}^{2} .

(3) This optimization problem can be solved for some types of dynamic data provided proper sizes of the temporal bins. Large time bins can mask the differences in the underlying TACs to be reconstructed. Overly small time bins can result in high relative noise in the image data, thus reducing the stability of the solution. Most of the modern scanners including the one used by our group can acquire projections in list-mode format (as time-stamped individual counts), so the temporal bins can be modified to the desired optimum durations that satisfy the described limitations.

The key challenges in FADS include non-uniqueness and the underdetermined nature of the optimization problem, ambiguity in the choice of the optimal number of factors J, and high likelihood of converging to local minima that is aggravated by the numerical complexity of the solution related to the large problem size. Typically, the number of voxels N is > 10⁶. The total number of the unknowns is (N+K)xJ, so one way to reduce the size of the problem is by selecting a small J, yet this can reduce the robustness of the method by affecting the accuracy of the approximation (1). In light of the outlined challenges, the crucial step in the solution is selecting the optimal number of factors and their initial values, or factor initialization. We propose to initialize the factors using a dedicated clustering approach.

Representative clustering (R-clustering) is used to estimate the number and the initial value of the factors f by partial grouping of the dynamic voxel data. We would like to formulate the clustering as a general procedure applicable to a set {V_n} of N noisy time series $V_{n} \equiv {V_{n 1}, V_{n 2}, \dots V_{n k}}$ . For the sake of convenience, we presume N to be significantly larger than K, the length of the time series, and consider only non-negative V_nk. We make two principal assumptions about the input data. First, the time series in {V_n} can be relatively accurately described using the FADS expansion (1). Second, we assume that for each of the factors there exists a subset of {V_n}, for which, in the expansion (Equation 1), the contribution of this factor is prevalent or at least significant. Mathematically, this means that for a fixed factor index j, we can identify a subset of {V_n}, such that for the FADS expansion of each time series from that subset, the following inequality holds true:

C_{n j} \geq C_{n i} for all i such that i \neq j .

(4) Time-series that satisfy (4) are representative of the j^th factor. Let these time series form a set {R}_j of length P, so that (4) is true for each n ∈ {R}_j. When applied to dynamic images, the described assumption translates into the requirement that for each tissue there are at least P voxels, where that tissue is prevalent, so the mean intensity in these voxels is similar to the corresponding tissue TAC. Even though the number of time series that satisfy the condition (4) may be different for different j, we select a common minimum cluster size P for the whole system. Note that not all V_n are classified into clusters. For an assumed number of factors J, we select the cluster size P so that:

J \times P ≪ N .

(5) The goal of R-clustering is to identify J subsets of {V_n} that are most distinct from each other. To perform the clustering, we need to define the characteristic property of a cluster and to define a method of quantifying the distinction between time series. Each time series cluster {R}_j is characterized by r_j, the L₂-normalized average of its elements:

r_{j k} = \sum_{n} {(V_{n k} / \sqrt{\sum_{l = 1}^{K} V_{n l}^{2}})}^{2}, n \in {R}_{j}, k = 1, \dots, K .

(6) The L₂-normalization is needed to ‘equalize’ the input of each of the time series within the cluster and to prevent the average from being determined only by the elements with the largest magnitude. The validity of the FADS approximation implies that all r’ _js can be approximated by linear expansions of factors f and the validity of condition (4) means that the j^th factor is prevalent or at least essential in forming r_j.

The distinction between a time series V_n and a cluster {R}_j is to be measured using a distance operator D_jn:

D_{j n} = \sum_{k = 1}^{K} {(V_{n k} - z * r_{j k})}^{2}, where z * = arg {min}_{z} [D_{j n} (z)] .

(7) This operator focuses on the ‘functional’ difference between time series, whereas D between a time series that are scaled versions of each other is zero. Geometrically, D_jn denotes the square of the difference between a vector V_n and its projection onto direction r_j in K-dimensional Euclidean space.

R-clustering is a recursive process of selecting time series with maximal functional distance from previously computed clusters. In terms of D_jn, the j^th cluster is defined as a set of P dynamic series Vn with the largest total distance from the previously computed clusters. That is, for n ∈ {R}_j the sum $(\sum_{i < j} D_{i n})$ is larger than for any n ∉ {R}_j The first cluster in the recursive procedure is computed as the most remote and is defined by its distance from <, the weighted ensemble mean of the complete set {V_n}. When all the clusters have been established, the clustering is refined self-consistently by maximizing the total distances $(\sum_{i \neq j} D_{i n})$ .

The process of clustering is best explained using a visual analogy. Let each V_n be a point in K-dimensional space where the distance metric defined in (6) and (7) is used. The first cluster includes P points that are farthest away from the mean <V>. Each following cluster j includes P points that are farthest away from the mean of the (j –1) known averages. The new clusters are added using an iterative self-consistent approach: after each update, the cluster averages r_j are recomputed; the distances D_jn are recomputed; and the voxels are re-grouped anew.

The clustering is performed until either the pre-defined number of clusters has been reached or until adding a new cluster (J+1) produces an average r_J+1 that is not essentially different from the averages for the previous J clusters (can be described as a linear combination of r_j, j≾J). Although this is a subjective condition, it is simple to implement in practice, especially since the expected types of tissue TACs are known from biophysical considerations. There is no formal rule for selecting the cluster size P. The upper limit for P is set by (5) and the lower limit is set by the noisiness of the elements of {V_n}: P has to be sufficiently large so that the cluster averages r_j in (6) are not dominated by the noise. The simplest way to establish the cluster size is to select the smallest P such that by changing it to P+1, we do not alter the R-clustering-computed sequences r_j significantly.

After the self-consistent calculation for J clusters has been completed, the cluster averages r_j are used as the initial values of the factors f_j in the FADS calculation.

CIFA: Algorithm

The CIFA procedure is performed in six steps:

Start of the procedure.

Input data: Decay-corrected and normalized to the same time frame duration dynamic images, represented as N time series of length K with non-negative elements V_nk.

Step 1. Pre-filtering low-intensity voxels. Only voxels with significant tracer concentration are included in the main calculation. To do that, for each n, we compute a time series mean X_n and maximum Y_n:

X_{n} = < V_{n} > and Y_{n} = max (V_{n}) .

(8) Then, we select thresholds ε and ε, for both X and Y, and determine voxel indices n, for which

X_{n} < ε_{X} \times max (X) and Y_{n} < ε_{Y} \times max (Y) .

(9) These voxels are excluded from all the calculations in steps 2 through 5 to prevent factor contamination by noise. We denote the number of the remaining voxels (those for which (9) does not hold) by N'.

Step 2. R-clustering initialization. The procedure described in the previous section is performed for V_n. The parameters J and P for the number of clusters and the number of voxels per cluster, respectively, can be either pre-set or adjusted during the clustering process taking into account the size requirements in equations (2) and (5) . The computed sequences rj are used to initiate the zeroeth iteration factors

f_{j k} = r_{j k} for all j and k .

(10) The coefficient arrays are initialized uniformly or randomly taking into account two basic constraints on C_nj and one on f_jk so that the coefficients are unitless fractions of the tracer-uptake tissue types in each voxel, which should be non-negative and scaled to unity or less:

f_{j k} \geq 0.

(11)

C_{n j} \geq 0.

(12)

\sum_{j = 1}^{J} C_{n j} \leq 1.

(13) When all the variables are initialized, we compute the initial estimate of the cost function X² defined in (3).

Step 3. Refining C_nj for fixed f_jk. As the spatial and temporal sums in X² are independent, we can consider a separate cost function for each voxel n (same as in (3) but without the index n in the outer sum). This prevents application of any spatial regularization terms, but makes the calculation numerically simpler. In step 3, we loop through all voxels, performing N' separate K × J optimizations with each having J unknowns.

Step 4. Refining f_jk for fixed C_nj. This time, a separate cost function for each time frame k is considered. We perform minimization for each of K separate time frames, each of dimension N' × J, and with J unknowns each. At the end of this step, we recompute the cost function X². If the change in X² exceeds the convergence tolerance, we return to step 3, otherwise, proceed to step 5.

Step 5. Factor adjustment. The previous steps contain neither temporal nor spatial regularization for the sake of computational convenience. Although spatial regularization in the pre-filtered image space is non-trivial, it may be desirable to incorporate prior knowledge about the tissue TACs by enforcing temporal smoothness and, if needed, by setting minimum values for f_jk that converged to non-physically small values. Once these changes have been made, the algorithm returns to step 3 and the computation continues. If the ultimate value of the cost function is larger than that before making the changes, the variables should be reversed to their pre-step 5 values. In practice, we observe that making one or two such adjustments in the calculation may be advantageous by ‘guiding’ the optimization problem out of an a posteriori obvious local minimum.

Step 6. Inclusion of low-intensity voxels. For better visualization of the results, it is useful to compute the FADS expansion even for the low-intensity voxels pre-filtered in step 1. We initialize C_nj for all voxels that satisfy condition (Equation 9) with random numbers of the order of ε_X/J and perform step 3 optimization once for these voxels. This step does not change the factors. Knowledge of the continuous tissue distribution is needed primarily for better display. Also, this information may change the subsequent tissue analysis by adding quantitative information about low concentrations of specific binding tissue, which may be manifestations of new amyloid plaques.

Output data: Factor arrays f_j that approximate the TACs in the blood and in the tissues and factor images corresponding to coefficient arrays C that represent the spatial distribution maps of the corresponding tissue types.

End of the procedure.

Step 3 is the most resource-consuming part of the calculation, its computation time is directly proportional to the number of voxels N. To reduce the overall computation time and to ensure that only the brain voxels are used to determine the tissue factors, at first steps 2 through 5 of CIFA are applied to a smaller subset of the brain voxels. Once the factors are computed for the brain tissues, one to two repetitions of steps 2 through 5 for the larger set of voxels are sufficient to refine the blood factors and to compute the factor images. The initial subset is defined by selecting only brighter voxels (with higher threshold values) that are localized within the brain.

Tissue Analysis

Upon completion of the CIFA algorithm, its results need to be analyzed to produce a quantitative variable that is proportional to the amount of β-amyloid in the brain. First, the f'_js are studied visually to establish which factors and tissues correspond to each other. In the ¹¹C-PIB dynamic brain PET data, we can identify the following factors:

A blood curve (can be more than one).

A nonspecific binding tissue curve, characterized by both visible tracer uptake and tracer wash-out. Let the index of the factor corresponding to this curve be denoted by q. Then the coefficients for nonspecific binding distribution are C_nq.

A specific binding tissue curve, characterized by the slower uptake and much slower wash-out. The factor index for this curve is denoted by s. The coefficients for the specific binding distribution are C_ns. Then, two variables are computed from the coefficient array C, i.e., integrals of concentrations of nonspecific and of specific binding tissues F_Q and F_S, respectively:

F_{Q} = \sum_{n = 1}^{N} C_{n q} / \sum_{n = 1}^{N} \sum_{j = 1}^{J} C_{n j} .

(14)

F_{S} = \sum_{n = 1}^{N} C_{n s} / \sum_{n = 1}^{N} \sum_{j = 1}^{J} C_{n j} .

(15)

In the above expressions, numerators denote the integrals of the selected tissue type coefficients over the reconstruction volume. The denominators denote the integral of all tissue coefficients, they serve as normalization factors that allow interstudy comparison. F_S is expected to be proportional to the amount of the β-amyloid, as it is the only specific-uptake tissue. As the nonspecific binding dynamics is expected to be present relatively uniformly in the brain, the F_Q in (14) is linked to the total brain volume. Therefore, the following ratio is expected to be proportional to the amount of the specific-uptake tissue in the brain:

[Fraction of β - amyloid] ~ F_{S} / F_{Q} .

(16)

Experiment

Clustering-initiated factor analysis was tested on dynamic PET images for n =13 subjects, including four Alzheimer's disease patients and nine healthy controls. All research participants provided written informed consent. The data collection was performed under protocols approved by the Lawrence Berkeley National Laboratory Human and Animal Research Committee Institutional Review Board in accordance with the US federal research guidelines and regulations. The images were acquired on a Biograph 6 True-point PET/computed tomography camera by Siemens (Erlangen, Germany). Each subject was imaged using the same protocol. The subject was placed in the camera with the head fixed in the middle of the camera field of view. A bolus injection of approximately 15 mCi of ¹¹C-PIB was injected into an antecubital vein. The list-mode PET acquisition started immediately before the injection and continued for 90 minutes. The X-ray computed tomography scan was performed immediately before the emission acquisition.

After the scan, the list-mode data were binned into 35 time frames: 4 × 15 seconds, 8 × 30 seconds, 9 × 1 minutes, 2 × 3 minutes, 10 × 5 minutes, 2 × 10 minutes and reconstructed using an ordered subset expectation maximization algorithm with attenuation correction (four iterations, 16 subsets), Gaussian smoothing of 4 mm, zoom of 2, and scatter corrected using the Siemens scanner software. The reconstructed images were normalized to correct for the radioactive decay and for the different time frame durations. The voxel size of the reconstructed images was 1.0182 mm in-plane and 2.027 mm axially. In Figure 1, we show sagittal maximum-intensity projection images for four of the time frames for one of the subjects, a recent Alzheimer's disease diagnosis.

Figure 1.

Sagittal maximum-intensity projection (MIP) for four time frames, left-to-right: 15-second frame starting 30 seconds after injection, 30 seconds/3 minutes after injection, 1 minute/10 minutes after injection, and 5 minutes/30 minutes after injection.

The dynamic images were analyzed using the Logan graphical analysis, cerebellar gray matter was used as the reference tissue ROI, and the distribution volume ratio (DVR) was computed for each patient data set using the Logan graphical analysis^20,21 with frames corresponding to 35 to 90 minutes post injection. This value is expected to be proportional to the amount of β-amyloid in the brain.

RESULTS

Thirteen dynamic PET data sets were processed using the CIFA algorithm. The algorithm was implemented in Matlab 2012b using the linear least squares fitting function lsqlin to implement step 3 and the Optimization Toolbox function fmincon to implement step 4. Each input data set contained 35 reconstructed volumes of 336 × 336 × 109 voxels. First, the algorithm was applied to the ~ 120 × 120 × 30-voxel region that contained approximately 75% of the brain. On reconstructing the factors, an ~ 200 × 200 × 100 volume that contained the complete patient head was analyzed. (The exact region sizes varied on the basis of the size of the subject's head.) The following adjustable parameters were used:

Pre-filtering thresholds: ε_X = ε_Y =0.02.

Number of clusters/factors: J = 4.

Cluster size in R-clustering: P = 50.

Cost-function X² convergence tolerance: 1%.

Of these free parameters, the exact value of importance was only J. Varying P between 50 and 100 and the filtering thresholds between 0.02 and 0.05 or reducing X² threshold did not produce significant variations in the results. Each data set was analyzed for approximately 8 hours on a 2009 MacPro desktop computer. Over 80% of the processor time was taken up by the step 2 optimization, a computationally trivial process that can be easily parallelized.

For all the 13 patients, CIFA produced meaningful TAC-like factors and the coefficient distributions. Figure 2 shows R-clustering-generated curves and the same curves refined by FADS for the same subject as shown in Figure 1. The algorithm recovered two similar but separate curves for the tracer concentration in the blood with a time lag of 15 seconds or less. As the blood circulation time in the head is approximately 10 to 12 seconds, these peaks are likely to capture the tracer concentration ‘upstream’ and ‘downstream’ in the head vasculature, possibly showing the arterial and venous blood vessels. Indeed, the color-matched surface renderings of the corresponding coefficient arrays in Figure 3 show that the first (red) factor is concentrated in carotid arteries, while the second (blue) factor image appears to show a mass of smaller vessels in the brain, venous sinuses, and facial/nasal vasculature. The artery-vein separation is approximate and varies slightly from study to study, likely influenced by the patient's anatomy, physical state, and by the relative timing of the injection and PET data acquisition. In Figure 2A, we can observe that the R-clustering curves that have become the blood factors are non-zero after the first minute or two, which can be attributed to the contribution of the tissue caused by the partial volume effect. It should be noted, that as with other image-based IF detection techniques, CIFA recovers the total amount of radioactivity in the blood, not only the plasma tracer concentration.

Figure 2.

(A) R-clustering curves for one of the AD patients. (B) Factor curves computed by CIFA. Two tissue curves and two blood curves can be distinctively identified. The two blood curves partially correspond to the arterial and venous blood. AD, Alzheimer's disease; CIFA, clustering-initiated factor analysis; R-clustering, representative clustering.

Figure 3.

Surface rendering of blood coefficient distributions for three of the patients, two views each: front-right and front-back. Rendering thresholds are selected for optimal viewing at ~25% of the maximum. The outline of the head is obtained by thresholding the temporal mean of the dynamic PET images at 2% of the maximum. (A) The same subject as in Figures 1 and 2. The ‘upstream’ factor surface clearly corresponds to the internal carotid arteries. The ‘downstream’ factor surface is inside the brain, likely to be the interior cerebral veins and capillaries with increased concentration at the base of the sinuses. (B and C) Blood factors in two other patients. Carotid arteries, superior sagittal sinus, and transverse sinuses can be identified. PET, positron emission tomography.

The distributions of the specific and nonspecific binding tissues are shown in Figure 4 along with the DVR images obtained using standard graphical analysis method. In the cerebral cortex, both tissues are present, whereas in the cerebellum region the concentration of the nonspecific binding tissue is the highest and the specific binding tissue is expectedly absent there. The distribution of the specific binding tissue within the brain closely matches the DVR images.

Figure 4.

Tissue coefficient distributions for the same patient as in Figures 1, 2, and 3A. (A and C) Surface rendering and tomographic slices of the nonspecific binding tissue. (B and D) Rendering and slices for the specific binding tissues. The surface rendering is performed at 66% of the maximum values. (E) Tomographic slices of the DVR image obtained using the graphical analysis approach. The tick marks at the edges of the tomographic slices indicate relative locations of the slices. DVR, distribution volume ratio.

For each of the 13 data sets, we computed the integral tissue estimates F_Q and F_S defined in equations (14) and (15) and used these variables to estimate the fraction of β-amyloid in the brain defined in equation (16). This estimate was compared with the DVR value obtained from averaging multiple cortical ROIs computed using the reference tissue model. The comparison plot is shown in Figure 5A. The Pearson correlation coefficient for the two quantities is r=0.681.

Figure 5.

(A) Correlation of F_S/F_Q vs. the DVR computed using the Logan graphical method with the reference tissue region. (B) The same correlation with F_S/F_Q computed after applying a tissue mask that isolates only the brain. The data point marked by a red circle corresponds to the patient whose PET-image data were used in Figures 1, 2, 3A, and 4. CIFA, clustering-initiated factor analysis; DVR, distribution volume ratio; PET, positron emission tomography.

As both the FS/FQ ratios and the DVRs were expected to be proportional to the distribution of the same compound, the correlation r=0.68 was found to be insufficient. One source of the mismatch was the fact, that when computing F_Q and F_S, the whole head volume was used. Examining the factor images shows that for some patients, a considerable amount of the specific binding tissue was present outside of the brain. To correct for this effect, we used a spatial mask based on thresholding the factor image for the nonspecific binding tissue, which refocuses the analysis only to the brain region. F_Q and F_S were recomputed using only the voxels where the nonspecific binding coefficient values were more than 0.1. Then, the correlation increased to r =0.784. The adjusted correlation plot is shown in Figure 5B.

DISCUSSION

To the best of our knowledge, this paper presents one of the first systematic applications of FADS-based image analysis to ¹¹C-PIB β-amyloid dynamic brain PET. The new CIFA technique provides a tool for evaluating the amount of specific binding tissue in the brain with minimal participation of a human operator. Currently, the operator interference is needed only for (a) selecting the volume fraction for the initial calculation, (b) identifying blood and tissue factors from the segmentation results, and (c) semi-manual smoothing of the blood factors after the first pass of the algorithm implementation. All these tasks can be implemented in computer code, thus providing a completely automated tool for the analysis of specific binding tissue in dynamic PET using ¹¹C-PIB or other tracers.

The main novelty of CIFA is adding the R-clustering algorithm to FADS, which reduces the computation time and increases the overall effectiveness of the method by providing a fast and accurate initial guess of the tissue TACs. Performing FADS with randomly or uniformly initiated factors on several test data sets required 8 to 10 iterations and sometimes converged to unrealistic local minima, whereas two to four iterations were sufficient for basic convergence after R-clustering. This is significant since clustering requires approximately 10 to 20 times less computational time than a single FADS iteration. Also, FADS initialized with R-clustering was noticeably less likely to converge to local minima. Finally, R-clustering was instrumental in making the decision to use J= 4 factors, thus allowing us to recover partial separation of arteries and veins.

Neither the clustering nor the factor analysis in CIFA were subjected to any spatial or temporal regularization, only basic non-negativity (11 to 12) and scaling (13) constraints were used. This simplified and accelerated the computation, and provided a strong evidence of the validity of our approach. Nevertheless, we feel that CIFA can be further improved by incorporating regularization in the temporal domain. A smoothness constraint may be applied to the factors almost everywhere excluding the sharp blood peaks during the injection. Although in some reported applications of FADS,^6,7,17 smoothness, tissue overlap, or other spatial regularization constraints were included in the optimization, we think that in amyloid imaging such constraints may not be productive. Unlike in cardiac imaging, large volumes of blood are not present. Moreover, as most blood vessels are smaller than the voxel size, very few voxels contain only the blood factors. The nonspecific tissue may occupy large volumes, but is not completely homogeneous. Finally, the specific-uptake tissue is expected to be highly nonhomogeneous by nature. For example, the β-amyloid is a distribution of plaques with size ranging from a tenth of a micron to tens of microns.

Even though knowledge of the TACs and of the input function can be used in kinetic modeling, at this stage, our claims are limited to tissue classification and to quantitation of the specific binding tissue. We believe that the fraction of this tissue computed using the coefficient concentration arrays is potentially a better measure of the amount of tracer uptake in the tissue than currently used variables of DVR or standardized uptake volume ratio. Unfortunately, at this stage, we can only compare our estimate of the amount of β-amyloid to these variables. Therefore, we cannot estimate how much of the discrepancy between the two approaches can be attributed to the imperfections of the CIFA implementation and how much can be attributed to the weaknesses of the reference tissue model and implementation of the graphical analysis methods. The primary focus of our future work on CIFA is validation of its accuracy in evaluating the amount of irreversible uptake tissues by comparison with independent metrics. More work should be done for incorporating physiologically meaningful constraints in the optimization portion of CIFA. Upon initial validation, the use of a priori knowledge about the spatial properties of the tissue distributions during reconstruction will be expanded. It has been shown that spatial de-noising techniques can lead to significant improvements in kinetic analysis of ¹¹C-PIB dynamic PET images²² even without applying factor analysis. Ultimately, CIFA can be used for direct reconstruction of the kinetic parameters.²³

Amyloid imaging has been widely applied in the study of aging and dementia. Currently three PET radioligands are approved by the US Food and Drug Administration. Although the clinical application of these tracers is still somewhat controversial in the absence of an effective treatment, the tracers are also used widely in clinical trials for the development of new anti-β-amyloid drugs.²⁴ These studies require precise measurements including the detection of longitudinal change related to disease progression and treatment. Such measurement techniques are also crucial for the development of new diagnostic and research tracers in other brain PET areas, a popular research direction in modern diagnostic imaging.^25,26 Technical problems related to changes in the reference region can introduce noise and bias in the results of longitudinal studies. While the dynamic nature of the scan acquisition required for our approach certainly introduces complexity into PET studies, this may be offset by greater precision that could benefit disease detection and drug development.

The main potential advantages of CIFA in comparison with the graphical analysis methods are related to its a priori independence of the kinetic model and a lesser dependence on the manual input. For example, factor analysis relies only on the fact that the two tissue TACs are different, whereas in the Logan plot analysis, both curves are assumed to have achieved and are sustained at a steady state. Also, our method is not sensitive to the precision and reproducibility of the choice of the reference region or other results of anatomic classification. At the same time, CIFA has its own potential weaknesses. Its accuracy is tightly related to the accuracy of the solution to the FADS optimization problem. Algorithm convergence to a local minimum or a bad choice of the number of factors can change its estimate of the specific binding tissue. Finally, CIFA has not been studied as extensively as the Logan analysis approach. Although it is likely to be a good alternative method for tissue classification from dynamic PET images, extensive further research and validation has to be performed before it can be adopted for regular application in the field.

CONCLUSIONS

We have developed the CIFA procedure, which combines basic factor analysis with a factor initialization technique based on a novel R-clustering method. Clustering-initiated factor analysis was successfully implemented and applied to a series of dynamic ¹¹C-PIB PET data sets. The ability of the method to extract the input function and to detect some aspects of the time lag between the arterial and the venous blood was demonstrated. The time–activity curves for the tissues with specific binding and nonspecific binding uptake were obtained. The factor images (coefficient distributions) for two blood factors exhibited strong visual correspondence to the anatomic features of the brain vasculature, providing indirect validation of the method. The factor images for the tissue factors were used to estimate the fractions of specific and nonspecific tissue types in the brain. The ratio of these two fractions was expected to be proportional to the amount of the specific binding tissue. On the basis of 13 patient scans, the CIFA-derived estimate was compared with the DVR evaluated using the Logan graphical analysis and the reference tissue model. The two variables exhibited significant proportionality. The proposed technique may prove to be a better way of quantitative estimation of the amount of specific binding tissue in dynamic PET than currently used techniques; however, it needs to be further improved and validated by direct comparisons with PET-independent methods.

AUTHOR CONTRIBUTIONS

The original development of CIFA, its software implementation, and its application to patient data processing, as well as the principal manuscript preparation were performed by RB. DM and GTG contributed to the algorithm development (R-clustering and FADS parts, respectively) and to manuscript preparation and revision. SLB and WJJ performed all the activities related to acquiring, reconstructing, and conventional processing of the patient dynamic PET data and prepared the portions of the manuscript related to the imaging experiments and to neuroscientific motivation of our work.

DISCLOSURE/CONFLICT OF INTEREST

WJJ has served as a consultant for Genentech and Synarch. The remaining authors declare no conflict of interest.

References

Cousins

O'Donnell

. Arterial cannulation: a critical review. AANA J 2004; 72: 267–271.

Zanotti-Fregonara

Chen

Liow

J-S

Fujita

Innis

. Image-derived input function for brain PET studies: many challenges and few opportunities. J Cereb Blood Flow Metab 2011; 31: 1986–1998.

Arbelaez

Benzinger

TLS

Snyder

Vlassenko

Mintun

. Noninvasive estimation of the arterial input function in positron emission tomography imaging of cerebral blood flow. J Cereb Blood Flow Metab 2013; 33: 115–121.

Schain

Benjaminsson

Varnäs

Forsberg

Halldin

Lansner

. Arterial input function derived from pairwise correlations between PET-image voxels. J Cereb Blood Flow Metab 2013; 33: 1058–1065.

Kao

Chen

Wernick

. Statistical analysis of dynamic sequences for functional imaging. SPIE 2000; 3978: 347–353.

Sitek

Gullberg

Huesman

. Correction for ambiguous solutions in factor analysis using a penalized least squares objective. IEEE Trans Med Imag 2002; 21: 216–225.

El Fakhri

Sitek

Guérin

Kijewski

Di Carli

Moore

. Quantitative dynamic cardiac 82Rb PET using generalized factor and compartment analyses. J Nucl Med 2005; 46: 1264–1271.

El Fakhri

Trott

Sitek

Bonab

Alpert

. Dual-tracer PET using generalized factor analysis of dynamic sequences. Mol Imag Biol 2013; 15: 666–674.

Hoh

Choi

Schelbert

Hawkins

Phelps

. Factor analysis for extraction of blood time–activity curves in dynamic FDG-PET studies. J Nucl Med 1995; 36: 1714–1722.

10.

Logan

. A review of graphical methods for tracer studies and strategies to reduce bias. Nucl Med Biol 2003; 30: 833–844.

11.

Ichise

Toyama

Innis

Carson

. Strategies to improve neuroreceptor parameter estimation by linear regression analysis. J Cereb Blood Flow Metab 2002; 22: 1271–1281.

12.

Zhou

Brasic

Wong

. Multigraphical analysis of dynamic PET. Neuroimage 2010; 49: 2947–2957.

13.

Salinas

Searle

Gunn

. The simplified reference tissue model: model assumption violations and their impact on binding potential. J Cereb Blood Flow Metab 2015; 35: 304–311.

14.

Mucke

. Alzheimer's disease. Nature 2009; 461: 895–897.

15.

Jack

Knopman

Jagust

Petersen

Weiner

Aisen

. Tracking pathophysiological processes in Alzheimer's disease: an updated hypothetical model of dynamic biomarkers. Lancet Neurol 2013; 12: 207–216.

16.

Wirth

Villeneuve

La Joie

Marks

Jagust

. Gene-Environment interactions: lifetime cognitive activity, ApoE genotype and beta-amyloid burden. J Neurosci 2014; 34: 8612–8617.

17.

Abdalah

Boutchko

Mitra

Gullberg

. Reconstruction of 4D dynamic SPECT images from inconsistent projections using a spline initialized FADS algorithm (SIFADS). IEEE Trans Med Imag 2015; 34: 1–13.

18.

Naganawa

Kimura

Nariai

Ishii

Oda

Manabe

. Omission of serial arterial blood sampling in neuroreceptor imaging with independent component analysis. Neuroimage 2005; 26: 885–890.

19.

Kimura

Senda

Alpert

. Fast formation of statistically reliable FDG parametric images based on clustering and principal components. Phys Med Biol 2002; 47: 455–468.

20.

Logan

Fowler

Volkow

Wang

Ding

Alexoff

. Distribution volume ratios without blood sampling from graphical analysis of PET data. J Cereb Blood Flow Metab 1996; 16: 834–840.

21.

Price

Klunk

Lopresti

Joge

Ziolko

. Kinetic modeling of amyloid binding in humans using PET imaging and Pittsburgh Compound-B. J Cereb Blood Flow Metab 2005; 25: 1528–1547.

22.

Floberg

Mistretta

Weichert

Hall

Holden

Christian

. Improved kinetic analysis of dynamic PET data with optimized HYPR-LR. Med Phys 2012; 39: 3319–3331.

23.

Wang

. An Optimization transfer algorithm for nonlinear parametric image reconstruction from dynamic PET data. IEEE Trans Med Imag 2012; 31: 1977–1988.

24.

Salloway

Sperling

Fox

Blennow

Klunk

Raskind

. Two phase 3 trials of bapineuzumab in mild-to-moderate Alzheimer's disease. New Engl J Med 2014; 370: 322–333.

25.

Naganawa

Nabulsi

Planeta

Gallezot

Lin

Najafzadeh

. Tracer kinetic modeling of [(11) C] AFM, a new PET imaging agent for the serotonin transporter. J Cereb Blood Flow Metab 2013; 33: 1886–1896.

26.

Mach

. New targets for the development of PET tracers for imaging neuro-degeneration in Alzheimer disease. J Nucl Med 2014; 55: 1221–1224.

Clustering-Initiated Factor Analysis Application for Tissue Classification in Dynamic Brain Positron Emission Tomography

Abstract

Keywords

INTRODUCTION

MATERIALS AND METHODS

CIFA: Clustering-Initiated Factor Analysis

CIFA: Algorithm

Tissue Analysis

Experiment

RESULTS

DISCUSSION

CONCLUSIONS

AUTHOR CONTRIBUTIONS

DISCLOSURE/CONFLICT OF INTEREST

References