Analysis of Dynamic Radioligand Displacement or “Activation” Studies

The basic idea

This section provides an intuitive idea of how the analysis proceeds. The general strategy is described, followed first by a more formal description of a “minimal” example. This example is then used in the analysis of empirical data.

In the next section, we consider a simple example of this general approach using a single-compartment model with a single form of time-dependent change in one kinetic parameter (k₂). This minimal model is the model used in the analysis of the empirical data.

The statistical model: An example

In this section, we present the statistical model used to assess the degree and significance of dynamic changes in kinetics during a ¹¹C-Flumazenil study. Because a single compartment (monoexponential) model is sufficient to describe the tissue dynamics of ¹¹C-Flumazenil (Koeppe et al., 1991), the corresponding statistical model was based on a single compartment. The observed variable F_ij corresponds to the tissue activity (per frame) measured at voxel i in frame j:

F i j = ∫ C i ( t ) + V b C b ( t ) d t + ( W j ) 1 / 2 ( r i j + o i )

(1)

where scan j starts at time t_j and finishes at t_j+1 and the integral is over this period. C_b(t) is the whole blood activity at time t, and V_bC_b(t) represents the small contribution to the measured activity from the blood compartment with volume V_b (fixed at 0.1 for all voxels). (W_j)^1/2 (r_ij + o _i ) is a residual term appropriate for Poisson statistics, where W_j is the whole brain activity for scan j, and r_ij is assumed to be identically and independently distributed as N(0, σ²); o _i is a small, constant background term. The time-dependent tissue activity C_j(t) is given by the differential equation:

d C i ( t ) / d t = κ 1 i C p ( t - δ i ) - κ 2 ( t ) i + λ ) . C i ( t )

(2)

κ_1i is the plasma to tissue transfer rate for voxel i. C_p(t – δ _i ) is the metabolite corrected plasma curve, or delayed input function, and δ _i is the delay incurred by its estimation in the periphery relative to brain region i. An example of C_p(t) is shown in Fig. 1 (upper panel). These data were obtained from the validation study described in a subsequent section. λ is the radiodecay parameter for the tracer (in this case, 0.0005663 s^–1). κ₂(t) _i corresponds to the rate of clearance from the tissue compartment for voxel i at time t. The clearance κ₂(t) _i is a time-dependent variable that may exhibit changes due to the activation [note that time-dependent changes in κ_1i would have little impact later in the study because C_p(t – δ _i ) will be small]. It is these activation-induced changes that we are interested in: Let κ₂(t) _i be expanded in terms of some regionally invariant “normal” value (κ₂*), a region or voxel-specific effect (Δκ_2i), and a voxel-specific, time-dependent component due to the activation Γ _i (t). The dynamic component due to changes in displacement Γ _i (t) can be expressed in terms of some appropriate function of time, say, h(t) so that Γ _i (t) = γ_i · h(t) (see Discussion for a more general expansion). The form of h(t) defines the time course of displacement dynamics about which one wants to make inferences. In this paper, we assume that the effect of the activation is to induce a fairly protracted increase in displacement κ₂ or equivalently a decrease in the volume of distribution (see below). The coefficient γ_i is the unknown that represents this effect of interest and under the null hypothesis of no time-dependent displacement γ_i = 0:

κ 2 ( t ) i = κ 2 * + Δ κ 2 i + γ i . h ( t )

and similarly

κ 1 i = κ 1 * + Δ κ 1 i δ i = δ * + Δ δ i

(3)

OPEN IN VIEWER

FIG. 1.

Upper panel: Whole blood C_b(t) and input function or metabolite corrected plasma activity C_p(t) as a function of time. These data were based on whole-blood activity monitored every second. The whole-blood data were transformed and metabolite corrected on the basis of discrete arterial samples. Lower panel: A typical curve for the canonical form of measured counts per frame. This curve was obtained by using C_b(t) and C_p(t) shown in the upper panel and assumes 64-s time frames throughout the study. Values for parameters were as follows; κ₁* = 0.0077/s, κ₂* = 0.0045/s, δ* = 16 s, V_b = 0.1 and λ = 0.0056/s.

For ¹¹C-Flumazenil, we could assume a value of κ₂* that is typical of gray matter (e.g., 0.0045) or estimate κ₂* by using the whole brain curve for the study in question and similarly for κ₁* and δ*. Combining Eqs. 1, 2, and 3, we can express expected activity at voxel i in frame j, (E _ij ) as a function of four unknown parameters:

E i j ( Δ κ 1 i , Δ δ i , Δ κ 2 i , γ i ) = ∫ C i ( t ) + V b C b ( t ) d t

(4)

Integrating over the time of the frame. If we make the simplifying assumption that the deviations from the canonical values (Δκ_1i, Δδ _i , Δκ_2i, γ_i) at each voxel are small, we can take a first-order approximation of a Taylor expansion of E_ij around E _j * = E_ij(0, 0, 0, 0) where

E i j ≈ E j * + Δ κ 1 i ∂ E j * / ∂ Δ κ 1 + Δ κ 2 i ∂ E j * / ∂ Δ κ 2 + Δ δ i ∂ E j * / ∂ Δ δ + γ i . ∂ E j * / ∂ γ

(5)

A typical curve for the canonical form is shown in Fig. 1 (lower panel). This curve was obtained by using C_b(t) and C_p(t), shown in the upper panel of Fig. 1, and assuming 64 second time frames throughout the study. Values for each of the parameters were as follows; κ₁* = 0.0077 per second, κ₂* = 0.0045 per second, δ* = 16 second, V_b = 0.1 and λ = 0.0056/s.

In practice, E _j * and the partial derivatives are most simply obtained by numerical integration of Eq. 3: For example, the partial derivative ∂E_j*/∂Δκ₁ is estimated with [E_ij(p, 0, 0, 0) – E_ij(0, 0, 0, 0)]/p, where p is some suitably small perturbation. Examples of these partial derivatives are presented in Fig. 2. The upper panel shows the partial derivatives with respect to Δκ_1., Δδ, and Δκ₂. These profiles can be thought of as the deviation from the canonical form in Fig. 1 that would ensue if Δκ_1i, Δδ _i , or Δκ_2i changed slightly. The fourth partial derivative ∂E_j*/∂γ modelling a time-dependent increase in κ₂ is seen in the lower panel of Fig. 2 (solid line). The assumed form for changes in κ₂ is superimposed (broken line) and corresponds to h(t). A transient increase in κ₂ produces an accelerated decay in observed counts; it is this perturbation of the dynamics that the model is designed to assess. The curves have been scaled to allow comparisons among them.

OPEN IN VIEWER

FIG. 2.

Upper panel: Examples of the partial derivatives with respect to Δκ, Δδ, and Δκ₂. These profiles can be thought of as the deviation from the canonical form in Fig. 1 (lower panel) that would ensure if Δκ_1i, Δδ _i , or Δκ_2i. changed slightly. Lower panel: The fourth partial derivative ∂E _j */∂γ modelling a time-dependent increase in κ₂ (solid line). The assumed form for changes in κ₂ is superimposed (broken line) and corresponds to h(t). A transient increase in κ₂ produces an accelerated decay in observed counts; h(t) is a Gaussian curve (of parameter 64 s) centred at 2,100 s and convolved with an exponential function with a half-life of 320 s (5 min).

Combining Eqs. 1 and 5, we can now define the response variable X_ij as:

X i j = ( W j ) - 1 / 2 ( F i j - E j * ) ≈ ( W j ) - 1 / 2 ( Δ κ 1 i ∂ E j * / ∂ Δ κ 1 + Δ κ 2 i ∂ E j * / ∂ Δ κ 2 + Δ δ i ∂ E j * / ∂ Δ δ + γ i ⋅ ∂ E j * / ∂ γ ) + r i j + o i

In matrix notation:

X i ≈ G ⋅ β i + r i

where

X i = diag ( W ) - 1 / 2 ⋅ ( F i - E * )

and

G = { 1 diag ( W ) - 1 / 2 ⋅ [ ( ∂ E * / ∂ Δ κ 1 ) × ( ∂ E * / ∂ Δ κ 2 ) ( ∂ E * / ∂ Δ δ ( ∂ E * / ∂ γ ) ] }

and

β i = [ o i , Δ κ 1 i , Δ δ i , Δ κ 2 i , γ i ] T

(6)

Eq. 6 is the general linear model with design matrix G (see Chatfield and Collins, 1980). X_i is a column vector of the response variables (counts per frame corrected for the canonical form and weighted) for each of the j scans. F _i is a column vector of counts per frame for voxel i. E* is a column vector of E _j *, similarly for W, ∂E*/∂Δκ_i and the remaining partial derivatives. The parameters to be estimated are the elements of the column vector β _i. Because of the highly constrained nature of the statistical model, G is of full rank and G^TG is invertible. It follows that the least-squares estimates of the parameters β _i (say, b _i ) are given uniquely by

b i = ( G T G ) - 1 G T X i

where

E { b i } = β i and Var { b i } = σ i 2 ( G T G ) - 1

(7)

b _i = [o _i , Δk_1i, Δd _i , Δk_2i, g _i ]^T is an estimate of [o _i , Δκ_1i, Δδ _i , Δκ_2i, γ _i ]^T. Var{b _i } is the variance-covariance matrix for the parameter estimates corresponding the ith voxel.

The parameter estimates

Before proceeding to statistical inferences, this section addresses the estimates of the kinetic parameters and time-dependent changes attributable to the “activation.” The form of Eq. 6 means that the solution for the parameter estimates can be effected simultaneously for all voxels. For any voxel i, the parameter estimates are b_i = [o _i , Δk_1i, Δd _i , Δk_2i, g_i]^T where, by Eq. 3 the rate constant κ₂ is given by:

k 2 ( t ) i = k 2 * + Δ k 2 i + g i . h ( t )

and similarly

k 1 i = k 1 * + Δ k 1 i

(8)

and the corresponding volume of distribution is given by:

V d ( t ) i = k 1 i / ( k 2 ( t ) i + λ )

(9)

The average relative change in k₂ is simply:

( g i / t J ) ∫ h ( t ) d t / k 2 ( 0 ) i

(10)

where t_J is the time over which the integral is taken. These estimates can be displayed as parametric images and provide a rough characterisation of the binding kinetics on a voxel by voxel basis. In this paper we use these results to validate the statistical model adopted (see below). Even if large effects are found, their significance is not assured. In the final theory sections, we address how to assess the significance of regionally specific activation-dependent changes in displacement.

Statistical inference

In this section, we address statistical inferences about the effects of interest (activation-induced changes in displacement or V_d) at a single voxel. Significance is assessed with the reduction in error variance that is seen when including the time-dependent effects or equivalently by testing the null hypothesis that g_i = 0. One could use the F ratio; but, for reasons developed in the next section, we suggest the t statistic obtained by using a vector that is 1 for the effect of interest (g_i) and zero elsewhere: c = [0, 0, 0, 0, 1]. It should be noted that, when testing one degree of freedom, the F statistic and the square of the t statistic are the same:

t i = c . b i / ɛ i

where the standard error of this compound, for voxel i, is estimated by ε _i ²:

ɛ i 2 = ( R i ( Ω ) / r ) c . ( G T G ) - 1 ⋅ c T

(11)

t_i has the Student's t distribution with degrees of freedom r = J – rank(G), where J is the total number of scans. The sum of squares due to error R_i(Ω) is obtained from the difference between the actual and estimated values of X_i:

R i ( Ω ) = ( X i - G ⋅ b i ) T ( X i - G ⋅ b i )

(12)

When displayed as an image, the t_i constitute a statistical parametric map, or SPM{t}. This SPM represents a spatially extended statistical process that reflects the significance of dynamic changes in displacement. We are now in a position to apply standard results from statistical parametric mapping (Friston et al., 1991; Worsley et al., 1992; Friston et al., 1994). The remaining analytic steps have been described in detail elsewhere but are summarised here for completeness.

The p value for a regional effect

In this section we address the problem of how to interpret the SPM{t} in terms of probability levels or p values. The problem here is that an extremely large number of nonindependent univariate comparisons have been performed, and the probability that any region of the SPM will exceed an uncorrected threshold by chance is high. In recent years, advances have been made that have solved this problem of multiple nonindependent comparisons. A local excursion of the SPM (a connected subset of voxels above some threshold) can be characterised by its maximal height (Z, the peak value) or by its size (n, the number of voxels that constitute the region). Both these simple characterisations have an associated probability of occurring by chance over the volume of the SPM. These two probabilities form the basis for making a statistical inference about any observed regional effect: The distributional approximations for Z and n derive from the theory of Gaussian random fields. To simplify the analysis the SPM{t} is transformed to a SPM{Z} using a probability integral transform or other standard device (although exact results are available for SPM{t}, using the maximal height inference, they have not yet been developed for the spatial extent inference).

The p value based on Z, the peak height of the region

The probability of getting at least one voxel with a Z value of, say, height u or more, in a given SPM{Z} of volume S is the same as the probability that the largest Z value in the entire volume is greater than u. It can be shown (see Adler, 1981):

P ( Z max > u ) ⩽ E { m } ≈ S ( 2 π ) - ( D + 1 ) / 2 w - D u D - 1 e - u 2 / 2

(13)

where E{m} is the expected number of maxima; wis a measure of smoothness and is related to the full width at half maximum (FWHM) of the SPM's resolution. In practice, w can be determined directly from the effective FWHM if it is known when w = FWHM/√(4log_e2) or estimated post hoc using the measured variance of the SPM's first partial derivatives. See Friston et al. (1991) and Worsley et al. (1992) for more details.

The p value based on n, the size of the region

The probability of getting one or more regions of, say, size k or more in a given SPM{Z} of volume S, thresholded at u, is the same as the probability that the largest region consists of k or more voxels P(n_max > k) where:

P ( n max ⩾ k ) = 1 - exp ( - E { m } ⋅ e - β k 2 / D )

and

β = [ Γ ( D / 2 + 1 ) ⋅ E { m } / S ⋅ Φ ( - u ) ] 2 / D

(14)

Φ(–u) is the error function (integral of the unit Gaussian distribution) evaluated at the threshold chosen (u). Eq. 14 gives an estimate of the probability of finding at least one region with k or more voxels in a SPM{Z}. In this paper, we present p values that are based on both the spatial extent and on peak height. The uncorrected p values are simply Φ(–u).

Data acquisition and experimental design

The data were obtained from a normal male subject in a manner approved by the local ethics committee and the Administration of Radioactive Substances Advisory Committee (UK). A radial arterial catheter was inserted with the subject under local anesthetic (bupivicaine 5%) into the nondominant arm and an i.v. cannula into the contralateral arm. Scans were obtained with a Siemens-CTI PET camera (model 953B CTI Knoxville, TN, U.S.A.) with performance characteristics as described by Spinks et al. (1992) and measured attenuation correction. ¹¹C-Flumazenil was prepared by a modification of the method of Maziere et al. (1984) and administered intravenously. About 6.7 mCi of activity was injected (248 MBq). Twenty-five emission scans were then collected. The duration of these frames ranged from 60 s in the initial scans to 900 s (3 × 1 min, 18 × 3 min and 3 × 15 min; in this particular study, the last frame was lost). Whole-blood radioactivity concentrations were monitored continuously with a bismuth germanate system connected to the arterial line. Discrete arterial samples were used to transform the whole-blood curve into a metabolite corrected well count calibrated plasma curve as described previously (Lammertsma et al., 1991).

Data preprocessing

The data were realigned and spatially normalised (Friston et al., 1995) (voxel size = 2 × 2 × 4 mm). The normalisation used a 12-parameter affine transformation that best matched the sum of all 25 frames to an appropriate template or model ¹¹C-Flumazenil image. This ¹¹C-Flumazenil template was taken from the SPM library of normalised images and was created by applying a nonlinear spatial normalisation (Friston et al., 1995) to an arbitrary normal ¹¹C-Flumazenil image (integrated over all frames) such that it conformed to the existing rCBF template. This normalising transformation solves simultaneously for a spatial and intensity transformation and accommodates differences in the regional distribution of tracers employed. Following normalisation the image was flipped and added to the un-flipped image and smoothed with a (8 mm FWHM) Gaussian kernel to provide a symmetrical ¹¹C-Flumazenil reference or template that conformed to the standard anatomical space (Talairach and Tournoux, 1988). Figure 3 shows the correspondence between the normalised integrated data (upper images) and the template (lower images). The data were then smoothed with an isotropic Gaussian kernel (FWHM of 8 mm). Only voxels exceeding 0.8 of the global activity for each frame separately were subject to further analysis. This is a standard threshold in SPM and selects voxels that are, predominantly, in gray matter.

OPEN IN VIEWER

FIG. 3.

Orthogonal sections through the ¹¹C-Flumazenil data (integrated over all 25 frames), upper set, and corresponding sections through a Flumazenil template image; lower set, used in spatial normalisation. These sections are through the anterior commissure [0,0,0] mm following spatial normalisation. These data demonstrate the efficacy and quality of the spatial normalisation procedure. The normalisation used a 12-parameter affine transformation as described in Friston et al. (1995).

Statistical analysis

The whole blood and metabolite corrected plasma activity curves C_b and C_p are depicted in Fig. 1 (upper panel). The response variables and design matrices were computed as described in the theory section using the following values: κ₁* = 0.0077/s, κ₂* = 0.0045/s, δ* = 16 s, V_b = 0.1, and λ = 0.0056/s, the same values used in Figs. 1 and 2. The form of h(t) used in the current analysis is seen in Fig. 2 (broken line). The analysis provided estimates of relevant parameters and t statistics. The SPM{t} was transformed to a SPM{Z} and the distribution of Z values evaluated. The first frame was added to the second to reduce sensitivity to delay and dispersion of the estimated input function.

Parameter estimates

Values for k₂(0) _i and V_d(0) _i were estimated using Eqs. 8 and 9 for all voxels. The k₂ estimates were converted to the corresponding half-life ln2/k₂. Maximum intensity projections were then constructed for ln2/k₂(0) _i and V_d(0) _i . Figure 4 (upper panel) shows the spatial distribution of ln2/k₂(0) _i. It can be seen that the half-life of ¹¹C-Flumazenil is fairly uniform within cerebral gray matter, but there are some interesting regional variations: the ventral aspects of the temporal cortex, from extrastriate to inferotemporal and temporopolar regions have relatively high values of ln2/k₂(0) _i. The highest values are seen in the left medial temporal lobes. The distribution of this parameter (Fig. 4, lower panel) shows the estimates to centre on 700 s, in good agreement with previous estimates of the half-life for Flumazenil of about 12 min (720 s) in gray matter. The equivalent results for V_d(0) _i are shown in Fig. 5. In this instance, the greatest values are seen in the occipital and posterior cingulate cortices in accord with the distribution of receptor density and binding potential. The highest values are in the striate cortex and are about 5. This value is slightly lower than would be expected using more conventional analyses (a value of about 7 is typical). This disparity is probably due to the fact that the data have been smoothed and the regional values contain both gray and white matter contributions.

OPEN IN VIEWER

FIG. 4.

Upper panel: Maximum intensity projection of the initial half-life associated with k₂, estimated as described in the text for the validation study. The display format is standard and provides three views of the brain from the front, back, and right-hand side. The gray scale is arbitrary, and the space conforms to that described in the atlas of Talairach and Tournoux (1988). Lower panel: Histogram of the relative frequency distribution of voxel values displayed in the upper panel.

OPEN IN VIEWER

FIG. 5.

Upper panel: Maximum intensity projection of the initial volume of distribution V_d(0) _i estimated as described in the text for the validation study. The format is the same as in the previous figure. Lower panel: Histogram of the relative frequency distribution of voxel values displayed in the upper panel.

In summary, the estimates appear to have some validity in terms of both the quantitative and anatomical distribution. This validity, however, does not imply that the model is statistically valid. To be confident that the model is appropriate for making statistical inferences, we must demonstrate that the statistics obtained have the expected distribution under the null hypothesis.

Identically and independently distributed error terms

Identically distributed error terms implies homoscedasticity. Homoscedasticity refers to the homogeneity of error variances over observations (i.e., frames) and is an assumption that should be verified given the Poisson nature of the underlying data. Clearly, we cannot estimate the error variance for each voxel in every scan (there is only one observation per voxel per scan); however, we can partially validate the assumption of homoscedasticity by showing that the variance in the error term over voxels is roughly constant over frames. The sums of squares resulting from error were calculated according to Eq. 12 and plotted as standard deviations for each of the frames analysed. The results (Fig. 6) demonstrate that, with the exception of the first frame, the error terms are reasonably well behaved (homoscedastic). The independence of the error terms in Eq. 1 is guaranteed because they come from different time frames.

OPEN IN VIEWER

FIG. 6.

Standard deviation of the error terms estimated over all voxels analysed in the validation study plotted as a function of frame number. Ideally, they should all be the same.

Distribution of the SPM{Z} values under the null hypothesis

The distribution of the Z values for all the voxels analysed is shown in Fig. 7. The histogram details the empirically derived distribution, and the solid line is that expected under the null hypothesis; there is a reasonable agreement. It is also evident that there is a slight but definite bias toward negative values. The reasons for this may include (1) forcing multicompartment kinetics into a single compartment model, (2) real decreases in κ₂ during the latter half of the study, or (3) the effect of smoothing. This last possibility is suggested by the observation that the identical analysis of “unsmoothed” data removed this bias (results not shown). It is possible that smoothing confounds tissue compartments with different values for κ₂ and other parameters. This would “smear” the true spectrum of parameters and render a single compartment (mono-exponential) model more approximate. If this is the case, it suggests some interesting extensions of the approach considered in this paper; however, it should be noted that although the distribution in Fig. 7 is less than perfect, the bias is in the conservative direction for detecting increases in κ₂.

OPEN IN VIEWER

FIG. 7.

Expected and observed distribution of Z values from the validation study SPM{Z}.

Study design and data acquisition

The data were obtained from a normal male volunteer and analysed in exactly the same way as in the previous sections. In this study, however, 50 μg/kg of midazolam was administered intravenously over 5 min between 30 and 35 min after beginning the study. The form chosen for the time-dependent change in κ₂[h(t)] was the same as that used in previous sections (the broken line in Fig. 2).

Parameter estimates

The parameter of interest in this study is g_i, the coefficient scaling the time-dependent displacement described by h(t). The maximum intensity projection of relative changes in k₂, calculated according to Eq. 10, is presented in Fig. 8 (upper panel). The corresponding distribution of relative changes (expressed as % change) shows these values to be largely positive with a maximum change of about 16%. These results pertain to the size of the effect but do not show whether these regional effects are significant or not. To make statistical inferences, one uses the corresponding SPM{Z}.

OPEN IN VIEWER

FIG. 8.

Upper panel: Maximum intensity projection of the average change in κ₂ estimated for the activation study. The display format is the same as in previous figures. Lower panel: The relative distribution of voxel values displayed in the upper panel.

Statistical inference

The SPM{Z} is depicted in Fig. 9. There was a remarkably significant effect (p < 0.01 uncorrected) in nearly all cerebral cortical gray matter regions. The highest Z score reached 5.25 in the left inferior occipitotemporal region. The next two most significant displacement effects were seen bilaterally in the left (Z = 5.06 p = 0.007 corrected using Eq. 13) and right (Z = 5.05 p = 0.008 corrected using Eq. 13) dorsolateral prefrontal cortices (DLPFC). The FWHM of the SPM{Z} were estimated to be 8.7, 8.9, and 9.9 mm in x, y, and z respectively. A total of 41,999 voxels were analysed.

OPEN IN VIEWER

FIG. 9.

Upper panel: Statistical parametric map SPM{Z}. This is a maximum intensity projection of the SPM{f} following transformation to the Z score. This SPM assesses the significance of the effects seen in the previous figure (Fig. 8). Lower panel: Tabular data are presented on regional effects in the upper panel. The locations of maxima within the largest region (10,918 voxels) are shown. For each maxima (Z), the significance is assessed in terms of P(Z_max > u) using Eq. 13 and its coordinates are provided.

The tissue curves for a voxel in the left dorsolateral prefrontal cortex or DLPFC (–36, 34, 24 mm) are shown in Fig. 10 superimposed on the curves predicted by the parameter estimates (solid line). The dotted line indicates what would have happened in the absence of any displacement. For comparison, the data for the same voxel from the validation study are shown in Fig. 11.

OPEN IN VIEWER

FIG. 10.

The observed (circles) and predicted activity (solid line) for a voxel in the left prefrontal cortex [–36, 34, 24 mm]. By omitting the component due to time-dependent changes, we can estimate the predicted activity in the absence of any dynamic changes in displacement (dotted line). The vertical broken line corresponds to the end of the midazolam injection.

OPEN IN VIEWER

FIG. 11.

As for the previous figure (Fig. 10), but using data from the validation study.

Figure 12 shows the SPM[Z} rendered onto the Flumazenil template in three orthogonal sections through this DLPFC voxel. The bilateral involvement of the DLPFC, temporopolar regions, and extrastriate regions is evident in these sections and highlights some of the regional specificity and symmetry of the more significant effects.

OPEN IN VIEWER

FIG. 12.

Orthogonal sections through the Flumazenil template (gray) with the statistical parametric map SPM{Z} (colour) rendered on top. The sections are through the voxel whose activity is depicted in the previous two figures. The SPM{Z} has been thresholded at p = 0.01 (uncorrected).

Limitations

The key limitation of the method described in this paper relates to the constraints and validity of the statistical model used. One key aspect of this model is that it has very few parameters (meaning many more degrees of freedom remain for statistical inference) and that it uses linear approximations. These linear approximations will fail if the observed kinetics deviate substantially from the canonical forms assumed for the tracer. This is a relatively safe limitation in the sense that a failure to model the real kinetics also implies a failure to model changes. This dual failure protects against false positives in regions where the model is inappropriate: the statistical model can be thought of as “focusing” on gray matter because the canonical curves assumed are typical of gray matter. This is an intuitive conjecture that will require a more formal analysis but, based on anecdotal observations, the conjecture seems reasonable.

A second limitation of the approach described here relates again to the validity of the model employed. In particular, a poor model may lead to bias. For example, the bias in Z scores observed in the validation study probably reflects the failure of a single-compartment model to fit the data properly from multicompartmental kinetics (either real or apparent, secondary to smoothing). Other explanations for this failure may include errors in the input function or potential interactions of the challenge with arterial measures. Many of these concerns would be addressed by using more elaborate kinetic models or by using spectral analysis to identify the number of compartments required. These and other ideas are the subject of current work (Vin Cunningham and Roger Gunn, personal communication).

Another limitation is that the form of the dynamic displacement changes is fixed and specified at the outset, thus precluding the possibility of regionally specific differences in the form of time-dependent displacement (e.g., phasic vs. sustained). This limitation and others are addressed in the next section on possible extensions of the current approach.

We have not presented any discussion on the timing of the acquisition frames in relation to optimising sensitivity to changes. Clearly, the challenge should be sufficiently late to render the evoked changes in kinetics relatively distinct from the kinetics themselves and yet not be so late that the changes occur in frames with poor counts. We do not address this issue here but acknowledge its importance for future work (see also Morris et al., 1995).

In this paper, we have focused on some of the methodological issues that arise in dynamic displacement studies. We acknowledge that many other pharmacokinetic issues will need to be addressed in the application of these sorts of techniques: For example, we have assumed that the changes in tracer clearance are due to interactions at the level of the receptor; however, other mechanisms may also contribute to the observed changes (e.g., changes in plasma protein binding of the tracer, peripheral metabolism of the tracer, tracer delivery, or changes in nonspecific binding).

Multicompartment models and complicated forms of time-dependent changes

In principle, the mathematical development described in the theory section can be extended easily to multicompartment models as described in the general approach section. Similarly, the expansion of the time-dependent effect on kinetic parameters [Γ _i (t)] can be in terms of many basis functions (e.g., polynomial, Fourier, and others), which would provide for a greater flexibility in the sense that maintained, phasic, or even biphasic time-dependent changes in displacement could be modelled. In this instance, the F ratio would be used to infer that some linear combination of the basis functions was sufficient to reduce the error variance significantly and the t statistic to test for specific linear combinations of these effects using a contrast with zero mean. Although we will be evaluating these extensions, it should be noted that the “minimal” analysis described above could be repeated for different forms of time-dependent changes in displacement [i.e., the shape of h(t) in Fig. 2 (broken line) could be changed to test for a more acute or even biphasic effect], resulting in a series of SPM{Z}, each testing a specific hypothesis about the form of dynamic displacement. Given that this facility exists, the usefulness of less constrained (and possibly less powerful) statistical models requires careful evaluation.