Recent developments in optimal experimental designs for functional magnetic resonance imaging

Abstract

Functional magnetic resonance imaging (fMRI) is one of the leading brain mapping technologies for studying brain activity in response to mental stimuli. For neuroimaging studies utilizing this pioneering technology, there is a great demand of high-quality experimental designs that help to collect informative data to make precise and valid inference about brain functions. This paper provides a survey on recent developments in experimental designs for fMRI studies. We briefly introduce some analytical and computational tools for obtaining good designs based on a specified design selection criterion. Research results about some commonly considered designs such as blocked designs, and m-sequences are also discussed. Moreover, we present a recently proposed new type of fMRI designs that can be constructed using a certain type of Hadamard matrices. Under certain assumptions, these designs can be shown to be statistically optimal. Some future research directions in design of fMRI experiments are also discussed.

Key Words: A-optimality; Blocked designs; Design efficiencies; D-optimality; Genetic algorithms; Hadamard sequences; M-sequences

Core tip: This paper provides an overview on recent developments in the design of functional magnetic resonance imaging experiments (fMRI). We discuss both analytical results and computational approaches that are currently available for selecting high-quality fMRI designs.

INTRODUCTION

Recent years have seen an upsurge of functional brain imaging experiments for a better understanding of how humans learn, remember and make decisions. Such experiments are also widely conducted by researchers to help provide paths to treat/prevent some terrifying brain disorders such as Alzheimer’s disease, and are thus very valuable. As in many scientific investigations, designing a high-quality experiment is an important first step for successful functional brain imaging studies. A carefully designed experiment allows experimenters to collect informative data to make precise inference on the goals/hypotheses at minimal cost. On the other extreme, data collected from a poorly designed experiment may fail to provide valid answers to the research questions of interest, resulting in a waste of resource. The importance of the use of a carefully selected experimental design (or data collection plan) cannot be overemphasized.

This paper provides a survey on some recent developments in experimental designs for functional magnetic resonance imaging (fMRI) experiments. Functional MRI is one of the most common functional brain mapping technologies. This pioneering, noninvasive technology helps to study experimental subjects’ brain activity when they are cognitively engaging with mental stimuli such as viewing pictures, tapping fingers, solving problems, recalling events, or making decisions. It is used in various research areas including psychology, economics, and cognitive neuroscience[1], and has great clinical potentials as highlighted in a special issue on clinical applications of fMRI in Neuropsychology Review, Vol. 7, No. 2, 2007. However, fMRI experiments are usually expensive, and the collected data is notoriously noisy, making it difficult to draw precise statistical inference on brain functions. We thus would like a high-quality experimental design to help us make the best use of the limited resources to collect informative fMRI data.

An fMRI design is a sequence of mental stimuli to be presented to an experimental subject in an fMRI experiment. While the subject is performing the tasks determined by the selected stimulus sequence, an MRI scanner repeatedly scans his/her brain to acquire fMRI data for making statistical inference about the brain activity. The quality of the collected data depends on the selected design. However, due to the complexity of fMRI, obtaining the “best” fMRI design suited to the goal(s) of the experiment is a challenging task. We usually need to consider not only the statistical efficiency in achieving one or more (competing) study objectives, but also some unwanted psychological effects that can contaminate the data. In addition, we may want the obtained design to fulfill some practical constraints. The large diversity of the fMRI experimental settings and protocols also contributes to the difficulty of design selection. In almost all cases, we deal with a very challenging combinatorial problem.

There are some advances in the selection of fMRI designs, but much more work is needed to move this new emerging research area forward. The purpose of this article is to provide a brief overview of stochastic and deterministic computational tools for designing efficient fMRI studies as well as recent insights obtained for such studies using analytical methods. We begin in the next section with background information on fMRI studies, and introduce terminology and notation used in this article. We then present the general linear models widely used for the design and analysis of fMRI studies and popular design criteria in this area. Some recently obtained results and guidelines for selecting fMRI designs are discussed. We close the article with a summary and discussion.

BACKGROUND

Terminology and notation

In a typical fMRI experiment, a sequence of mental stimuli (e.g., pictures) of one or more types interlaced with periods of rest or, say, visual fixation is presented to each experimental subject. These stimuli give rise to neuronal activity at some brain regions that triggers an increased inflow of oxygenated blood, leading to a decrease in the concentration of deoxygenated blood. This change in the ratio of oxy- to deoxy-blood can influence the strength of the magnetic field, and results in a rise and fall in the intensity of signals collected by the MRI scanner. Specifically, the MRI scanner collects MRI measurements by repeatedly scanning each of the, say, 64 × 64 × 30 brain voxels, which are volumetric image elements that cover (part of) the subject’s brain. Some voxels may fall outside the brain; see also Subsection 2.1.1 of Lazar[2]. At each voxel, MRI measurements are collected every τ_TR (e.g., 2) seconds to form a blood oxygenation level dependent fMRI time series. The pre-specified time τ_TR is called the time to repetition. These time series serve as surrogate measurements of the underlying neuronal activity, and are analyzed to make inference about how the brain reacts to the stimuli; see also Lazar[2].

The inference on brain activity is mainly based on some characteristics of the hemodynamic impulse response function (HRF). The HRF is a function of time describing the rise and fall of the noise-free MRI measurements following a brief neuronal firing that occurs at a voxel. Previous studies suggest that the HRF may increase from baseline in about two seconds after the onset of a brief stimulus, reach the peak in five to eight seconds, and possibly fall down below baseline before its complete return to baseline[1,3]. This process may take about 30 s, counting from the onset of the brief stimulus to the HRF’s complete return to baseline. If there are other neuronal firings (e.g., due to the onset of other stimuli) before the cessation of the previous HRF, the evoked HRFs overlap and their heights accumulate. Since fMRI time series is typically very noisy, identifying the characteristics of the HRF by visual inspections is difficult, if not impossible. Statistical methods are thus needed to help extract useful information from the data. As an integral part of the statistical process, we would like to select a “good” fMRI design that helps to make valid inference.

An fMRI design is a sequence of mental stimuli of one or more types. When the sequence is presented to an experimental subject, each stimulus may last as brief as several milliseconds or as long as, say, a minute. Stimuli with extended presentation duration, e.g., 10-60 s, are used in traditional blocked designs, which are also termed as boxcar designs. In such a design, the stimulus of the same type can appear at multiple time points during the experiment, but each long stimulus is immediately followed by a long stimulus of another type or by a period of control (e.g., rest). It also is not uncommon to replace each long stimulus by a short sequence of separate but brief stimuli of the the same type. The resulting designs are still called blocked designs. For experiments with Q stimulus types, a typical blocked design may be the repetitions of {A₁A₂…A_QA₀}, where, for q =1, ..., Q, A_q represents a presentation of a long stimulus (or a sequence of brief stimuli) of the q^th type, and A₀ is a period of control. At a brain voxel responding to the q^th-type stimulus, neuronal firings can be expected throughout the time span of each “on-period” A_q. This leads to an accumulation of overlapping HRFs. With a long on-period of the stimulus, the MRI signal intensity increases to a high level, and may reach a plateau before dropping down to baseline following the cessation of the stimulus. The large contrast between the elevated signal intensity and baseline facilitates the detection of brain voxels (or regions) that respond to the stimulus. Blocked designs are thus often recommended for detecting brain voxels that are activated by the stimuli; see the Results on Design Selection section for a further discussion.

Moving away from blocked designs, some studies showed that an individual stimulus that is as brief as several tens of milliseconds can evoke a detectable change in the MRI measurements; see Rosen et al[3] and references therein. In addition, the heights of overlapping HRFs following multiple brief stimuli tend to be (roughly) additive when the time between stimulus onsets is not overly short (e.g., at least 2 s); see also Friston et al[4]. These observations make it possible to consider event-related (ER-) fMRI designs that consist of brief stimuli whose order may be randomized. An ER-fMRI design of Q stimulus types is often written as a finite sequence of elements 0, 1, ..., Q, and may look like d = (1012021…1). A positive integer q in d represents an onset of a q^th-type stimulus, and 0 means no stimulus onset. Specifically, when the i^th element of d is d_i = q (> 0), a q^th-type stimulus appears briefly at time (i-1)τ_ISI for a pre-determined τ_ISI; time 0 may be synchronized to the first valid MRI scan. For example, when d₃ = 1 and τ_ISI = 4 s, a stimulus of the first type (e.g., a picture of a familiar face) will occur briefly at the (3-1)τ_ISI = 8^th second after the first valid MRI scan. With d₄ = 2, a stimulus of the second type (e.g., a picture of an unfamiliar face) will appear at the 12^th second after the first valid MRI scan. When d_i = 0, there is no stimulus onset at time (i-1)τ_ISI. With these 0’s in the design, time between stimulus onsets may be “jittered”[5], and thus, may not be fixed to τ_ISI. Typically, the control (e.g., a visual fixation or rest period) fills in the time between the offset of a brief stimulus to the onset of the next stimulus. Due to its flexibility, ER-fMRI designs have gained much popularity[6]. However, a typical design can easily contain tens or hundreds of elements, making it very challenging for selecting good designs. In this paper, we discuss some recently developed approaches for finding high-quality fMRI designs, including both blocked and ER-fMRI designs. Most of these approaches are built upon the popular general linear model framework. This framework is described below.

The general linear model framework

The fMRI time series, {y(t): t ≥ 0}, of a brain voxel is typically modeled as the sum of (1) the convolution of the stimulus function and the HRF, (2) a nuisance term allowing for a trend or drift of y(t); and (3) noise[7,8]. We consider the foxllowing continuous-time model:

y (t) = ∑^Q_{q = 1}∫^t₀χ_q (t-τ)h_q (τ; β_q)dτ + s(t; γ) + e(t) (1)

where x_q(t) is the stimulus function for the stimuli of the q^th type, h_q(τ; β_q) is the HRF evoked by the q^th-type stimulus, β_q is an unknown parameter vector, q = 1, ..., Q, s(t; γ) is a nuisance term approximating the drift/trend of the time series, γ is the corresponding unknown parameter vector, and e(t) is noise. The stimulus function x_q(t) indicates the appearances of the q^th-type stimuli, and may be a sum of boxcar functions or a sum of (shifted) Dirac delta functions; see also Henson and Friston[9]. Boxcar functions are often employed in experiments with blocked designs. In this case, x_q(t) takes a positive value during the “on-periods” of the q^th-type stimulus of a block design, and is 0, otherwise. The resulting model is sometimes referred to as the epoch model[10]. In an event-related model, the x_q(t) is a sum of (shifted) Dirac delta functions that indicates the onset times of the brief stimuli of the q^th type.

The most commonly used fMRI data analysis method is probably the general linear model approach[1]. Partly due to this popularity, existing studies on fMRI designs mainly focus on linear models such as models (2) and (3) that are extensions of (1), and are linear in the parameters β_q’s and γ. In the fMRI literature[11,12], dual models are commonly considered for two popular study objectives, namely the detection of brain activations (or detection) and the estimation of the HRF (or estimation). The main difference between the two models is that they used different sets of basis functions to describe h_q(τ; β_q) of model (1); see also Friston et al[4].

For detection, the HRF h_q(τ; β_q) is typically approximated by θ_qh^*(τ), where h^*(τ) is an assumed shape of the HRF, and θ_q is the unknown amplitude (or maximum height) of the HRF. Thus, β_q contains only one parameter θ_q that signals the strength of brain activation due to the q^th-type stimulus. Since the MRI measurements y(t) is collected every τ_TR seconds, we consider the following discrete-time model:

y = ∑^Q_{q = 1} z_qtheta;_q + Sγ + e. (2)

Here, y=(y₁, ...., y_T)’ with y_t = y((t-1)τ_TR). The vector z_q is obtained by subsampling the convolution of x_q(t) and h^*(τ) with a sampling rate of τ_TR seconds. Sγ corresponds to s(t; γ) of (1) with S being a specified matrix. For example, the t^th element of Sγ might be γ₀+γ₁t+γ₂t². The vector e in model (2) represents the noise. The focus of model (2) is typically on C₁θ for a given matrix C₁ whose rows contain coefficients of linear combinations of θ₁,…, θ_Q; here, θ = (θ₁,…, θ_Q)’. When C₁= I_Q is the identity matrix of order Q, the focus is on the strength of brain activation due to each stimulus type. It is also common to study (θ_p - θ_q) for p ≠ q. In such a case, the rows of C₁ contain the coefficients of the pairwise comparisons between the HRF amplitudes.

The estimation of the HRF is a study objective that has gained much popularity with the advent of ER-fMRI. A widely used model for this objective is:

y = ∑^Q_{q = 1} X_qh_q + Sγ + e. (3)

Here, h_q = (h_1q, ..., h_Kq)’ is an unknown parameter vector representing the heights of the HRF that contribute to y. Specifically, h_kq= h_q((k-1) ∆T; β_q) is the HRF height at (k-1)∆T seconds after the onset of a q^th-type stimulus, where ∆T is the greatest real value making (τ_ISI/∆T) and (τ_TR/∆T) integers; k=1, ...,K. The value of K is selected so that h_q((k-1) ∆T; β_q) becomes negligible when τ > (K-1)∆T. For a commonly considered 32-second HRF, K = [32/∆T] with [a] being the integer part of a. The consideration of h_q is equivalent to modeling h_q((k-1) ∆T; β_q) of (1) by a linear combination of K shifted Kronecker delta functions; i.e., h_q((k-1) ∆T; β_q) = ∑^K_{k = 1} h_kqδ_k(τ), where δ_k(τ) = 1 when τ = (k-1)∆T, and δ_k(τ) = 0, otherwise; β_q thus contain all the K coefficients (HRF heights) h_1q, ..., h_Kq. X_q = [x_1q, ..., x_Kq] in model (3) is the 0-1 design matrix of size T-by-K for the q^th-type stimuli. The t^th element of x_kq is 1 when h_kq contributes to y_t. The remaining terms in (3) are as in (2). In contrast to model (2) for detection, model (3) does not assume a known shape for the HRF. The goal is to estimate all the unknown HRF heights h_kq or to study some linear combinations C₂h of these heights with h = (h₁’, ..., h_Q’)’ and a given linear combination coefficient matrix C₂.

Design selection criteria

With models (2) and (3) respectively for detection and estimation, the main design goal is to select an fMRI design that yields the most precise parameter estimates of the parametric functions of interest. Some statistically meaningful optimality criteria have been proposed for evaluating the goodness of competing designs. Two popular criteria in the fMRI literature are A- and D-optimality criteria. For detection problems with model (2), the A-optimality criterion can be defined as the following ‘larger-the better’ criterion:

φ^A_θ (d) = r₁/trace {C₁[Z’V’(I_T-ω{VS}VZ)]-C^’₁ } = r₁/trace {C₁[M₁(d)]-C^’₁}’ (4)

Here, r₁ is the number of rows of C₁, and Z = [z₁, ..., z_Q]. V is a whitening matrix such that cov(Ve) = σ² VRV’ = σ²I_T, where σ² is the error variance, and R = corr(e) is the correlation matrix of errors. The matrix ω{A} = A (A'A)^-A’ is the orthogonal projection matrix onto the column space of A. A^- is a generalized inverse of A, and M₁(d) = Z’V’(I_T - ω{VS}VZ is the information matrix of theta;. We note that V may be obtained by, e.g., the Cholesky decomposition of R^-1, and, depending on the assumptions made at the design stage, it may or may not contain unknown parameters; see also the Results on Design Selection section and Maus et al[13]. The criterion in (4) depends on the selected design d through the design matrix Z, and is inversely proportional to the average variance of the least-squares estimates of the parametric functions defined by C₁theta;. For estimating the HRF with model (3), the A-optimality criterion can be written as:

φ^A_h (d) = r₂/trace {C₂[X’V’(I_T-ω{VS}VX)]-C^’₂ } = r₂/trace {C₂[M₂(d)]-C^’₂}’ (5)

where r₂ is the number of rows of C₂, X=[X₁, ..., X_Q] is the design matrix depending on the selected design d, M₂(d) is the information matrix for h, and all the remaining terms are as in (4).

The D-optimality criterion seeks to minimize the volume of the (asymptotic) confidence ellipsoid of C₂h. For the detection of brain activation with model (2) and the estimation of the HRF with model (3), D-optimal designs are found by maximizing the following two criteria, respectively:

φ^D_θ (d) = det{C₁[M₁(d)]-C^’₁}^-1/r1 (6)

φ^D_h (d) = det{C₂[M₂(d)]-C^’₂}^-1/r2 (7)

All the terms in (6) and (7) are as in (4) and (5), respectively. For the D-optimality criteria, the coefficient matrices C₁ and C₂ are required to be full row rank. The selection between the A- and D-optimality criteria depends on the need and preference of the experimenter. As indicated in Maus et al[14], while early works on fMRI designs mainly focused on the A-optimality criterion, there is no obvious reason to generally prefer one criterion over the other. In the subsequent sections, we discuss some results on fMRI design selection. Most of these results are based on the A- or D-optimality criterion.

RESULTS ON DESIGN SELECTION

Blocked designs for detecting brain activations

There is some guidance on selecting blocked designs for detecting brain activations in the literature. For example, Henson[15] advocated the use of blocked designs having a 15-s-on-15-s-off pattern. For such a blocked design formed by {A₁A₂…A_QA₀}, the duration of each A_q is fixed to 15 s. This suggestion is based on the Fourier transformations of the convolution in (1) by assuming that the HRF has the form of the double-gamma function:

g*(τ) = τ^se^-τ/S! - 1/6 ×τ^1se^-τ/1S! (8)

The double-gamma function is widely used as the HRF shape, and is built in a software package, called SPM (http://www.fil.ion.ucl.ac.uk/spm/), for fMRI data analysis. In the frequency domain, this HRF acts as a low-pass filter that `passes’ low-frequency signals and reduces the amplitude of high-frequency signals. As demonstrated in Henson[15], after the Fourier transformation, a large proportion of the signal energy of a 15-s-on-15-s-off blocked design is retained by the selected HRF shape. In addition, the use of an on-period A_q that is longer than 50 seconds is not recommended. This is because the signal energy of the resulting blocked designs may be lost after accounting for the low-frequency nuisance signals such as heartbeats or respirations which is modeled by s(t; γ) in (1).

Setting the block length (or duration of A_q) to 15 s may not be optimal for an HRF shape that is different from (8). For example, Liu et al[16] considered cases with one stimulus type (Q = 1), and evaluated the performance of designs with the A-optimality criterion. They observed that the blocked design with a block length of 64 s tends to have a high statistical efficiency in detection when a single gamma density function is used to model the HRF shape. This latter HRF shape is also not uncommon, especially for cases where the HRF does not fall below baseline when returning from its peak. In addition, Liu et al[16] suggested that the selection of block length also depends on s(t; γ). In particular, they demonstrated that the blocked design with a 64-s block length can yield a smaller φ^A_θ -value than designs with a shorter (e.g., 32 s) blocked length when the statistical model also allows for a second- or third-order Legendre polynomial drift.

To provide additional information on design selection, Maus et al[14] studied blocked designs of two stimulus types (Q = 2) with selected block lengths (10, 15, 20, 30 or 60 s), and patterns (repetitions of {A₁A₂ }, {A₁A₂A₀}, or {A₁A₀A₂A₀}). Each block A_q is formed by a sequence of 1-second stimuli of the q^th type; q = 1, 2, and the time between the onsets of consecutive stimuli in the same block is τ_ISI = 1, 2, or 3 s. They compare the statistical efficiencies of these blocked designs in detecting brain activations via model (2). In their model, the nuisance term Sγ corresponds to a linear trend, and the HRF shape used to construct z_q is set to the double-gamma function of (8). The errors are assumed to have one of the three possible structures, including uncorrelated errors, first order autoregressive (AR1) process, and an AR1 process plus a measurement error (AR1+ME).

Considering both φ^A_θ (d) and φ^D_theta; (d), Maus et al[14] suggested to keep τ_ISI as short as possible. In addition, they recommended to use the design pattern {A₁A₂A₀} for studying the HRF amplitudes θ₁ and θ₂. When the focus is on comparing the amplitudes (i.e., θ₁ - θ₂), blocked designs formed by {A₁A₂} are recommended. The results of Maus et al[14] also indicate that the selection of block length may hinge on the assumed error correlation. When the focus is on θ_q’s, a block length of 15 s is recommended for both uncorrelated and AR1 errors. As for AR1 + ME errors, a block length of 10 s is the best among the selected blocked lengths. For studying the contrast between the HRF amplitudes, the suggested block lengths are 20 s and 15 s for uncorrelated errors and correlated errors (AR1 or AR1 + ME), respectively.

These previous studies provide some guidelines on selecting blocked designs for detecting brain activations. It can also be seen that the selection of blocked designs depend on a few factors. These factors include the parametric function C₁θ of interest, the selected HRF shape, the model for capturing the drift/trend of the fMRI time series, and the error correlation structure. For cases that are not covered by these guidelines, we may obtain a good design for detection by using a computer algorithm. Some algorithms have already been proposed in the fMRI literature. Most of these computational approaches can be employed for cases considering the detection of brain activations, the estimation of the HRF, or when both detection and estimation are of interest. Some practical constraints may also be imposed when using these computational tools. In what follows, we first describe some guidelines for selecting ER-fMRI designs for estimating the HRF. We then discuss computer algorithms for obtaining good fMRI designs.

ER-fMRI designs for estimating the HRF

The estimation of the HRF helps to make inference about some characteristics of the underlying neuronal activity as also described in Lindquist et al[17]. For this objective, model (3) may be considered, and the goal is to obtain a design yielding the most precise parameter estimates of C₂h for a given C₂. By considering the A-optimality criterion of (5), Dale[18] suggested to allow for variable time intervals between onsets of consecutive stimuli, and the average of these time intervals should be kept small. This suggestion can also be applied to the D-optimality criterion of (7). However, one should take caution that if the time between stimulus onsets is overly short (e.g., < 2 s), the accumulated heights of the overlapping HRFs may saturate at a certain level. Consequently, the assumption of the additivity of the HRF heights can be violated. For such a case, the nonadditive HRF heights should be taken into account when evaluating the goodness of designs; see also, Wager et al[19] and Wager et al[20]. However, current methods for accounting for the nonadditive HRF heights tend to be ad hoc, and additional investigations are needed.

While rendering useful information, Dale[18] did not provide a systematic way for design construction. Buračas and Boynton[21] worked on the same design issue, and advocated the use of maximum length shift-register sequences (or m-sequences). Such a design can be generated by a primitive polynomial over a Galois field GF(Q+1) consisting of Q+1 elements, where Q+1 is a prime power. To construct an m-sequence, one may select a primitive polynomial f(x) = x^r - ∑^r_{i = 1}α_ix^r-i from, e.g., Table 3.5, 3.6 or 3.7 of Golomb and Gong[22]. The m-sequence d = (d₁, ..., d_N) is then determined by the relation, d_{n + r} = ∑^r_{i = 1}α_id_{n + r - j} (mod Q+1) with a nonzero initial r-tuple (d₁, ..., d_r); see also Lidl and Niederreiter[23], and MacWilliams and Sloane[24]. Such a design can also be obtained via an MATLAB program developed by Liu[11]. For an m-sequence of length N = (Q+1)^r - 1, every nonzero r-tuple appears exactly once in the set {(d₁, ..., d_r), (d₂, ..., d_r+1), ..., (d_N,d₁ ..., d_r-1)}.

Buračas et al[21] and Liu[11] reported the high performance of m-sequences in terms of the φ^A_h -value when C₂ = I_QK is the QK-by-QK identity matrix. However, when Q > 1, the frequency of the appearance of each stimulus type of an m-sequence can be different from the optimal stimulus frequency approximated by Liu et al[12] for A-optimality. In particular, Liu and Frank[12] indicated that the optimal stimulus frequency of an A-optimal design for estimating the HRF h is about 1/(Q + √Q) for each of the Q stimulus types. The optimal number of 0 is thus approximately N/(1 + √Q). Since the stimulus frequency of m-sequences is about 1/(Q+1), these designs may not be A-optimal; see also Kao et al[25]. For φ^D_h with C₂ = I_QK, the optimal stimulus frequency approximated by Maus et al[13] is 1/(Q+1), and is close to that of m-sequences. Maus et al[13] thus suggested that the optimality of m-sequences may depend on the selected criterion.

However, attaining the (approximated) optimal stimulus frequency does not guarantee an optimal design. To derive additional insightful results, Kao[26] also studied model (3) with the following assumption: Assumption 1. (a) The number of MRI scans T equals the length N of the design d and τ_TR = τ_ISI; (b) S = j_T is the T-by-1 vector of ones, and cov(e) proportional I_T; and (c) the last K - 1 elements of design d are also presented to the subject before the first valid MRI scan.

Assumptions 1(a) and 1(c) are mild, and can often be controlled by the experimenters. Assumption 1(b) is mainly for mathematical simplicity, and is also considered in some previous studies such as Liu et al[16] and Maus et al[27]. Following an argument in Kushner[28], for the results to be discussed in the remaining of this subsection, Assumption 1(b) can be relaxed to include cases with cov(e) = αI_T + λj_T’ + j_Tλ', where α is a constant and λ is a vector of constants. The results thus hold for a compound symmetric covariance matrix with cov(e) = αI_T + λJ_T, where λ is a constant, and J_T is the T-by-T matrix of ones. For estimating the K-by-1 HRF parameter vector h₁ with one stimulus type (Q = 1), Kao[26] showed that a design of length N having n₁ = N/2 and n_r⁽¹¹⁾ = (n₁)²/N for all r = 1, ..., K - 1 is universally optimal. Here, n_q is the frequency of the q^th-type stimuli in the design d, and n_r^(pq) is the number of times (d_n-r, d_n) = (q,p) for n = 1, ..., N; d_n-r = d_N+n-r when n ≤ r. We also note that an universally optimal design can be shown to be optimal in a large class of optimality criteria, including A- and D-optimality[29]. For Q > 1, a similar sufficient condition for an ER-fMRI design to be D-optimal can also be found in Kao[26]. In particular, if all the symbols 0, 1, ..., Q appear equally often in a design d of length N, and that n_r^(pq) = n_pn_q/N for all p, q = 1, ..., Q and r = 1, ..., K - 1, then the design d maximizes φ^D_h of (7) under Assumption 1 and C₂ = I_QK.

As described in Kao[26], designs satisfying the previously mentioned sufficient conditions can be constructed by inserting an additional 0 to any (K - 1)-tuple of zeros in an m-sequence of length (Q+1)^K - 1. The resulting design is a de Bruijn sequence[22,30]. Aguirre et al[30] proposed to use de Bruijn sequences for estimating the HRF. The results of Kao[26] help to establish the optimality of such designs.

Clearly, m-sequences do not satisfy the sufficient conditions provided by Kao[26]. Additional results are thus needed for establishing the optimality of these popular designs. Kao[31] worked on this direction, and proved that a binary m-sequence of length N ≥ 2K - 3 is D-optimal for estimating the HRF h₁ under Assumption 1 with Q = 1. He also proposed a new type of ER-fMRI designs for estimating the HRF. This new type of designs, which are termed as Hadamard sequences, can be constructed by a normalized Hadamard matrix, H, having a circulant core. Specifically, the elements of the first row and column of H are 1 and all the other entries are +1 or -1 with HH’ proportional I. After deleting the first row and column of H, we have a circulant matrix called the circulant core. As described in Kao[31], a D-optimal design can be achieved by replacing +1 and -1 in any column of the circulant core by 0 and 1, respectively. It is noteworthy that binary m-sequences can also be generated using this same method and are thus special cases of Hadamard sequences. Nevertheless, Hadamard sequences exist in many different lengths for which a binary m-sequence is unavailable. These newly proposed designs are thus much more flexible than m-sequences and the previously mentioned de Bruijn sequences in terms of design length.

Kao[31] also conducted some case studies on the performance of Hadamard sequences when Assumptions 1(b) and 1(c) are violated. Based on empirical results, Hadamard sequences tend to remain efficient when the nuisance term Sγ in model (3) corresponds to a second-order polynomial drift, the noise follows an AR1 process, and/or no stimulus is presented before the first valid MRI scan. This result is especially true when the autocorrelation coefficient of the AR1 noise is not as high as ρ = 0.5 or when the design is not too short (e.g., N < 100). We also note that a violation of Assumption 1(a) can have a great impact on the performance of Hadamard sequences. For cases with τ_TR≠τ_ISI, we may consider efficient computational methods for obtaining good designs. Some computational approaches are introduced in the next subsection. These approaches are also applicable when both estimation and detection are of interest.

Computational tools for obtaining fMRI designs

In the fMRI literature, some computer algorithms are proposed for finding an ER-fMRI design of the form d = (d₁, ..., d_N) with d_n belong to {0, 1, ..., Q} that optimizes a specific single- or multi-objective optimality criterion. To efficiently search over the enormous space of ER-fMRI designs for good designs, Wager and Nichols[19] advocated the use of the genetic algorithm (GA) technique. Due to their versatility, GAs can accommodate various experimental settings to find designs suited to individual fMRI experiments. Following Wager and Nichols[19], Kao et al[25] put forward an efficient GA that takes advantage of knowledge on the performance of some ER-fMRI designs to improve the efficiency of the GA search. Some well-known designs such as m-sequences, blocked designs, and their combinations are employed in the algorithm of Kao et al[25] to increase the diversity of the designs being explored, and to maintain a supply of good traits (or building blocks) that help to form good designs during the GA search. As demonstrated in Kao et al[25], this strategy is very effective.

With the previously mentioned GAs, one can find a (near-)optimal design for user-specified number of stimulus types Q, design length N, τ_ISI, τ_TR, and model assumptions, including the model for drift/trend of the time series, error correlation structure, and, if model (2) is considered, the HRF shape. Depending on the study objective(s), the optimality criterion for evaluating the quality of designs may be φ^A_θ , φ^A_h , φ^D_θ , φ^D_h or a weighted sum of some of these criteria; weights are user-selected to reflect the relative importance of detection and estimation. In a weighted sum criterion, one may also include other individual criteria to account for quantifiable constraints/requirements of the study. For example, Wager et al[19] included a counterbalancing criterion for avoiding psychological confounds such as anticipation and habituation. By optimizing this criterion, the order of the stimuli in the resulting design cannot be easily predicted by the experimental subject. Moreover, we may include an additional individual criterion to measure the departure from a target frequency of appearances of each stimulus type; see also, Kao et al[25]. Such a customized requirement on the stimulus frequency may help to increase the subject’s engagement in the presented mental tasks[32].

The GA of Kao et al[25] has been applied for studying several fMRI design issues. For example, this algorithm was used to obtain designs for cases where both individual stimulus effects (h and θ) and pairwise comparisons (h_p - h_q and θ_p - θ_q for p ≠ q) are of interest. Maus et al[13] used the GA to work on cases where the autocorrelation coefficient ρ of the AR1 noise is uncertain. The GA is also adapted in Kao et al[33,34] for finding designs suited to experiments with multiple scanning sessions.

In addition, Maus et al[35] and Kao et al[36] utilized the GA to tackle the design problem concerning an uncertain HRF shape. The need for considering the uncertainty of the HRF shape is manifested in some previous studies[37,38]. These studies pointed out that the HRF shape may vary across brain voxels, and that specifying a wrong HRF shape in, say, model (2) for detection may lead to an incorrect conclusion. To accommodate different HRF shapes, Kao[39] considered at the design stage the following nonlinear model:

y = ∑^Q_{q = 1} X_q h(u)θ_q + sγ + e (9)

where h(u) is a K-by-1 vector representing the shape of the HRF, u is an unknown parameter vector that needs to be estimated from data, and all the remaining terms are as in (2) and (3). The vector h(u) may be determined by the double-gamma function of (8) with free parameters for accounting for the variability in the HRF shape; see also Wager et al[20]. In particular, the k^th element of h(u) is g((k - 1) ∆T; u)/max_s g(s; u) with u = (u₁, u₂)’ and

g(τ; u) = [(τ-u₂)^{u1-1e-(τ-u2)}]/Γ(u₁) - 1/6 × [(τ-u₂)^{1S-1e-(τ-u2)}]/15! (τ≥ u) or 0 (otherwise) (10)

Here, u₁ is the time-to-peak parameter, which mainly determines the time for the HRF to reach the peak, counting from its onset time. The time-to-onset parameter u₂ determines the time when the HRF starts to increase from baseline, counting from the onset of a stimulus. As indicated by Wager et al[20], these two parameters are the most influential, although some additional free parameters may also be included in (10). For example, one may use a free parameter to replace the coefficient 1/6 in the second term of the non-zero part of (10). The function Γ(u) = ∫^∞₀ t^u-1e^-1dt = (u-1) Γ(u - 1) in (10) is the gamma function. We note that the function g^*(τ) in (8) is a special case of (10) with u = (6, 0)’. Specifically, the HRF shape h^*(τ) in model (2) depends on g^*(τ), and is fixed. By contrast, the HRF shape in model (9) is determined by g(τ; u), and involves unknown parameters to be estimated from the data. The latter model is thus more flexible.

When making inference about θq for detecting brain activations, model (9) allows for an uncertain HRF shape. However, obtaining a good design for such a flexible model is quite challenging. Again, we would like a design optimizing some function (e.g., the A- or D-optimality criterion) of the information matrix of θ. For model (9), this information matrix, denoted by M(d; θ, u), can be approximated by first-order Taylor approximation. In contrast to M1(d) and M2(d) in (4)-(7), M(d; θ, u) depends not only on the design d, but also on the unknown model parameters θ and u; see Kao[39] and Kao et al[36] for details. By treating θ and u as random variables, and assuming the availability of a (prior) distribution of θ and u, Kao[39] targeted a (pseudo-)Bayesian design that maximizes E{(M(d; θ, u))} for a larger-the-better criterion , where the expectation E{.} is taken over the (prior) distribution of the parameters.

When a prior distribution of the parameters is unavailable, it is common to consider to maximize the minimum of (M(d; θ, u)), where the minimum is taken over the possible values of θ and u. It also is popular to maximize the the minimum of the relative efficiency, which is defined as

min (θ ∈ u, u ∈ U) [M(d; theta;, u)]/[M(d^*_θ,u; θ, u)]

Here, Θ and U contain the possible values for theta; and u, respectively; and d^*_{θ, u} is a locally optimal design that maximizes (M(d; θ, u)) for given θ and u. Designs maximizing the former criterion are termed as maximin designs, whereas those optimizing the latter criterion are maximin-efficient designs. Both criteria are popular in the literature; see also Kao et al[36] and references therein. However, obtaining maximin-type designs is computationally very expensive. Kao et al[36] proposed an efficient shortcut. Building on some analytical results, they showed that the size of the parameter space of Θ can be greatly reduced when obtaining maximin-type designs. Specifically, when Q=1, we may find a very efficient maximin (or maximin-efficient) design by focusing on θ₁ = 1 (or θ₁ belong to {0,1}). For Q > 1, instead of setting Θ to the entire Q-dimensional space, we may focus on a subspace consisting of (1/Q!) of the surface of the Q-dimensional unit hemisphere centered at the origin when obtaining a maximin design; the origin needs to be included in the subspace for finding a maximin-efficient design. To further reduce computing time, Kao et al[36] focused on a restricted class, Ξ₀, of designs when using a search algorithm to find maximin-type designs. Specifically, each design of length N in Ξ₀ is formed by a short design of length [N/Q], where [a] is the smallest integer greater than or equal to a. For any short design, a full-length design is constructed by cyclically permuting the labels of the Q stimulus types with 0's staying intact, and then leaving out the excess elements, if any. The stimulus frequencies in the resulting design are thus (nearly) equal across stimulus types. Kao et al[36] showed that their approach is quite efficient and effective when obtaining maximin-type designs when the HRF shape is uncertain.

In addition to the GA technique, a deterministic optimization algorithm for obtaining optimal fMRI designs has recently been proposed and studied by Kao and Mittelmann[40]. Without stochastic explorations, this latter approach has been demonstrated to be efficient for some cases for which the GA requires much CPU time in finding a good design. The main idea is to combine a greedy hill-climbing algorithm with the previously mentioned cyclic permutation method for constructing designs of Ξ₀. In particular, the algorithm first systematically perturbs a small fraction (e.g., the first four elements) of a short design d_s of length [N/Q] to create some neighboring short designs that are close to d_s in terms of Hamming distance. The search then moves to the neighboring short design d_s that yields the best full-length design via the cyclic permutation method. After this movement, the algorithm continues to work on perturbing another small fraction (e.g., the fifth to eighth elements) of d_s. This process is repeated until no improvement can be achieved. Based on our experience, this approach tends to lead to very efficient designs with greatly reduced CPU time, although the obtained design might not be optimal. Kao and Mittelmann[40] demonstrated the usefulness of their algorithm by finding maximin designs that are robust to mis-specified error autocorrelation coefficients when stationary AR2 errors are assumed. For this case, the GA approach can be very challenging in terms of CPU time.

The algorithms described so far are used to optimize a single objective function. For experiments with two or more study objectives, these previous studies mainly considered weighted-sum criteria that are convex combinations of all the individual criteria of interest. However, selecting appropriate weights for such a weighted-sum criterion might be challenging for some cases, and the assigned weights may not guarantee a satisfactory design. For example, assigning equal weights does not always lead to a design with equal relative efficiency across all the study objectives of interest. To address this fMRI design issue, Kao et al[41] proposed a multi-objective optimization algorithm by modifying the nondominated sorting GA II (NSGA II) of Deb et al[42]. With a single run of the algorithm, the experimenter can obtain not one, but a class of diverse designs for approximating the Pareto frontier; a Pareto frontier is formed by the best possible solutions in a multi-objective optimization problem. A design best suited to the needs of the experiment can then be selected from the obtained design class. The algorithm can also be used to find fMRI designs when there is a constraint such as a required stimulus frequency. This algorithm is recommended when weights on the multiple study objectives are hard to determine.

CONCLUSION

Design of fMRI experiments is an exciting research area. Several analytical and computational approaches have been proposed for obtaining designs that attain high efficiencies in terms of certain practically meaningful design selection criterion. As demonstrated in Jansma et al[43], among others, fMRI designs with theoretically superior performance are often very useful in real-world experiments. The designs obtained in the previous studies are thus valuable. However, much work remains to be done in this area. As indicated by Lindquist[1] in his recent survey on statistical methods for fMRI studies, “as research hypotheses ultimately become more complicated, the need for more advanced experimental designs will only increase further.”

One possible direction of future research is on developing designs for cases with compound stimuli, each containing two or more components; e.g., each stimulus is formed by a cue followed by a task. To our knowledge, there is no systematic study on this important design issue. In addition, fMRI is also widely considered for studying the functional connectivity between brain regions. High-quality experimental designs for this type of studies are also in a great demand. Moreover, developing powerful computational approaches, and insightful analytical results for optimal fMRI designs should always be helpful. For example, the analytical results described in the Results on Design Selection section are mainly for cases where Assumption 1 holds and C₂ = I_QK. It is also useful to consider the case where C₂ is not the identity matrix when contrasts between the HRFs are of interest. Developing novel, insightful analytical results by relaxing Assumption 1 can also help to move this new research field forward.