Multicenter fixel based analyses: is this possible?

Dear MRtrix experts,

I was wondering if it is, in principle, possible to perform multicenter fixel-based analysis? And if yes, do I need to perform some extra steps?

I am working on a dataset of children at high-familial risk of developing severe mental illness.
Participants were scanned at two different sites with respectively a Siemens Prisma fit (64-channel head coil) or Skyra scanner (32-channel head coil). See below for some details on preprocessing, wm response functions retrieved for Prisma and Skyra data, and DWI sequence parameters.

I notice that the first 4 Prisma values are approximately a factor 3 larger than those for Skyra.

I am looking forward to your response
Best wishes

Data were preprocessed with qsiprep1.6 :

  1. Denoising consisted of four steps:
    • MP-PCA (Marchenko-Pastur distribution – Principal Component Analysis) denoising, as implemented in MRtrix3’s dwidenoise (Veraart et al., 2016), using a 5-voxel window.
    • Gibbs unringing using Mrtrix3’s mrdegibbs (Kellner et al., 2016).
    • B1 field inhomogeneity correction using dwibiascorrect from Mrtrix3 with the N4 algorithm (Tustison et al., 2010).
    • b=0 image-based intensity normalization
    • Head motion, eddy current, and susceptibility by movement distortion correction
  2. b=0 reference image creation
  3. Rigid 6-parameter (rotations and translations only) co-registration to the skull-stripped T1w image and DWI image resampling (1 mm isotropic) and gradient rotation.

I generated (and checked the fit) DWI masks using synthstrip-singularity (-b -2) on generated FA images

I used dhollander to retrieve wm, gm, and csf response functions:

dwi2response dhollander ${indata} \
-fslgrad ${inbvec} ${inbval} \
${OUT_DIR}/${bids_id}_dwi_${res}_response_wm.txt \
${OUT_DIR}/${bids_id}_dwi_${res}_response_gm.txt \
${OUT_DIR}/${bids_id}_dwi_${res}_response_csf.txt \
-mask ${newmask} \
-voxels ${bids_id}_dwi_${res}_response_voxel_selection.nii.gz \
-info \

Example: wm response functions:


Shells: 5,1000

command_history: /home/wimb/anaconda3/envs/mrtrix3/bin/amp2response dwi.mif voxels_sfwm.mif safe_vecs.mif response_sfwm.txt -shells ‘5,1000’ (version=3.0.4)

21307.7996795693 0 0 0 0 0
11488.0266075202 -4509.69162685159 978.658571954375 -137.330096590635 11.2965465825901 4.40098708965145


Shells: 5,999

command_history: /home/wimb/anaconda3/envs/mrtrix3/bin/amp2response dwi.mif voxels_sfwm.mif safe_vecs.mif response_sfwm.txt -shells ‘5,1000’ (version=3.0.4)

7053.85982998796 0 0 0 0 0
3723.35267227721 -1498.52515490046 328.62205129044 -54.0957827205385 7.7072439141182 6.4925266058699

DWI sequence parameters:
A diffusion-weighted imaging (DWI) sequence was setup for diffusion tensor imaging (DTI) and tractography and diffusion kurtosis imaging (DKI), using a pulse gradient spin echo (PGSE) echo planar imaging (EPI) sequence with uniformly distributed gradient directions. The DTI-tractography part consisted of 14 non-DW images (b = 0) and 80 DW images (b=1000 s/mm2); TR = 2541 ms; TE = 64 ms; flip angle = 90 degrees; matrix = 112 x 114; FOV = 200 x 204; phase partial Fourier = 6/8; echo spacing = 0.73 ms; bandwidth = 1540 Hz/Px; phase encoding = Anterior > Posterior, 75 transversal slices, multiband acceleration factor = 3; GRAPPA: acceleration factor = 2, reference lines = 24; 1.8 x 1.8 x 1.8 mm3 voxels.

Diffusion-weighted imaging (DWI) consisted of two sequences. One for diffusion tensor imaging (DTI) and tractography and one for diffusion kurtosis imaging (DKI), both using a pulse gradient spin echo (PGSE) echo planar imaging (EPI) sequence with uniformly distributed gradient directions. The DTI-tractography part consisted of 13 non-DW images (b = 0) and 57 DW images (b=1000 s/mm2); TR = 2656 ms; TE = 77.2 ms; flip angle = 90 degrees; matrix = 100 x 100; FOV = 200 x 200; phase partial Fourier = 6/8; echo spacing = 0.68 ms; bandwith = 1665 Hz/Px; phase encoding = Anterior > Posterior, 72 transversal slices, multiband acceleration factor = 3; GRAPPA: acceleration factor = 2, reference lines = 24; 2 x 2 x 2 mm3 voxels

1 Like


I did a multicenter fixel based analysis, described here: (code: GitHub - smeisler/Meisler_Reading_FBA: Code used in the Meisler and Gabrieli 2022 eLife paper on fixel based analyses of reading abilities.).

One extra thing you have to do is to average fiber response functions within each center, NOT across everyone. I additionally added a linear regressor for site in my final model to capture remaining variance, but there are probably better ways of accounting for site in the final model.


1 Like

Many thanks Steven, much appreciated!


Following on this, this year, it has been presented an extension of combat for the fixel-based analysis. I don’t know when it will be available. I hope this helps.

Best regards,



Hi Steven,

I just finalized a grant application and went through your article. A very comprehensive read.
I have a few questions I hope you can clarify.

  1. I might have missed it but after determining the response functions, averaging response functions separately for each site, generating the individual FODs, how did you deal with the different sites in the FOD template generation?
    It seems you included individuals from different sites, but it is not clear if you balanced across sites? I would think the latter is important

  2. Registration of individual FOD to template FOD
    Did you use the default settings of mrregister or did you for example employ the scheme used in the FOD template generation?

  3. Quality control of FOD warps. It seems you only warped and overlaid brain masks warped into template space to check it the individual FOD to template FOD warp was successful

I was a bit surprised by this, as I would have thought that warping and overlaying for example FA images would give a better idea, as this would visualize how well white matter (“tracts”) across subjects overlap.

Some time ago I warped individual FODs generated by tractseg 2.8 to a study FOD template ( based on tractseg FODs) including data from two sites (48 subjects balanced across sites, sex and group (not age, as subjects were all around 12 years of age)). Although, this might be not ideal with respect to response functions, I noticed when applying warps to FA data, that while the overall fit seemed ok (brains occupying the same space, a la warping brain masks), their was considerable variation in the fit of white metter. (this was true for mrregister default warps and warps using the same scheme as in the FOD template generation)

So my question is if my surprise is due to incorrect expectations from my side? i.e., are relative rough/ low dimensional non-linear warps sufficient for fixel analyses?

Many thanks in advance

Hi @WilliamFCB sorry for delayed response, hopefully this is still useful to you:

I did not consider balance sites when making the template (I focused more on balancing age and reading skill). Maybe I should have, but the template turned out fine.

I did this as described in the usual fixel based analysis pipeline, e.g.,

for_each * : mrregister IN/wmfod_norm.mif -mask1 IN/dwi_mask_upsampled.mif  \
../template/wmfod_template.mif -nl_warp IN/subject2template_warp.mif IN/template2subject_warp.mif

I did this because the final analysis mask is the intersection of all of the warped masks. Given the ~1000 images I had to quality control, this seemed like the most efficient way to do this, but there probably is a more comprehensive way to QC this stage.

I don’t think this is necessarily surprising, in that internal white matter shapes can vary substantially across individuals. The fixel based pipeline has the fixelcorrepsondance command to account for this, attempting to correspond fixels across subjects.


Hi Steven,

Many thanks for your reply.

I finished the processing some time ago. Just went back to the data after having to focus on other stuff. Basically, I used the default pipeline, with a few exceptions ( e.g., for example, using synthstrip singularity to optimise some masks). I experimented with warping, using the default mrregister settings and those used in the population template generation step, for warping individual wmFODs to the study wmFOD template. The template warp settings did not add that much, so I decided to go for the default. For the -fmls_peak_value in the fod2fixel step, I experimented with values between 0.03 and 0.13 in steps of 0.1. After initial evaluation, I went further with 0.06, 0.07 and 0.08. Thresholds 0.06 and 0.07 retrained a bit more crossing fibres. However, compared to 0.08, fixels extended more/too much into GM. Threshold 0.08 retains most crossing fibres and covers primarily WM. However, as this is subjective, I generated data for all 3 thresholds (including tractography etc) to perform post-hoc analyses on the 0.06 and 0.07 thresholded data.

Below some figures

figure 1. Mean of 276 FA images in study wmFOD template space

Figure 2a,b and c: Showing fixels with a threshold of 0.06 in colour on top of fixels with a threshold of 0.08 in black

figure 3
tractography with cutoff 0.06

figure 4
tractography with cutoff 0.08

1 Like

Hi all,

I’m having a similar problem, however, the data I’m using has ~100 subjects acquired through 6 acquisition sites, expanding from 35 to 5 subjects per site.

I wonder whether the site-specific average response function will benefit my case. As an alternative, I am also considering the ComBat application for the fixel data (the ISMRM abstract posted earlier; their paper with application I assume: Still, my cases are all patients (though two different groups) and no healthy controls are available.

Also, are there any differences on the fixel-level if these harmonization approaches are applied before any preprocessing, i.e., the raw data rather than on the post-processed data?

I am looking forward to your comments.