By Tor Wager, 7/2/02
Wager, T. D. & Nichols, T. E. (2003) Optimization of Experimental Design in fMRI: A
General Framework Using a Genetic Algorithm.
, 18, 293
The idea behind the genetic algorithm is that it is
an efficient way to search through a
large and complicated space. In this case, the complicated space is the space of possible
fMRI designs. In order to use a GA, you have to 1) be able to parameterize a design as a
list of numbers, and 2) evaluate the
fitness of a design using some kind of metric testable
in computer simulation. This set of Matlab routines will help you to run a GA to
optimize experimental designs (hopefully) a little easier than if you didn’t have them. In
the best case, it’s “plug a
just create (or modify an example script) a script with
the experiment parameters you want and run it.
Text script vs. GUI
This set of functions is called from a single script, which is unique to each experiment. I
chose to implement it this
way instead of creating a graphical user interface because the
text version has several advantages. First, in the text version, you can see and save all the
input parameters, and you can easily save each script you run for each experiment as a
cond, text is easier to use (since there are a number of input parameters that
must be specified), and easier to transfer across platforms.
The principal resources you can use to help you are
This readme file
The paper (currently under review
) describing the GA
When you should use the GA (and when you shouldn’t)
Using the GA is good if you have to optimize a randomized event
related design, and you
want to detect one or more contrasts across different event types. (A contrast is a
combination of signal estimates for each event type; for example, 2 different “visual”
events versus two different “auditory” events.) Or, it could be good if you want to
recover hemodynamic response estimates (HRF shapes) for each event type. Or,
be good if you want to counterbalance the order of event types up to any number of time
steps back (i.e., each trial follows each other one equally often).
If you really want to maximize detection power for a single contrast, consider using a
blocked design. That will give you the best detection power. If you can’t use a blocked
design because you want to isolate particular events (and avoid confounds), the task can’t
be done in a blocked fashion, or for any other reason, the GA may be able t
o help. Some
block designs are better than others, especially with more than two conditions, but the
problem can be solved analytically
the strong suit of the GA is in cases where the
solution is too hard to predict a priori. ER designs with multiple con
trasts generally fall
under this category.
If you optimize for contrast detection, what the GA will try to do with your event
design is to create mini
blocks, or clumps of similar event types. This may make the
design somewhat predictable to subj
after all, it’s not truly random once you start
selecting specific trial orders. You can counteract this by specifying maximum numbers
of the same trial types that can occur in a row (hard constraint) or optimize for HRF
shape estimation and/or count
erbalancing in addition to contrast detection.
The original GA was designed to work with single events as models, not sustained
epochs. If you have an epoch
related design, there is an alternative GA function that has
worked for me in the past, but hasn’
t been extensively tested. The normal event
GA function (that does most of the work) is called optimizeGA.m. The epoch
function is called optimizeGA_epochs.m. Using this second function might be
frustrating at this point.
What to do t
o get started
Here are some simple steps to get you going.
Download the zipped tar file from the website. Unzip it using:
tar zxvf OptimizeDesign10.tar.gz
The resulting folders should be included in your matlab path. Some ways to do
this are with th
e addpath command (try help addpath), pathtool, or specifying the
GA folders in the matlab path in your .cshrc file.
Choose one of the example scripts stored in the example_scripts subfolder
(ga_example_script.m) and open it in a text editor. The script
with each line you should edit for your experiment that explains some of the
options. Type the script name to run it. Here are some of the example scripts:
A basic 2
condition design with probe periods,
pass filtering, and
condition difference detection
factor design optimizing main effects and the
mmented example script for an epoch
design with 3 epochs per trial (long trials). Use of 3
epochs/trial is mandatory with this version.
An uncommented example script for a GA that uses
ISI as a variable parameter, and attempts to
mize ISI as well as trial order.
If your script terminates successfully, a variable called M (for “Model”) will be
created in the workspace. M is a structure that contains information about your
experiment. Refer to the fields of M with the . operator;
e.g., typing M.stimlist
will show the optimized stimulus presentation sequence. You can use this to
program your experiment. Try typing:
M.ga is a nested structure contains all the input parameters to the GA, and some
s of performance across generations. More explanation on the fields can be
obtained by typing
modelDiagnostics2.m is the function that creates the M structure given a stimulus
list and a set of filtering and other experiment p
optimality vs. d
optimality is currently the default option in the GA. It minimizes the standard error of
each contrast or predictor you’re trying to estimate. This is good for testing an a prior
power to look for an A
B effect. However, it doesn’t
take into account the correlations between estimated coefficients, and so isn’t as good if
you want to know WHICH of several trial types are producing a change in a particular
brain area. For these t
ypes of questions, d
optimality may be more appropriate. D
optimality maximizes the determinant of the Fisher Information Matrix X’X, and imposes
a more severe penalty for multicolinearity than a
optimality. I recommend using this
option for estimating H
To use it, make sure GA.dflag = 1 in your script.
Optimizing the spacing of events within multi
The GA optimizes the order of discrete events or trial types in a sequence. However,
many experiments have trials with multiple parts (e.g
., encoding, delay, probe, and
feedback periods), and we’d like to separate responses to each part. The way to do that is
to vary the intervals between each part of the trial. A way to optimize the design by
selecting the exact placement of events relati
ve to one another is provided in ga2.m, in
the ga2 subfolder. Rather than optimizing the order, GA2 optimizes the timing
parameters, including presentation length and jitter. It may be desirable to optimize both
the order and timing in sequences with dep
endencies, but that’s not quite implemented
yet. GA2 has its own help files, located in the GA2 subfolder.
“Jitter” in designs
There are two easy ways to introduce random jitter into GA designs. The first is by
setting the field GA.freqConditions to b
e a vector that sums to less than 1 (this is set in
the GA script that you run, e.g., ga_example_script.m. freqConditions specifies the
proportion of the total number of events for each trial type. For example, if you have two
conditions, then [.5 .5] me
ans 50% event type 1 and 50% event type 2. [.4 .6] signifies
40% event type 1 and 60% type 2. [.4 .4] signifies 40% for type 1 and 40% for type 2;
the remaining 20% will be blank rest intervals.
The second way is to create an extra condition that will
be the “baseline” condition. If
you have two event types plus rest, you might set GA.conditions = [1 2 3] to specify three
trial types: type 1, type 2, and rest. GA.freqConditions = [.4 .4 .2] signifies 40% for type
1, 40% for type 2, and 20% rest (condi
tion 3), and is equivalent to the example above.
How designs are coded
Designs are coded as a column vector of integers. Each integer represents a different
condition, or “event type,” and the elements represent the position in time. Zeros
sts. So the vector represents the stimulus presentation state at each time
interval throughout the run, where the intervals are defined at an arbitrary resolution
(stored in GA.ISI). ISI stands for inter
stimulus interval; in practice, that’s a pretty go
time resolution to use, as design vectors too fine
grained take a long time to process.
Another way of saying this is that you must specify some stimulus
presentation state at
some regular interval (ISI). Zeros indicate rest intervals, and they are no
t modeled in the
model matrix or considered in counterbalancing or other fitness measures.
Much of the GA script and the main GA function (optimizeGA.m) are for the purpose of
constructing a number of randomized lists of event codes according to your
ecifications. These vectors, concatenated into a matrix of column vectors stored in the
variable listMatrix, forms the starting point for the GA. The rest of optimizeGA.m
creates model matrices (the X matrix in a statistical GLM analysis) from each vecto
tests them with fitness (or “objective”) functions, and performs crossover on the best half
of each generation. This is illustrated graphically above.
GA timeline: Sequence of operation
Here’s a slightly more detailed look at what the GA does at w
hat point during its
creates GA structure variable with all necessary input values
For 1 to nmodels (all these names refer to variables defined in the
script), runs the function optimizeGA.m
the main GA function;
can also substitute optimizeGA_epochs.m
(discretion is advised). This function does the following:
SETUP: set up contrasts, nonlinear saturation, periodic rest
intervals, filtering and intrinsic autocorrelation matrices, task
switching (trial order effec
ts) and mixed block/ER setup (if
specified), create initial population of randomized designs (the
constraint criteria setup, approximately in that
order. Prints output so you can check it.
#1 #2 #3
1 2 1
2 3 4
4 2 3
2 1 4
1 4 3
1 2 4
3 1 2
#1 #2 #3
1 2 1
4 3 2
1 4 3
1 2 4
3 1 2
iterate over generations
Loop through generations; for each generation, and e
ach list in
each generation runs these functions in order (if necessary
given your input specifications):
dna2rna.m: Create intermediate design vector (“rna”)
from the original event list (“dna”)
getCounterBal.m: Counterbalancing fitness calculation
reqDev.m: Calculate deviation from input
frequencies for each event
type (fitness measure)
getMaxInARow.m: Calculate maximum repeats for the
designvector2model.m: Contrast detection model
pute contrast detection
tor_make_deconv_mtx2.m: Construct model matrix for
HRF shape estimation efficiency.
CalcEfficiency.m: Run again for HRF shape.
After each generation
GetOverallFitness.m: calculates overall fitness score
given results o
f all fitness tests.
Save the best of this generation.
Breedorgs.m: Take the best half (or roughly so,
depending on value of “alpha” input parameter) of
designs (“dna”) and crossover, exchanging upper
halves. Returns a new population of design vectors.
rint output for the generation.
Every 10 generations, save intermediate results in
listMat_working. If the GA quits partway through, you
can use the listMatrix (the population list) variable in
this file to start again where the GA left off.
After all g
Create M structure using the best overall design, with
all diagnostic measures.
modelDiagnostics2.m: adds resampled and filtered
design matrices to M structure, as well as efficiency,
variance inflation, condition number, predictor
atrix, average time between repeats, and
other design fitness/evaluation measures. Also plots
predictors and contrasts graphically, if specified.
M structure is saved in model#_<date>.mat, where # is
the #’th iteration of the GA (out of nmodels).
nal functions not further explained here
A number of parameters coded in the scripts and functions are “add
ons” designed to
perform specific functions. If you generally leave these at their default values when
getting started, you should be OK. See the
simple ga_example_script.m for defaults.
Some of these functions (and their defaults) are:
GA.doepochs = 0;
Build epochs instead of events from design vector
Use with optimizeGA_epochs.m
GA.numhox = 0;
“Hox” genes (by analogy) allow you to op
timize things like ISI
directly with the GA. This is experimental. It’s worked for me,
but it hasn’t been widely tested, and probably only works under
specific conditions. The idea is simple; include ISI, or any other
input parameter, in the GA with a r
andom value between two
limits. Then optimize that parameter value with the GA. In
practice, some input parameters interact with how model matrices
are constructed, and some are very computationally intensive to
implement. The current ones are hard
d to represent [ISI TR
cuelen cuerest stimlen stimrest resplen resprest (in s)] in that order.
ISI should work, TR does not work (specify and leave at whatever
your TR is), and the rest code for epoch lengths for a long single
trial epoch design with mult
iple periods (cue, stimulation,
response) within the trial.
GA.hoxrange = ;
Starting ranges for each Hox parameter to optimize.
GA.alph = 2.1;
"selection pressure": higher is more extreme selection.
GA.lowerLimit = ;
From an older version, not used;
don’t worry about it.
GA.restevery = ;
GA.restlength = ;
These two parameters insert rest intervals (0’s) of length restlength
(in units of the ISI) every restevery ISIs. For example, restevery =
5, restlength = 3 would insert 3 zeros into each des
ign vector (3
ISIs) every 5
element. This is a convenient way to model
instruction and probe periods, which occur periodically and
regularly, but are not of interest in the study. These periods do
affect the design efficiency, so including them makes t
simulation a bit more realistic.
GA.trans2switch = 0;
This parameter is set to 1 when you have a history
design. For example, if your predictors of interest are switches
between two event types [1 2], and you set trans2switch to 1,
ns 1 and 2 will represent events A and B, where A is
preceded by A and B is preceded by B (“no switch”). Conditions 3
and 4 will be created automatically, and they will represent A and
B, where A is preceded by B and B by A (“switch”) condition.
s will then have 4 conditions, potentially representing
main effects of A vs B, main effects of switching, and interactions.
This special feature works only when GA.conditions = [1 2].
GA.trans2block = 0;
This parameter is for mixed blocked
d designs. It takes
related conditions you specify, and doubles the number
of conditions. If you specify [1 2 3 4] as the conditions, it will
create [5 6 7 8] as well, and alternate between blocks of 1:4 and
5:8 in an ABAB fashion. Contrast le
ngths should double, so that if
you were interested in a [1 0 0
1] contrast, you would specify that
as [1 0 0
1 1 0 0
1]. The contrast for block A vs. block B would
be [1 1 1 1
1]. The blocks alternate between A and B
every # elements, where
# is specified by GA.restevery. If you
want to use a mixed blocked/ER design but don’t want to include
rest intervals, set restevery to the number of ISIs in a block, and
restlength to 0. If you use both trans2switch and trans2block, you
would start wit
h [1 2] in the GA.condition vector, and specify eight
values in the contrasts. [1 1
1 1 1
1] would be a main
effect of switching, and [1
1] would be a main
effect of event 1 vs. event 2.
GA.dofirst = 0;
This was included because i
n designing mixed block/ER switching
experiments, it may be desirable to code the first trial of every
block as a special trial type, neither a switch or a non
this when using trans2switch and trans2block, but do not lengthen
the contrast vect
ors to account for this; this special “first trial”
predictor is automatically added to the end of the contrast with a
weight of 0.
What usually goes wrong
Most of the time, an error will be caused by one of several common problems relating to
ification of input parameters.
The contrasts are the wrong length. There should be one column in the contrast
matrix for each condition (not including the intercept). If you’re using the
trans2switch option (history
dependent designs) or the trans2blo
ck option (mixed
block/ER design), the number of columns in the contrasts should double
your conditions are [1 2], 8 if they’re [1 2 3 4]. If you’re using both options,
they’ll double twice. The intercept shouldn’t be included in the contrast weight
If you’re using the dofirst option, the contrast should not reflect this; an extra 0 is
added onto the end.
The contrastweights vector is not the same length as the number of contrasts
(rows in the contrast matrix, GA.contrasts) you specified.
using trans2block, trans2switch, or dofirst when you didn’t mean to. A 1
indicates “yes, use this” and a 0 is “off.”
Restarting the GA with a partially completed simulation
If, for some reason, the GA completes a number of generations, and then the si
stopped partway through, you can re
start it approximately where it left off. Every 10
generations, the workspace of the main GA function is saved to listMat_working.mat. To
Load that mat file
Run the first part of your script and b
reak partway through, or paste the GA.xxxx
definitions from your script into the workspace.
Type GA.listMatrix = listMatrix;
Run the bottom part of the script (paste into the workspace) after “End User
Input.” The function optimizeGA will see that ther
e is an existing listMatrix to
start with (the one from your old run), so it will use that instead of creating a new
The “truth” of the simulations depends on the validity of the assumption that the fMRI
response is a linea
r function of the input. One of the things that means is that the
response to each event should be the same as to all other ones of the same type,
regardless of the history of past events or the identities of the proximal events. One of
iolations of this assumption is that the fMRI signal saturates when
stimuli of the same type are presented close together in time (less than about 2 s): After
the first stimulus, subsequent repeated events within the next few seconds (at least)
s additional activation.
One simple way to account for this effect in the design is to model the saturation effect as
an activation “ceiling.” The response is allowed to increase linearly up to some multiple
of the unit response; after that, it is “cap
ped” at that value and not allowed to increase.
This type of model is incorporated into the GA under the “nonlinthreshold” parameter.
The value entered by the user is the multiple of the unit response at which the predicted
signal is capped. This nonlin
earity model effectively implements a penalty for designs
that place stimuli too close together. Predictors that are “capped” suffer from restricted
range and possibly higher collinearity with other predictors. Without this model,
efficiency would increa
se without bound as events are spaced closer together, particularly
for difference contrasts.
This model is not as accurate as more complex models of nonlinearity in BOLD signal.
Its main advantage in this context is that it’s very computationally effici
ent to implement,
so it’s feasible to test the hundreds of thousands of designs necessary for GA operation
with this model.
Other functions included in the package (e.g., parameter
There are one or two other functions that are not documented.
…a function for creating surface maps of exhaustive search over parameters,
particularly ISI and jitter (created by inserting 0’s into the design vector).
…generates a population of random designs and tests your de
sign (in the M
structure), both in contrast efficiency and estimated t
scores, against this random
population. Graphical output.
…and a couple of other ones in there, too.