Stimulus counter-balance is important for many experimental designs. This command-line software creates pseudo-random sequences with arbitrary levels of counter-balance. "Path-guided" de Bruijn cycles may also be created. These sequences encode a hypothesized neural modulation at specified temporal frequencies, and have enhanced detection power for BOLD fMRI experiments. These sequences are particularly well-suited for carry-over fMRI experiments.

This paper describes the use of de Bruijn cycles in neuroscience experiments in general, and fMRI adaptation experiments in particular:

• GK Aguirre, MG Mattar, L Magis-Weinberg. (2011) de Bruijn cycles for neural decoding. NeuroImage 56: 1293-1300

The sections below describe the use of our open source software for the creation of these sequences.

The software contains Hamiltonian Cycle code that was developed by Basil Vandegriend and Joseph Culberson at the University of Alberta (Copyright 1998). Please see their Conditions of Use page regarding the further use and dissemination of this code. In a personal communication (July 4, 2013), Dr. Vandegriend has provided authorization to release the code under the BSD-3 license.

Calling with either ''-h'' or with no parameters produces this help text:

<code - debruijn -h>
Usage: debruijn [-t | -v] <k> <n>
	[<B> <guide function> <neural model>]
	[<neural model #2>] [<neural model #3>] [-eval <SOA>]
Generate random and "path guided" de Bruijn cycles
  
Default output has one line of labels, separated by commas
	-t: terse output (no delimiters)
	-v: verbose output in "necklace" format
  
  	k: number of letters (maximum is 36)
	n: length of the word (i.e., level of counterbalance)

Specify parameters for path-guided cycle
	B: number of bins for available transitions (integer between 1 and k^2)
	guide function: can be specified in three different ways -
		HRF: guide with power spectrum of the BOLD HRF  -OR-
		path to .txt file containing guide function -OR-
		[Tmin,Tmax]: range of the periods of sinusoids
		(in units of labels) to use as a guide function
	neural model: path to .txt file containing neural model for optimizing
	neural model(#2,#3): optional; detection power only reported
  	-eval: Evaluate detection power
	SOA: stimulus-onset asynchrony (in ms)

Code by: Marcelo Mattar ([email protected])
	Dongbo Hu ([email protected])
	Hamiltonian circuit code from Basil Vandergriend
	(http://webdocs.cs.ualberta.ca/~joe/Theses/vandegriend.html)

Full counterbalance is present with the de Bruijn cycle; for a linear sequence to exhibit complete counterbalance, the last n elements of the sequence must be presented prior to the start of the sequence to establish the context for the first stimulus. Moreover, it is optimal in BOLD fMRI studies to present 8-15 seconds of stimuli from the end of the sequence at the start, to allow the hemodynamic response to reach steady-state. The neural data corresponding to these initial stimuli may be dropped, resulting in a final fMRI signal that represent the entire sequence of stimuli, with the delayed and dispersed BOLD fMRI response to the neural sequence “wrapping around” from end to start. A similar "repeat and clip" technique may be used at the boundaries of scans if a lengthy de Bruijn sequence must be broken into portions for subject comfort.

The de Bruijn cycle will be given as a sequence of labels. The characters ''0-9'' are used for the first 10 labels; the characters ''A-Z'' are used thereafter. The maximum k permitted is 36.

Path-guiding selects amongst available nodes whose transition match that called for by the guide function. To provide stochasticity, the available nodes are cast into bins, and the next node in the sequence is selected from within the best matching bin. A large number of bins produces a path-guided sequence of greater initial fidelity, but with decreasing precision as the sequence progresses. Conversely, a small number of bins produces a consistently less precise approximation of the guide function.

Better performance is provided when the number of nodes in each bin is roughly equivalent. For this reason, we recommend selecting B such that [k² modulo B] = 0.

The neural matrix represents the expected neural response to stimuli (or transitions between stimuli). Using this information, "path guided" sequences may be created which are optimized for the detection of the hypothesized neural effect.

The parameter provides the path to a k x k matrix of positive floating point values in a plain text file. Columns of the matrix are delimited by spaces, and rows by carriage returns. Example neural matrices are included in the repository in .txt files.

NOTE The code currently requires that the neural matrix file end with the last row of the matrix. Ending the file with an additional newline or carriage return character causes the program to fail with the message: ''The neural model file must have a nxn matrix of floating points''.

For some stimuli (or transitions), the entry in the neural matrix may be left undefined. For example, from a target stimulus or for the stimulus following a "null" trial. For these transitions, a value of ''-1'' may be entered in the neural model matrix. These undefined transitions are randomly distributed throughout the sequence.

The direct effect is the neural response to a stimulus itself (as opposed to stimulus context). Classically, fMRI designs are optimized to maximize the direct effect. For example, the classic "block design" of alternating between two stimulus conditions every 30 seconds positions the expected neural modulation at an ideal temporal frequency with respect to the BOLD signal and noise properties.

To guide de Bruijn sequences for the detection of direct effects, construct the neural model to contain columns of expected relative magnitude of neural response. For example, an ever larger amplitude of response to the A, B, C, and D stimuli may be represented as:

The resulting sequence will order the stimuli (within the constraints of counter-balance) to modulate the direct effect.

The transitions between stimuli may be expected to modulate neural response. These are carry-over effects:

• GK Aguirre (2007) Continuous carry-over designs for fMRI Neuroimage, 35, 1480-1494.

Neural adaptation (or habituation) is an example of a (symmetric) carry-over effect.

To design a sequence optimized for these effects, create a neural model which represents the effect of pair-wise transitions between the stimuli. The row index indicates the prior stimulus, and the column index indicates the current stimulus. This model, for example, contains both a symmetric (adaptation) effect, and an asymmetric bias effect:

In this example, the carry over effects are predicted to be larger for the transition of ''[A ⇒ B]'' as compared to ''[B ⇒ A]'' (for example, transitions towards the "A" end of a linear stimulus space are more salient and perceived as larger). The transition values represented are:

Transitions between identical stimuli may be modeled as an undefined transition (''-1'') or as having a transition of zero (''0''), depending upon the hypotheses of the study and the modeling approach of the analysis.

Any valid de Bruijn cycle will provide the specified level of counter-balance amongst the labels, and thus stimuli, of the experiment. Not all orderings of stimuli, however, are equally useful for neuroimaging experiments. Because of the signal and noise properties of fMRI, some temporal frequencies of neural modulation are preferentially detected by the method. The path-guided approach to de Bruijn cycle generation encodes in the ordering of the stimuli a hypothesized neural modulation at these preferred frequencies.

To do so, we define a guide function in one of three (mutually exclusive) ways:

• Enter the word ''HRF''. A guide function will be internally generated that has the same power spectrum as the BOLD hemodynamic response function, although with no power below 0.01 Hz. This will generally produce a sequence with good detection power and a stochastic variation in stimulus transitions, although further improvements can be obtained by using a guide function that is positioned solely at a few or a single frequency. NOTE An SOA must be specified using the ''-eval'' option for this guide function to be defined.
• Enter a range of periods (in units of labels), as described by [Tmin,Tmax]. A guide function will be internally generated as the sum of sinusoids of random phase, each having a period equal to an integer between T_min and T_max. Enter the same value for T_min and T_max to guide the modulation at a single frequency. Generally, when ''(1e5 / SOA in ms) > T > (1e4 / SOA in ms)'' encoded modulations will be within a detectable range for the BOLD system.
• Provide the path to an external file, which contains kⁿ space-separated floating point values. Each element will correspond to a relative transition in the output sequence. The elements in the guide function may be either positive or negative; the vector will be normalized.

For fMRI, the optimal range of temporal frequencies is ~0.01-0.1 Hz (E Zarahn, GK Aguirre, M D'Esposito. (1997). Empirical analyses of BOLD fMRI statistics. I. Spatially unsmoothed data collected under null-hypothesis conditions Neuroimage, 5, 179-197). Higher frequencies are attenuated by the dispersed hemodynamic response, while lower frequencies are lost within the pink (1/f) noise of the system. Perhaps surprisingly, it is not the case that the best performing sequence is always given by the lowest temporal frequency, even ignoring the presence of 1/f noise. For some neural models, the available sets of transitions more readily fit higher temporal frequencies. This seems to be true in particular for one-dimensional stimulus spaces with a small number of stimuli. We recommend a search over a range of guide functions.

The interpretation of the path-guided sequence in Hz requires a specification of the stimulus-onset asynchrony, described next.

Given a particular hemodynamic response function, and a model of the 1/f noise of the BOLD fMRI system, detection power (Thomas Liu. (2004) Efficiency, power, and entropy in event-related fMRI with multiple trial types. Part II: design of experiments. Neuroimage 21: 401-13.) is the proportion of neural variance which appears in the imaging signal. It is calculated as the variance of the hypothesized neural modulation proportional to sequential stimulus distance as the denominator, and the numerator as the variance of that modulation after passing through the BOLD fMRI system. In our implementation, the BOLD system is modeled using a standard, population averaged hemodynamic response (GK Aguirre, E Zarahn, M D'Esposito. (1998). The variability of human BOLD hemodynamic responses NeuroImage, 8, 360-369.) and the elevated 1/f noise range by a 0.01 Hz, high-pass notch filter.

When the ''-eval'' flag is set, the stimulus-onset asynchrony (SOA) parameter is specified in units of milliseconds (i.e., the time that elapses between the start of one stimulus and the start of the next stimulus). Along with the path-guided de Bruijn sequence, the routine then returns both the detection power, and the correlation coefficient of the guide function (input) with the sequence of distances between stimuli generated (output).

file name	image	description
k=17_16stim_dioct.txt		16 stimuli in a di-octagon arrangement (Euclidean geometry), stimulus zero as a null-trial, perfect repetitions have undefined distance
k=6_5stim_linear.txt		5 stimuli in a linear array, stimulus zero as null-trial, perfect repetitions have zero distance
k=9_8stim_circular.txt		8 stimuli with a circular similarity, stimulus zero as a null-trial, perfect repetitions have zero distance

This shell script may be used to search for and retain the best debruijn sequence, or to explore the effect of changes over the parameters of the search.

#######################################
# Finds the best de Bruijn sequence
#######################################

clear

#######################################
# PARAMETERS:
#######################################
k=11
n=3
numbins=11
distfile=k=11_10stim_linear.txt
isi=1100
tmax=80
numiter=100
#######################################

maxdetecpow=0

for (( tmin = 10; tmin <= 80; tmin++))
do

for (( i = 1; i <= $numiter; i++ ))
do
	echo i=$i
        ./debruijn -t $k $n $numbins $distfile [$tmin,$tmax] -eval $isi -debug > myoutput.txt

	detecpow=`cat myoutput.txt | grep DETECTION | awk '{ print $3 }'`        #saves the detection power from the output text in a variable
	correlation=`cat myoutput.txt | grep CORRELATION | awk '{ print $3}'`    #saves the correlation from the output text in a variable

        if [ $detecpow ]                                                         #if the detection power was succesfully saved in a variable
        then
	        detecpow=`echo $detecpow | bc`                                   #converts from str to float
                compare_result=`echo "$detecpow > $maxdetecpow" | bc`            #compares with the current max value
                if test $compare_result -gt 0                                    #if the current one is greater
	        then
		        correlation=`echo $correlation | bc`                     #converts from str to float
                        echo CORRELATION = $correlation                          #outputs a text showing the correlation value for the sequence with max detection power
			echo DETECTION POWER = $detecpow                         #outputs a text showing the maximum detection power so far
			
			maxdetecpow=$detecpow                                    #saves the maximum detection power in a variable for future comparisons
			cp myoutput.txt bestoutput.txt                           #saves the output with the best detection power in bestoutput.txt
                fi
        fi

done

done