MUX EPI

Simultaneous Multi-Slice (SMS) EPI at CNI (a.k.a., MUX EPI)

NOTE: if nothing else, be sure to read and understand the Calibration Images section!

The CNI, in collaboration with GE, has implemented simultaneous multi-slice EPI (also known as multiband EPI, multiplexed EPI, or, as we like to call it, "mux EPI"). We have functioning sequences for both gradient-echo (GRE) and spin-echo (SE) EPI. GRE-EPI is useful for BOLD imaging and SE-EPI is useful for diffusion imaging and other quantitative anatomical methods.

Calibration Images

There are two ways of collecting the calibration images for the mux reconstruction. One is the internal calibration where the calibration and mux acquisition are done in a single scan, the other is the external calibration where the calibration and mux acquisition are two separate scans.

Internal Calibration

If the mux scan contains the internal calibration, the first few acquired images will be fully-sampled and used as reference images for the reconstruction. It is critical that the subject not move during this calibration phase, or the whole dataset may be impossible to reconstruct without significant artifacts (if at all). The number of reference images depends on the mux and arc factors as well as the CV that controls the number of reference scan repetitions (num_mux_cycle, exposed in the User CVs as "Number of calibration scans"). The number of TRs needed for the reference images is: arc * mux * num_mux_cycle. (num_mux_cycle default is 2). The NIFTI file that you get will include these extra images at the beginning, so you should exclude them from your analysis.

Here are some useful equations to keep in mind if you are setting up a mux sequence for an fMRI experiment that is n TRs long per run, not including calibration time (i.e., you want to have n usable TRs from a scan):

• total number of TRs acquired = n + calibration TRs
• number of calibration TRs = mux * arc * num_mux_cycle
• number of TRs in nifti = n + num_mux_cycle
• number in the "phases per location" field: n + mux * num_mux_cycle

For example, in order to have 100 usable TRs from a run with mux=3, arc=2, and num_mux_cycle=2, you would manually enter 106 in the phases per location field, the total scan time would automatically adjust to be 112 TRs, and your experiment should begin 12 TRs after you trigger the scan to begin. If you are using the internal calibration, the nifti for this run would contain 102 volumes, the first two of which are the reconstructed calibration scans that should be skipped for BOLD analysis.

External Calibration

External calibration uses a separate EPI scan of the same spatial coverage to collect the calibration data. It can be either a single-band EPI scan or a mux EPI scan of the same mux factor. The single-band external calibration is usually recommended over the mux external calibration because the single-band image has a better brain tissue contrast therefore it can perform better registration between the EPI and anatomical images.

When using external calibration, the number of reference scan repetitions in the mux scan (num_mux_cycle) can be set to 0, so that the actual data acquisition starts the same time as the scan starts. Therefore if you want to set up an fMRI experiment of n TRs, you would set the number of phases to n. For single-band external calibration, the calibration images are NOT included in the reconstructed nifti of the mux scan; whereas for mux external calibration scan, the calibration images ARE included in the reconstructed nifti.

Sample Scanning Protocols

For each SMS EPI scan, you can either use internal calibration or set up a separate calibration scan (usually non-SMS, or single band, marked as 'sbref' in series description) for the SMS reconstruction, and another short SMS scan with reversed phase encoding direction (marked as 'pe1' in series description) for geometric distortion correction. The scan coverage and matrix size are the same for the calibration and the SMS scans.

SMS fMRI

BOLD EPI

• fMRI mux8 1.6mm: SMS factor 8, 1.6 mm^3 voxel size, TR 1.1s, number of muxed slices 10 (80 unmuxed slices, 12.8 cm)
• fMRI mux8 2mm: SMS factor 8, 2 mm^3 voxel size, TR 720ms, number of muxed slices 8 (64 unmuxed slices, 12.8 cm)
• fMRI mux8 2.8mm: SMS factor 8, 2.8 mm^3 voxel size, TR 400ms, number of muxed slices 6 (48 unmuxed slices, 13.4 cm)
• fMRI mux6 2.4mm: SMS factor 6, 2.4 mm^3 voxel size, TR 710ms, number of muxed slices 10 (60 unmuxed slices, 14.4 cm)

If you want to use in-plane acceleration to reduce the EPI distortion, we also have a couple of scans set up for that purpose:

• fMRI mux4r2 2.2mm: SMS factor 4, 2x in-plane acceleration, 2.2 mm^3 voxel size, TR 950ms, number of muxed slices 16 (64 unmuxed slices, 14.1 cm)
• fMRI mux3r2 2.4mm: SMS factor 3, 2x in-plane acceleration, 2.4 mm^3 voxel size, TR 1s, number of muxed slices 18 (54 unmuxed slices, 12.9 cm)

We have one additional option, a quiet fMRI scan. This is enabled by the Acoustic Reduction Technology (ART) option in the SMS fMRI sequence. The ART scan is quieter, but the acquisition will be slower and the geometric distortion will be worse.

• fMRI mux3r2 ART medium 2mm: SMS factor 3, 2x in-plane acceleration, 2 mm^3 voxel size, TR 2s, TE 30ms, number of muxed slices 23 (69 unmuxed slices, 13.8 cm)
• fMRI mux3r2 ART high 2mm: SMS factor 3, 2x in-plane acceleration, 2 mm^3 voxel size, TR 2.25s, TE 30ms, number of muxed slices 23 (69 unmuxed slices, 13.8 cm)

Multi-echo EPI

We also have a multi-echo SMS sequence where you can measure up to three echoes per TR. This allows studying the T2* effect of the BOLD signal.

• fMRI mux3r2 me3 3mm: SMS factor 3, 2x in-plane acceleration, 3 mm^3 voxel size, TR 1s, TE1 13ms, TE2 29.2ms, TE3 49ms, number of muxed slices 13 (39 unmuxed slices, 11.7 cm)

Advanced scan parameters

• CV3 - Inplane acceleration (arc factor): determines the ky undersampling. Currently limited to an integer. E.g., 2 will acquire every other ky line.
• CV9 - Phase encoding gradient polarity: determines the acquisition direction in ky, usually labelled as pe0 and pe1 in series description.
• CV17 - Number of mux calibration scans: how many internal calibration cycles to include at the beginning of the mux scan, i.e. num_mux_cycle explained above.
• CV22 - Number of simultaneously-acquired bands (mux factor): how many slices to acquire simultaneously.
• CV23 - Slice separation (mm): for a contiguous set of slices, this must equal (slice_thickness + gap) * num_slices. (Note that setting this to zero will automatically compute the separation needed for contiguous slices.)
• CV24 - CAIPI z-blip flag: enables CAIPI z-blips (see Setsompop et. al. 2012)
• CV16 - CAIPI FOV shift: controls the CAIPI field shift. In general, this should be set to 3, or maybe 4 for higher mux factors. (Note that this CV only shows up when CV24 is set to 1.)

SMS DWI

Diffusion Spectrum Imaging (DSI) (Magn. Reson. Med., 2005, 54: 1377–1386) and multi-shell diffusion (Magn. Reson. Med., 2013, 69: 1534–1540) scans can be realized by designing gradient tables that specify direction and amplitude of the b-vectors. We set up a several customized gradient tables that are optimized for DTI, HARDI, 2 or 3-shell diffusion scans.

Consult with us if you would like to set up your own diffusion gradient scheme.

HARDI

• DTI pe0 (pe1) 80dir 1.7mm (5:50): SMS factor 3, axial slices, 1.7mm^3 voxel size, number of muxed slices 21 (63 unmuxed slices, 10.7 cm), b-value 2500, 80 diffusion directions, 8 b=0 images. The pe1 scan is a repeat of the pe0 scan with reversed phase encoding direction.
• DTI pe0 (pe1) 80dir 2mm (4:45): SMS factor 3, axial slices, 2mm^3 voxel size, number of muxed slices 23 (69 unmuxed slices, 13.8 cm), b-value 2500, 80 diffusion directions, 8 b=0 images. The pe1 scan is a repeat of the pe0 scan with reversed phase encoding direction.
• DTI pe0 (pe1) 96dir 2mm (5:50): SMS factor 3, axial slices, 2mm^3 voxel size, number of muxed slices 23 (69 unmuxed slices, 13.8 cm), b-value 3000, 96 diffusion directions, 10 b=0 images. The pe1 scan is a repeat of the pe0 scan with reversed phase encoding direction.
• DTI 80dir arc2 2mm (4:05): SMS factor 3, 2x in-plane acceleration, axial slices, 2mm^3 voxel size, number of muxed slices 23 (69 unmuxed slices, 13.8 cm), b-value 2500, 80 diffusion directions, 8 b=0 images. This scan uses 2x in-plane acceleration so there's no need for a separate pe1 scan for distortion correction.

Multi-shell diffusion

• Multi-shell with different TE: this set of scans are set up for a series of different b-values and number of diffusion directions to collect multi-shell diffusion data. The optimized TE for different b-values are different, so to get the best SNR possible we set up separate scans for each b-value. All scans have the same following acquisition parameters: SMS factor 3, axial slices, 2.4mm^3 voxel size, number of muxed slices 21 (63 unmuxed slices, 15.1 cm).
  • DTI pe1 b1k 2.4mm (2:49): b-value 1000, 64 diffusion directions, 6 b=0 images.
  • DTI pe0 b1k 2.4mm (2:47): b-value 1000, 64 diffusion directions, 5 b=0 images.
  • DTI pe1 b3k 2.4mm (3:20): b-value 3000, 64 diffusion directions, 6 b=0 images.
  • DTI pe0 b3k 2.4mm (3:18): b-value 3000, 64 diffusion directions, 5 b=0 images.
  • DTI pe1 b5k 2.4mm (0:39): b-value 5000, 6 diffusion directions, 1 b=0 images.
  • DTI pe0 b5k 2.4mm (6:46): b-value 5000, 128 diffusion directions, 10 b=0 images.
  • DTI pe1 b7k 2.4mm (0:42): b-value 7000, 6 diffusion directions, 1 b=0 images.
  • DTI pe0 b7k 2.4mm (7:15): b-value 7000, 128 diffusion directions, 10 b=0 images.
• Multi-shell with single TE: the following two scans are simplified multi-shell diffusion scans where diffusion directions of different shells are acquired in a single scan with the same TE (optimized for the highest b-value).
  • DTI g103 b2k 2-shell (4:50): 2-shell with 9 b=0 images, 30 directions at b=700, 64 directions at b=2000. SMS factor 3, axial slices, 2x in-plane acceleration, voxel size 2mm^3, number of muxed slices 25 (75 unmuxed slices, 15 cm)
  • DTI g150 b3k 3-shell (6:15): 3-shell with 10 b=0 images, 30 direction at b=1000, 45 direction at b=2000, 65 direction at b=3000. SMS factor 3, axial slices, 2x in-plane acceleration, voxel size 2mm^3, number of muxed slices 21 (63 unmuxed slices, 12.6 cm)

Reconstruction

Reconstruction on Flywheel

The mux epi recon is started automatically on Flywheel at the end of a mux scan. It can take anywhere from a few minutes to an hour per scan, depending on the resolution and number of TRs acquired. The default reconstruction method is 1D-GRAPPA for fMRI scans and split-slice-GRAPPA plus SENSE1 coil combination for diffusion scans. For a more detailed discussion on the difference between these methods, please also refer to this CNI blog post. The reconstruction on Flywheel also generates a metadata file (.json) for each mux scan. The file contains many useful parameters about the scan, such as acquisition parameters, timing information, calibration method, etc. Many of these information are not available in the nifti header, so it is often very useful to refer to this file.

Useful Parameters

Timing parameters in JSON

• mux_cal_duration: duration of the internal calibration (sec)
• num_acquired_timepoints: total number of acquired TRs
• num_mux_cal_cycle: number of reference image repetitions, i.e. num_mux_cycle
• num_mux_cal_timepoints: number of TRs used for internal calibration
• num_mux_cal_volumes_in_nifti: number of calibration volumes saved in the reconstructed nifti
• num_usable_volumes_in_nifti: number of usable volume in nifti
• num_volumes_in_nifti: total number of volumes in nifti, i.e. num_usable_volumes_in_nifti + num_mux_cal_volumes_in_nifti

Translating to BIDS parameters

• TotalReadoutTime: EffectiveEchoSpacing * (ReconMatrixPE - 1). EffectiveEchoSpacing is the "ec" in the "descrip" field of NIFTI header. You can get it by running "fslhd NIFTI | grep ec=". ReconMatrixPE is the y resolution. You can also look for "effective_echo_spacing" and "acquisition_matrix" in the json file for these parameters.
• MultibandAccelerationFactor: "num_bands" in the json file.
• DwellTime: 1/bandwidth. Bandwidth is the "pixel_bandwidth" in the json file, in kHz.

Slice acquisition order

The multiband sequences acquire slices in either ascending or descending interleaved order, depending on which direction the slices were prescribed (the nifti header should indicate whether it's ascending or descending-- see the nifti faq for a description of the two orders). And, of course, for mux>1, multiple slices are acquired simultaneously. E.g., for a mux 3 scan with 14 muxed slices (that will result in 14*3=42 unmuxed slices), the acq. order is [1,3,5,7,9,11,13,2,4,6,8,10,12,14]. And mux slice 1 will get reconstructed as unmuxed slices [1,15,29], which are all acquired simultaneously. Likewise, mux slice 2 will get deconvolved into unmuxed slices [2,16,30]. And so on.

Here's an algorithm to compute the relative slice acquisition time within the TR:

   mux = 3
   nslices = 5
   tr = 1.0
   mux_slice_acq_order = range(0,nslices,2) + range(1,nslices,2)
   mux_slice_acq_time = [float(s)/nslices*tr for s in xrange(nslices)]
   unmux_slice_acq_order = [nslices*m+s for m in xrange(mux) for s in mux_slice_acq_order]
   slice_time = {slice_num:slice_time for slice_num,slice_time in zip(unmux_slice_acq_order, mux_slice_acq_time*3)}
   for slice,acq in sorted(slice_time.items()):
       print "    slice %02d acquired at time %.3f sec" % (slice+1,acq)

   slice 01 acquired at time 0.000 sec
   slice 03 acquired at time 0.200 sec
   slice 05 acquired at time 0.400 sec
   slice 02 acquired at time 0.600 sec
   slice 04 acquired at time 0.800 sec
   slice 06 acquired at time 0.000 sec
   slice 08 acquired at time 0.200 sec
   slice 10 acquired at time 0.400 sec
   slice 07 acquired at time 0.600 sec
   slice 09 acquired at time 0.800 sec
   slice 11 acquired at time 0.000 sec
   slice 13 acquired at time 0.200 sec
   slice 15 acquired at time 0.400 sec
   slice 12 acquired at time 0.600 sec
   slice 14 acquired at time 0.800 sec

To generate an FSL slice time file:

   mux = 3
   nslices = 15
   descending = False
   filename = 'slicetime.txt'
   mux_slice_acq_order = range(0,nslices,2) + range(1,nslices,2)
   if descending:
       mux_slice_acq_order = mux_slice_acq_order[::-1]
   unmux_slice_acq_rel_time = [float(s)/nslices - 0.5 for s in xrange(nslices)] * mux
   unmux_slice_acq_order = [nslices*m+s for m in xrange(mux) for s in mux_slice_acq_order]
   slice_time = {slice_num:slice_time for slice_num,slice_time in zip(unmux_slice_acq_order, unmux_slice_acq_rel_time)}
   with open(filename,'w') as fp:
       for slice,acq in sorted(slice_time.items()):
           fp.write("%.4f\n" % acq)

Sharing

Pulse sequence

We are happy to share our SMS-EPI pulse sequences (and associated reconstruction code) with other GE research sites running a compatible platform (e.g., MR750 and related variants). We have two versions, one for BOLD fMRI and one for diffusion imaging. To request access to the pulse sequence (PSD) binaries and associated image reconstruction software, please send a note to Bob Dougherty (bobd@stanford.edu). You will also need to complete the sharing agreement and verify that you have an active research agreement with GE. Once the paperwork is out of the way, we will share the PSDs and recon code via a private Github repo, so you will need to send us your Github username and we will grant you access to the current code and future updates. The current versions of the PSDs run on both DV24 and DV25.

Please note that this is very much "research grade" software and as such, it is not a polished "turn-key" package. You will need someone with experience installing research PSDs on a GE platform to get the sequences working, and will also need to set a mechanism for collecting the raw data files from the scanner and configuring the off-line reconstruction software on a Linux machine.

Offline reconstruction

The standalone reconstruction package is shared via a private GitHub repository. Our recon should run on any modern linux system. However, it may be pretty slow (hours per fMRI scan), unless you have a beefy Intel CPU (e.g., a 12-core Sandy Bridge will get you down to 30 minutes or so) and/or good GPU cards (a Sandy Bridge CPU with three Nvidia Titan-X cards will get the recon down to a few minutes and allow it to run in real-time). We have spec'ed a linux workstation that will cost under $8k to build and can do real-time recons. All components are available on Amazon and NewEgg:

• ASRock X99 Extreme6 LGA 2011-v3 Intel X99 SATA 6Gb/s USB 3.0 ATX Intel Motherboard
• Intel Xeon E5-2697 v3 Haswell 2.6 GHz 35MB L3 Cache LGA 2011-3 145W BX80644E52697V3 Server Processor (14 core)
• G.SKILL Ripjaws 4 Series 64GB (8 x 8GB) 288-Pin DDR4 SDRAM DDR4 2800 (PC4 22400) Intel X99 Desktop Memory
• Noctua NH-U14S 140x150x25 (NF-A15 PWM) SSO2-Bearing CPU Cooler
• EVGA SuperNOVA 1600 G2 80+ GOLD, 1600W Power Supply
• 3 X EVGA GeForce GTX TITAN X 12GB (superclocked)
• Samsung 850 Pro 1 TB 2.5-Inch SATA III Internal SSD
• 2.5” to 3.5” mounting bracket
• Corsair Obsidian Series 750D Performance Full Tower Case
• keyboard & mouse

Notes:

• The 14-core CPU (currently $2700) could probably be replaced by a 10 or 12 core 2011-v3 CPU to reduce the cost (e.g., Intel Xeon E5-2680 v3 Haswell 2.5 GHz 12 x 256KB L2 Cache 30 MB L3 Cache LGA 2011-3120W BX80644E52680V3 Server Processor). We haven’t tested this, but we suspect such a processor would be adequate for real-time reconstruction based on the spare capacity of the 14-core CPU when running a recon.
• The power supply is also bigger than needed (we got this one to allow for upgrade to a dual-processor motherboard and a fourth GPU card, if needed). Total system power should be well under 1000 Watts at full load, so a 1000W supply would be fine (e.g., EVGA 120-G2-1000-XR 80 PLUS GOLD 1000 W).
• A smaller case could also be used (again, we got this one to allow for an upgrade with an an extended ATX motherboard). Just be sure there is enough room between the motherboard and the power supply to intall the third GPU card.
• Do NOT populate the Ultra M.2 socket on this ASRock motherboard, as doing so disables the third PCIe 3.0 x16 slot. (This is mentioned in the manual, but we didn't read that until after many hours of frustrated debugging!)