The weird sequence generated by @stim_analyze

Hi afni experts,
I have three conditions in my new event-related fMRI experiment: exp1,exp2 & control. I am interested in two t-test pairs among these conditions: exp1-control & exp2-control. So I search the best experiment sequence use the method provided by @stim_analyze: generating sequence using make_random_timing.py, calculating NSD using 3dDeconvolve and choosing the sequence using the criteria with the sum of the NSD of the two t-test pairs.
But the best sequence is very weird. Some events in this sequence are too close in time and some other events’ ISI are too long.


The ISI statistics :

[size=small]data              min      mean     max     stdev
-----------     -------  -------  -------  -------
pre-rest          4.000    9.000   14.000    4.163
post-rest        14.000   28.500   44.000   12.369

run #0 ISI        0.000   12.000   56.000   13.156
run #1 ISI        0.000   13.130   68.000   17.277
run #2 ISI        0.000   12.870   62.000   15.630
run #3 ISI        0.000   12.870   68.000   16.807

all runs ISI      0.000   12.717   68.000   15.542[/size]

Visualizing the sequence:
https://ww2.mathworks.cn/matlabcentral/answers/uploaded_files/116192/55.png
The code used was copied from @stim_analyze without big changes:


set num_stim    = 3
set num_runs    = 4
set pre_rest    = 0
set post_rest   = 6
set min_rest    = 0 
set tr          = 2.0 
set stim_durs   = 2
set stim_reps   = 8 # This specifies the number of repetitions of each stimulus type, per run
@ run_lengths   = $post_rest + $pre_rest + 186 * `printf '%.f' "$tr"` # set the total time, per run (in seconds)
set iterations  = 1000000   # number of iterations to compare
set seed        = 123       # initial random seed
set outdir      = optimal_results
set LCfile      = NSD     # file to store norm. std. dev. sums in
set datapath    = optimal_results
set labels      = "exp1 exp2 control"
# recreate $outdir each time
if ( -d $outdir ) then
   echo "** removing output directory, $outdir ..."
   \rm -fr $outdir
endif
echo "++ creating output directory, $outdir ..."
mkdir $outdir
cd $outdir
echo -n "" > $LCfile
echo -n "iteration (of $iterations): 000000000"
# HRF
#3dDeconvolve -num_stimts 1 -polort -1 -nodata 51 0.5 \
#   -stim_times 1 '1D: 0' 'BLOCK5(1,1)'       \
#   -x1D stdout: | 1dplot -stdin -THICK -del 0.5 -jpeg bk5
# run the test many times
set numlen=`expr length $iterations`
foreach iter (`count -digits $numlen 1 $iterations`)
@ seed = $seed + 1
make_random_timing.py -num_stim $num_stim -stim_dur $stim_durs  \
                -num_runs $num_runs -run_time $run_lengths              \
                -num_reps $stim_reps -prefix stimes.$iter               \
                -pre_stim_rest $pre_rest -post_stim_rest $post_rest     \
                -min_rest $min_rest                                     \
                -stim_labels $labels                                    \
                -seed $seed                                             \
                -tr $tr                                                 \
                -tr_locked                                              \
                -max_consec 2                                          \
                -show_timing_stats                                      \
                -save_3dd_cmd cmd.3dd.$iter                             \
                        >& out.mrt.$iter
3dDeconvolve                                                     \
    -nodata 756 2.000                                            \
    -polort 3                                                    \
    -concat '1D: 0 189 378 567'                                  \
    -num_stimts 3                                                \
    -stim_times 1 stimes.${iter}_01_exp1.1D 'BLOCK5(1,1)'    \
    -stim_label 1 Global_constant                                \
    -stim_times 2 stimes.${iter}_02_exp2.1D 'BLOCK5(1,1)'         \
    -stim_label 2 Scr_nature                                     \
    -stim_times 3 stimes.${iter}_03_control.1D 'BLOCK5(1,1)'       \
    -stim_label 3 Scr_constant                                   \
    -jobs 10                                                     \
    -gltsym 'SYM: exp1 -control'                      \
    -glt_label 1 locals                                           \
    -gltsym 'SYM: exp2 -control'                 \
    -glt_label 2 global                                         \
    -x1D X.xmat.1D  >& out.3dD.$iter
set nums = ( `awk -F= '/LC/ {print $2}' out.3dD.$iter` ) 
# make a quick ccalc command
        set sstr = $nums[1]
        foreach num ( $nums[2-] )
            set sstr = "$sstr + $num"
        end
set num_sum = `ccalc -expr "$sstr"`
echo -n "$num_sum = $sstr : " >> $LCfile
echo    "iteration $iter, seed $seed" >> $LCfile
foreach i (`count -digit 1 1 $numlen`)
echo -n '\\b' >>tem.1D 
end
set xxx=`cat tem.1D`
echo -n "${xxx}$iter"
rm -f *.1D cmd.3dd* out.* 
end
echo ""
echo "done, results are in '$outdir', LC sums are in '$LCfile'"
sort -n $LCfile | head -20 > opt_results

I want to know what’s wrong with my code.
Thank you in advance!

Hi Zhang,

Looking over the timing information, for a single run, say,
I see:

  • ignoring pre/post rest, run lengths of 372 (186*2) seconds
  • 3 stim classes
  • 8 stim repetitions per run, with stim durations of 2 seconds
  • min_rest = 0 s

So out of a 372 s run, there is 48 s (3 * 8 * 2) of stimulation,
leaving 324 s of rest. Since there are 24 events (3*8), that
means an average of approximately 324/24 ~= 13 s of rest
between events (or the start/end of the run).

That is a lot of rest time. Do you intend to have so much?

  • rick

Hi Rick,
Thanks for replying me!
I thought the longer ISIs the better the deconvolved results. And I thought we needed a slow event design to get better signal change curve of each stimulus. So the long rests were expected. But I am not sure about this. I thought the distribution of events should look evener but some ISIs are too small and the others are too big.
I also tried to reduce run_lengths to 180 includes pre_rest and post_rest:


[size=small]set num_stim    = 3
set num_runs    = 4
set pre_rest    = 0 #pre_rest will be added manually later in the formal experiment
set post_rest   = 10 
set min_rest    = 10
set tr          = 2.0 
set stim_durs   = 2
set stim_reps   = 8 # This specifies the number of repetitions of each stimulus type, per run
@ run_lengths   = 180 * `printf '%.f' "$tr"`[/size]

The optimal sequence looks evener. But the post_rest is still too long to take for me:


[size=small]timing_tool.py -multi_timing stimes.${iter}_0*                  \
                 -run_len $run_lengths -multi_stim_dur $stim_durs \
                -multi_show_isi_stats

ISI statistics (3 elements) :
                        total      per run
                       ------      ------------------------------
    total time         1440.0       360.0    360.0    360.0    360.0  
    total time: stim    192.0        48.0     48.0     48.0     48.0  
    total time: rest   1248.0       312.0    312.0    312.0    312.0  
    rest: total isi    1142.0       290.0    280.0    280.0    292.0  
    rest: pre stim       18.0         0.0     12.0      6.0      0.0  
    rest: post stim      88.0        22.0     20.0     26.0     20.0  
    num stimuli          96          24       24       24       24  

                         min      mean     max     stdev
                       -------  -------  -------  -------
    rest: pre-stim       0.000    4.500   12.000    5.745
    rest: post-stim     20.000   22.000   26.000    2.828
    rest: run #0 ISI    10.000   12.609   36.000    5.606
    rest: run #1 ISI    10.000   12.174   20.000    2.480
    rest: run #2 ISI    10.000   12.174   16.000    2.329
    rest: run #3 ISI    10.000   12.696   22.000    3.548
    all runs: ISI       10.000   12.413   36.000    3.674[/size]

And the NSD of optimal sequence becomes bigger (from 0.332600 to 0.354300).
Do you have suggestions for me to improve?
Yu Zhang

Hi Yu Zhang,

If you actually want to deconvolve to get the response
shape (e.g. with TENTs), then yes, it is useful to have
the longer ISIs. But even for that, it is not necessary
to exceed about 15 seconds (depending on how long you
expect the brain to actually respond to a stimulus, which
might be different from the presentation duration).

So it is highly dependent on the experiment design. But
having ISIs of 20-36 seconds seems very long. Again, the
maximum can be limited (see -max_rest or the advanced
usage), but that would probably mean similarly limiting
the average.

Decide up front what the pre-, post- and average ISI rest
should be.

If the general parameters are changed, the NSD values
cannot be directly compared.

  • rick

Hi Rick,
Thanks a lot!
Besides the ISIs, the post stimuli rests are also too long. Are there any parameters in make_random_timing.py to make the rest at the end of each run changing at a specific range (e.g. from 10 to 16s)?
Yu