Multi-Objective Optimal Experimental Designs for ER-fMRI Using MATLAB

Designs for event-related functional magnetic resonance imaging (ER-fMRI) that help to eﬃciently achieve the statistical goals while taking into account the psychological con-straints and customized requirements are in great demand. This is not only because of the popularity of ER-fMRI but also because of the high cost of ER-fMRI experiments; being able to collect highly informative data is crucial. In this paper, we develop a MATLAB program which can accommodate many user-speciﬁed experimental conditions to eﬃciently ﬁnd ER-fMRI optimal designs.


Introduction
Event-related functional magnetic resonance imaging (ER-fMRI) is one of the leading technologies for studying human brain activity in response to mental tasks or stimuli. Before conducting an ER-fMRI experiment, a sequence of stimuli of one or more types interlaced with the control (rest or fixation) is prepared. This sequence of stimuli is presented to an experimental subject while the MR scanner scans his/her brain every few seconds. The blood oxygenation level dependent (BOLD) time series is collected from each brain voxel (a small region of the brain) for statistical analysis. See Josephs et al. (1997), Rosen et al. (1998), Dale (1999), and Bandettini and Cox (2000) for overviews of ER-fMRI.
One important design problem of ER-fMRI is to find an optimal sequence of the stimuli best suited to the researcher's needs. However, this problem is difficult due to the following reasons. First, the design space, consisting of all possible sequences of stimuli, is enormous. Searching over this space for a good design is hard. In addition, the flexibility of ER-fMRI allows researchers to consider two popular statistical objectives, namely estimation and detection. Estimation refers to the estimation of the hemodynamic response function (HRF), a function of time describing an effect of a single, brief stimulus. Detection is to identify brain regions that are activated by the stimuli. Considering both objectives in one study is not uncommon (see also Wager and Nichols 2003), but this requires good multi-objective designs that efficiently achieve these two competing goals. Moreover, statistics is not the only concern for ER-fMRI. Psychology plays an important, even crucial, role. When a design sequence is patterned or easy to predict, psychological effects such as habituation or anticipation can occur to contaminate the data (Dale 1999). A good design should help to avoid these confounds. Furthermore, customized requirements such as the required number of stimuli for each stimulus type can also arise to make the problem even more complicated. As a result, searching for a good multi-objective design is inevitable. We need well-defined multi-objective design criteria (MOcriteria) for evaluating the quality of competing designs, an efficient search algorithm and a program that implements such an algorithm.
In this paper, we develop a program using MATLAB (The MathWorks, Inc. 2006) for finding multi-objective optimal ER-fMRI designs. Our program implements the approach proposed by Kao et al. (2009a), which includes rigorously formulated models, well-defined MO-criteria and a genetic algorithm (GAs). They incorporate knowledge about the performance of wellknown ER-fMRI designs to increase the effectiveness and efficiency of their approach. As demonstrated in their paper, this approach is more efficient than the previous methods and is flexible enough to accommodate different experimental conditions and assumptions. To make the best use of this approach, our program allows users to specify the experimental conditions based on their needs. The designs that we obtain can help researchers to achieve efficient statistical inference.
The rest of the article is organized as follows. Section 2 reviews the approach proposed by Kao et al. (2009a). Section 3 illustrates our computer codes. An example for using our program is in Section 4, followed by conclusion and discussion in Section 5.

Methodology review
A typical ER-fMRI design can be viewed as an alignment of events, including the stimuli and the control. For convenience, the symbols 0, 1, ..., Q are used to represent the events with 0 indicating the control and i a type-i stimulus, i = 1, ..., Q; Q is the total number of stimulus types. A design, denoted by ξ, looks like ξ = {101201210...1}.
While being presented to an experimental subject, each stimulus lasts for a short period of time relative to the inter-stimulus interval (ISI), the fixed time interval between the onsets of consecutive events. We note that 0s in the sequence are "pseudo-events"; they help to calculate the onset times of stimuli. For example, with a 0 in between, the first, second and the third stimuli (1, 1, and 2) of ξ occur, respectively, 1ISI, 3ISI, and 4ISI seconds after the outset of the experiment. The control fills up the time period between the end of a stimulus and the start of the next one.
Our goal is to find a best sequence of the events to efficiently achieve four popular objectives, which are 1) estimating the HRF, 2) detecting brain activation, 3) avoiding psychological confounds and 4) maintaining the desired frequency for each stimulus type. To define the design criteria for the first two statistical objectives, we need to specify the underlying models.

Models
Following previous approaches (e.g. Liu 2004a; Liu and Frank 2004;Wager and Nichols 2003), two popular linear models are considered for the two statistical objectives; see also Friston et al. (1995), Worsley and Friston (1995), and Dale (1999). These models are: (1) where Y is the voxel-wise BOLD time series, h = (h 1 , ..., h Q ) is the parameter vector for the HRFs of the Q stimulus types, X = [X 1 · · · X Q ] is the design matrix, θ = (θ 1 , ..., θ Q ) represents the response amplitudes, Z = Xh 0 is the convolution of stimuli with an assumed basis, h 0 , of the HRF, Sγ is a nuisance term describing the trend or drift of Y , and e and η are noise. We assume a known whitening matrix, V , such that V e and V η are white noise. The whitening matrix can be obtained empirically from previous experiments; see also Wager and Nichols (2003). Model (1) is typically used for estimating the HRF and model (2) is for detecting activation. Note that, for detection problems, a basis h 0 for the HRF needs to be assumed.
To enable the use of a finite set of interpretable parameters to capture the fluctuation of the continuous HRF over time, the discretization interval (Dale 1999) is utilized for parameterizing the HRF in model (1). The length of the discretization interval, denoted by ∆T , is set to the greatest value dividing both the ISI and TR; the TR is the time interval between consecutive MR scans. The HRF parameters, captured in the vector h, then represent the heights of the HRF for each stimulus after every ∆T seconds following the stimulus onset. This parametrization is explained in detail in Kao et al. (2009b).

Design criterion
For the two statistical objectives, two popular optimality design criteria, namely A-and Doptimality criteria (Atkinson et al. 2007), are considered. A-optimality aims at minimizing the average variance of estimators of parametric functions. A D-optimal design minimizes the generalized variance of estimators of linearly independent parametric functions, or, under normality, it minimizes the volume of simultaneous elliptical confidence regions for these parametric functions at any specified confidence level. The value of the design criterion for estimating the HRF, referred to as "estimation efficiency", is denoted by F e . Likewise, the term "detection power" and the notation F d are used to indicate the value of the design criterion for detecting activation. These criteria have the forms of A is the orthogonal projection on the vector space spanned by the column vectors of A, A − is a generalized inverse matrix of A, C is a matrix of estimable linear combinations of the parameters, and r c is the number of rows of C. We note that F e and F d are "larger-the-better" criteria.
The third objective is to avoid psychological confounds. We would like a sequence that makes it difficult for a subject to anticipate future stimuli based on past stimuli. Designs minimizing the following criterion help to achieve this objective.
Here, the sub-design excluding all 0s but retaining all the stimuli of the original design is considered. In this design criterion, n is the length of the sub-design, and n (r) ij is the number of times that i and j are r elements away in the sub-design; i.e., they are, respectively, the tth and the (t + r)th elements, t = 1, ..., (n − r). P i is the specified proportion for the type-i stimulus in the sub-design which may be taken as 1/Q if there is no preference, and |a| is the integer part of the absolute value of a. R is a given integer; the unpredictability of a design increases when R increases. Therefore, F c aims at having a design with each pair appearing a number of times that is proportional to the product of the specified proportions for the stimuli. Designs minimizing this criterion are said to be Rth order counterbalanced (cf. Wager and Nichols 2003).
The fourth design criterion helps to maintain the desired stimulus frequency and is also defined on the sub-design.
where n i is the number of the type-i stimulus in the sub-design. Designs achieving the desired stimulus frequency minimize F f .
With these four individual design criteria, the family of MO-criteria is then defined as where w i s are weights selected based on the researcher's emphasis in a given study, Since these minimal values are zero, they are omitted here. Each member of this family can serve as an objective function of the search algorithm.

Search algorithm
The search algorithm of Kao et al. (2009a) is built upon the genetic algorithm technique (Holland 1975(Holland , 1992. This technique is popular for solving optimization problems, in which good solutions (parents) are used to generate better ones (offsprings). To efficiently apply this technique to our problem, well-known results about good fMRI designs are incorporated so that the search over the huge design space can be carried out more efficiently. The outline of the algorithm is as follows: Step 1. (Initial designs) Generate G initial designs consisting of random designs, an msequence-based design, a block design and their combinations. Use the objective function to evaluate the fitness of each initial design.
Step 2. (Crossover) With probability proportional to fitness, draw with replacement G/2 pairs of designs to crossover -select a random cut-point and exchange the corresponding fractions of the paired designs.
Step 3. (Mutation) Randomly select q% of the events from the G offspring designs. Replace these events by randomly generated ones. Here, an event is a stimulus or the control.
Step 4. (Immigration) Add to the population another I designs drawn from random designs, block designs and their combinations.
Step 5. (Fitness) Obtain the fitness scores of the offsprings and immigrants.
Step 6. (Natural selection) Keep the best G designs according to their fitness scores to form the parents of the next generation. Discard the others.
Step 7. (Stop) Repeat steps 2 through 6 until a stopping rule is met; e.g. after M g generations. Keep track of the best design over generations. In Step 1, m-sequences or m-sequence-based designs are generated following Liu (2004a); see also Buračas and Boynton (2002). These designs are well-known for their high estimation efficiencies. Good designs for detection are block designs. A block design is a sequence where stimuli of the same type are clustered into blocks. For example, a two-stimulus-type block design with a block size of four can consist of repetitions of {111122220000}. Repetitions of {1111000022220000} and {11112222} are other possible patterns. In steps 1 and 4, block designs with various block sizes and patterns are considered. A fraction of an m-sequencebased design or a random design is combined with a fraction of a block design to form a mixed design. These mixed designs along with random designs are also included as part of the initial designs and immigrants.

Code description
In this section, we describe our MATLAB program, including the input parameters, the way to run the program and the output variables.

Nuisance term
The nuisance term S in (1) and (2) is specified in Inp.Smat. By default, S corresponds to a second-order Legendre polynomial drift. The degree of the polynomial can be changed in 'Par_Assign.m' through the PolyOrder variable. While polynomial drift is popular (e.g. Worsley et al. 2002;Liu 2004a), other nuisance term can also be considered.

Whitening matrix
The square of the whitening matrix, V 2 , described after (2) is specified in Inp.V2. By default, the following matrix, which corresponds to a stationary AR(1) process, is considered:  (1) and (2) 2 nd -order polynomial V2 real matrix the square of the whitening matrix see (

Linear combinations of parameters
Linear combinations of the parameters of interest are specified in Inp.CX for model (1) and in Inp.CZ for model (2). These fields are, by default, set to identity matrices, allowing the study of individual stimulus effects. The number of columns for Inp.CX equals to the length of h. For a K-second HRF (by default, K = 32), the length of h is Q(1 + K/∆T ). The number of columns for Inp.CZ is Q, corresponding to the length of θ.
In addition to setting Inp.CX and Inp.CZ to identity matrices, researchers might also be interested in pairwise contrasts between stimulus types. Kao et al. (2009b) provides a systematic study of designs for convex combinations of these two interests. The example code, 'Exp_combinedInterest.m', is provided here for the case where equal weights are assigned to both interests.

Maximal values of the individual criteria
Values of max(F e ) and max(F d ) are assigned through Inp.MaxFe and Inp.MaxFd, respectively. These values are used to standardize F e and F d before combining them into an MO-criterion. Theoretical values of max(F e ) and max(F d ) are generally not available. We therefore obtain numerical approximations by performing "pre-runs" of our program. For approximating max(F e ), we use Inp.MaxFe = 1 and Inp.MOweight = [0 0 1 0]; i.e., w e = 1 in (3). This is equivalent to using the non-standardized F e as the objective function. The F e -value achieved by the design that we obtain approximates max(F e ). Similarly, we can use Inp.MaxFd = 1 and Inp.MOweight = [0 1 0 0] to find the optimal design for detection and to numerically approximate max(F d ). These approximates can then be specified in Inp.MaxFe and Inp.MaxFd for further searches for multi-objective optimal designs. 'PreRun_Fe.m' is an example code for approximating max(F e ) and 'PreRun_Fd.m' is for max(F d ). The maximal values of the other two criteria, F c and F f , are automatically calculated in our program.

Basis for the HRF
We need to assume a basis for the HRF when using model (2). A popular choice is the canonical HRF of SPM2 (The Wellcome Trust Centre for Neuroimaging 2003), a popular software for fMRI. This basis is a combination of two Gamma distributions. In our program, we use this canonical HRF, scaled to have a maximum of one, as the default setting for h 0 . The parameters, such as the time-to-peak and time-to-onset, used to create the canonical HRF can also be altered in 'Par_Assign.m'. By changing Inp.basisHRF, other basis functions can be considered as well.

Stopping rules
We consider two stopping rules. The first stopping rule terminates the search after M g generations. The second stopping rule is inspired by Liefvendahl and Stocki (2006). This second method considers the accumulated improvement of the design efficiency from the n + 1st generation to the ( + 1)nth generation; = 0, 1, 2, .... Denote the accumulated improvement by [F * (ξ * ( +1)n ) − F * (ξ * n+1 )], we stop the search at the ( + 1)nth generation when the following condition is met (for given n and δ): To use the first stopping rule, we set Inp.StopRule = 1, and Inp.numITR = M g . By setting Inp.StopRule = 2, the second stopping rule is considered. In this case, we set Inp.numITR to n (say, 100) and Inp.improve to δ (e.g., 10 −7 ).

Running the code
The m-file 'Par_Assign.m' can directly be used to perform the search for optimal designs. With user-specified parameter values, this m-file calls the subroutine fMRIMOD(Inp) to start the search. Our program requires subroutines from SPM2 (The Wellcome Trust Centre for Neuroimaging 2003) and mttfmri (Liu 2004b, see also Liu 2004a). From SPM2, we need 'spm_Gpdf.m' and 'spm_hrf.m' to calculate the canonical HRF. From mttfmri, we apply 'gen_mseq.m ', 'qadd.m', 'qmult.m', 'mseq2.m' and 'return_mtaps.m' to generate m-sequencebased designs, and 'stimpatch.m' for plotting the final design. These m-files are freely downloadable from The Wellcome Trust Centre for Neuroimaging (2003) and Liu (2004b), respectively.

Output variables
The output variables of our program are listed in Table 2. The best design achieved by our program is contained in Out.bestList and its design efficiency is Out.bestOVF. Our program also tracks the best designs over generations (Out.Lists) and their design efficiencies (Out.OVF). The value for each individual criterion F * i is also provided. Time spent on the search is recorded, too.

An example
An illustrative example is described in this section. We consider ISI = TR = ∆T = 2s, so that Inp.ISI = 2.0, Inp.TR = 2.0 and Inp.dT = 2.0. The number of stimulus types is set to Q = 2 (Inp.nSTYPE = 2). A total of 242 events (stimuli plus the control) are included in the design sequence; i.e., Inp.nEvents = 242. A second-order polynomial drift and the The canonical HRF, scaled to have a maximum of one, is used as the basis function h 0 of model (2); see also Wager and Nichols (2003); Wager et al. (2005). After discretization using ∆T , this basis is assigned to Inp.basisHRF. Following the default setting of the canonical HRF, the duration of the HRF in model (1) is K = 32s, so that Inp.durHRF = 32.0. The number of parameters contained in each h i of model (1) is therefore 17(= 1 + K/∆T ), and the length of h is 34. To investigate individual stimulus effects, we set Inp.CX = eye(34) and Inp.CZ = eye(2); they are identity matrices.
We implement our simulations by using MATLAB (version 7.3) on a Linux cluster with 64-bit AMD Opteron, dual-processor, mix of single-core node and dual-core node; each core has 2GB RAM and the Linux operating system is 2.6.9-78.0.5.ELsmp.
We first find the (near-)optimal design for estimating the HRF by setting Inp.MOweight = [0 0 1 0] and Inp.MaxFe = 1. The resulting design is presented in Figure 1 along with the curve of the achieved efficiencies over generations. In Figure 1 control. This design looks rather random and its F e -value is 39.2715. The CPU time spent for this search is 0.23 hours.
We then search for the best design for detecting activation. We set (Inp.MOweight = [0 1 0 0]) and Inp.MaxFd = 1. As presented in Figure 2, the resulting design looks like a block design. This design starts with five 0s, followed by eight stimuli of the first type and nine stimuli of the second type. The F d -value achieved by this design is 132.0670. We spend 0.13 hours of CPU time on this search.
We can also assign equal weights to the four objectives to search for a multi-objective optimal design; i.e., Inp.MOweight = [1/4 1/4 1/4 1/4]. The maximal values of F e and F d are approximated numerically; i.e., Inp.MaxFe = 39.2715 and Inp.MaxFd = 132.0670. For F c , we require a third-order counterbalancing property (R = 3), so that Inp.cbalR = 3. For F f , equal frequencies for the two stimulus types are required; i.e., P i = 1/2, i = 1, 2, and Inp.stimFREQ = [1/2 1/2]. Note that, when assigning Inp.stimFREQ, we do not take into account the number of the control event. Therefore, the actual frequency of each stimulus type in this example is less than 1/2. The number of the control event is decided by the GA search based on other requirements (design criteria). In our experience, the number of the control event is greatly influenced by the linear combinations of parameters of interest; see also, Kao et al. (2009b). The actual stimulus frequency of our designs agrees with the approximated optimal stimulus frequency of Liu and Frank (2004).
The parameter values for this last search are the same as those listed in Table 1 and those in 'Par_Assign.m'. This search requires 0.42 hours of CPU time.

Conclusion and discussion
Optimal designs are crucial to the success of ER-fMRI experiments. Due to the nature of ER-fMRI, planning a good design is very complicated. Therefore, an efficient program that helps to find such good designs is called for. In this paper, we develop a program using MATLAB to search for multi-objective optimal experimental designs for ER-fMRI. The algorithm utilized in our program is proposed by Kao et al. (2009a), which is shown to outperform the previous approaches. We describe the use of our program. An example is provided for illustration. In addition to default settings, we allow the users to assign the parameter values so that our program can achieve designs best suited to the researcher's needs.
The approach considered in this article is built upon popular linear models. However, the assumption of linearity may be invalid when the stimuli are very close, and a 'saturation' in the accumulated BOLD response is observed (Wager et al. 2005). To take into account such a nonlinear effect, Wager and Nichols (2003) propose to use two to replace the elements of the matrix Z in model (2) that are greater than two. This is also allowed in our program by setting Inp.Nonlinear to 1. Developing a more sophisticated method for incorporating such nonlinear effects can be useful.