mrHARDIflow : A pipeline tailored for the preprocessing and analysis of Multi-Resolution High Angular diffusion MRI and its application to a variability study of the PRIME-DE database

The validation of advanced methods in diffusion MRI requires finer acquisition resolutions, which is hard to acquire with decent Signal-to-Noise Ratio (SNR) in humans. The use of Non-Human Primates (NHP) and anaesthesia is key to unlock valid microstructural maps, but tools must be adapted and configured finely for them to work well. Here, we propose a novel processing pipeline implemented in Nextflow, designed for robustness and scalability, in a modular fashion to allow for maintainability and a high level of customization and parametrization, tailored for the analysis of diffusion data acquired on multiple spatial resolutions. Modules of processes and workflows were implemented upon cutting edge and state-of-the-art MRI processing technologies and diffusion modelling algorithms, namely Diffusion Tensor Imaging (DTI), Constrained Spherical Deconvolution (CSD) and DIstribution of Anisotropic MicrOstructural eNvironments in Diffusion-compartment imaging (DIAMOND), a multi-tensor distribution estimator. Using our pipeline, we provide an in-depth study of the variability of diffusion models and measurements computed on 32 subjects from 3 sites of the PRIME-DE, a database containing anatomical (T1, T2), functional (fMRI) and diffusion (DWI) imaging of Non-Human Primate (NHP). Together, they offer images acquired over a range of different spatial resolutions, using single-shell and multi-shell b-value gradient samplings, on multiple scanner vendors, that present artifacts at different level of importance. We also perform a reproducibility study of DTI, CSD and DIAMOND measurements outputed by the pipeline, using the Aix-Marseilles site, to ensure our implementation has minimal impact on their variability. We observe very high reproducibility from a majority of diffusion measurements, only gamma distribution parameters computed on the DIAMOND model display a less reproducible behaviour. This should be taken into consideration when future applications are performed. We also show that even if promising, the PRIME-DE diffusion data exhibits a great level of variability and its usage should be done with care to prevent instilling uncertainty in statistical analyses.

as simple as possible for the end-user. 87 To test the pipeline, we processed 32 datasets from 3 sites of the PRIME-DE [Neff (2019)   . 4 scopes composing the Nextflow library and how they interface with each others. A pipeline definition (Pipeline X for example) is placed at the root of the library. To load the data, it uses one or more convention from the input scope (convention 1 and the bids convention). Then, to process it, it calls utilities in either or both the process scope (ANTs registration, Concatenate images and others) and the workflow scope (execute Topup, then register T1 to DWI). A call to an utility in the workflow scope will trigger calls to processes in the process scope they include to process the data (execute Topup firsts calls a concatenation of two opposite phase encoding DWI and extracts their b0 volumes, to pass the result to Topup for phase distortions correction). The input, workflow and process scopes form the core of the Nextflow library. Their content can be used by any pipeline defined in the pipeline scope, allowing high levels of code reusability. prescribed to process raw anatomical and diffusion weighted images, their execution is by default optional 154 and can be individually turned on and off using the Nextflow configuration file, for instance to allow the 155 computation of the models and measures on already preprocessed data or to shorten the execution time by 156 turning off steps deemed unnecessary after input data quality control.  Figure 2. DWI preprocessing steps. (Background denoising) Input 4D DWI from both phase encoding directions are denoised to remove background noise contamination. (b0 to b0 normalization) The signal is normalized so all b0 volumes in both phase encoding directions exhibit the same mean value. Each 3D diffusion volume is normalized using as reference either the b0 volume coming before it, after it or a linear interpolation of both. (Topup) b0 volumes from both phase encoding directions are concatenated and a deformation field is estimated. (Apply Topup) This field is applied on the forward phase encoded DWI to produce an undistorted 4D DWI. (Extract b0 + BET) From the undistorted 4D DWI is extracted a mean b0 volume, upon which is computed a brain mask. All voxels outside the brain mask are set to 0 in the mean b0 volume and the undistorted 4D DWI. ( 253 Once segmentation is computed, the upsampled T1 image is registered to the upsampled DWI space using     To evaluate correctly the potential variations in metrics between the several subjects of the PRIME-DE 299 databases, it is a requirement that the processing carried on the data be as reproducible as possible. To do 300 so, all algorithms used by the pipeline that require the execution of a random process have seen their seed 301 enforced to the same number across all subjects and sites. To evaluate the reproducibility, we did a test-test 302 analysis, using the 4 subjects from the aix-marseille database, on which the pipeline was run 3 separate

Registration of T1 to DWI space
Here, I ss refers to a modality, mask or measurement to evaluate for a specific session ss belonging to a 306 subject s,Î s is the mean image for a specific subject s taken over all its sessions andÎ is the mean image 307 taken over all sessions and all subjects.  with either isotropic or rectangular voxels, and anatomical images at voxels sizes ranging from 0.3 mm 3 to 319 0.8 mm 3 isotropic.  and manually fixed to prevent exclusion of brain voxels and inclusion of skull and eyes. The inclusion of 331 this technology in the pipeline was considered, but was ultimately rejected due to the weight in memory of 332 its dependencies, practically doubling the size of the container required to run the pipeline. Moreover, the 333 lack of validation on databases other than the PRIME-DE, which was used to generate the segmentation models, and its insertion near the beginning of the processing chain was deemed too risky, since its failure 335 at providing a good brain mask would impede on most of all subsequent steps executed.
CNR was computed on the b0 by comparing white matter with gray matter regions using equation 4.

CNR(I, WM mask
We also performed a variability study in tissue maps (WM and GM) of several metrics extracted from 359 the diffusion models available in the pipeline. Using the average over sessions of subjects, we analyzed CV intra-vendor ({Ŝ 1,D 1 , . . . ,Ŝ n,D 1 , . . . ,Ŝ 1,Dv , . . . ,    Table 2. Reproducibility measures on images generated by the pipeline.  variable at 6% and fractional anisotropy scores at 10%, 2% bellow the fascicle based measure (max fFA).

421
Results shows that performance from those measure is approximately two times more variable in GM.

422
Free-water fraction is the most variable measure estimated at 28% in WM and 21% in GM. For gamma 423 distribution parameters, values range widely depending on the parameter and on the tissue. For WM, kappa 424 and kappaAD are the less variable at 7% and 10%, and hei and heiAD score at 33% and 18% respectively.

425
For GM, the more variable parameter is kappaAD at 32%, followed by kappa and hei at around 17% and 426 heiAD at 13%.

503
Another trend observed in our study hinting to the need of guidelines and protocols is the high mean 504 intra-subject variability observed for UC-Davis. While analyzing this site, we noticed a strong discrepancy 505 of diffusion measurements between the two sessions images were acquired at. This behaviour could be the 506 result of equipment upgrades, software or hardware, as well as poor sequence configuration, resulting in an 507 ill-conditioning of the tissues required for diffusion acquisition. The only quantities were this behaviour 508 was not detected are fractional anisotropies (FA and max fFA) and total apparent fiber density (AFDt).

509
Nevertheless, this makes the usage of separate sessions from UC-Davis impossible in a study of diffusion 510 measurements as the variability would erase any statistical power of an analysis. We hope in the future, 511 good quality control could be integrated in the scanners to prevent or correct the acquisition of images 512 presenting this kind of behaviour.

LIMITATIONS AND FUTURE WORK
In the current version of our pipeline, we implemented a robust and efficient preprocessing chain of 514 state-of-the-art algorithms to correct for well know artifacts affecting the MRI signal, such as background 515 noise, gibbs ringing, susceptibility distortions, motion, signal dropout and intensity disparities within 516 tissues. We are well aware that our technique is not universal and that some use-cases may require the 517 execution of different algorithms than the one we considered. MRI preprocessing still need to be thoroughly 518 validated; the effect of applying specific techniques and algorithms, as well as their ordering relative to one 519 another, requires to be quantified in order to attain an adequate level of consensus in the community. We 520 designed our pipeline and its modular framework taking those facts into account, to ease the process of 521 swapping and upgrading its different steps, in order to fit specific study cases. We also plan to release in 522 the near future a version of the pipeline integrating alternative denoising techniques that are popular in the 523 diffusion MRI research community based on recommendations from the ISMRM diffusion study consensus 524 effort, once their work becomes available. Inasmuch, we intent to include brain extraction workflows, 525 for both anatomical and diffusion weighted images, adapted to the morphology of input subjects. Good 526 solutions have been developed that display good resiliency on human subjects, but their extension to other 527 primates and to small animals is still to be validated and can lead to segmentation lacking specificity. As 528 the field continues to evolve and those issues become addressed, we will consider the addition of those 529 solutions in our pipeline.

540
It is part of future work to adress better computing resources management, there is room for improvement.

541
The Nextflow scripting language offers extensive capabilities to define per process resources prerequisites.

542
Using them, we were able to define clearly the number of processing cores (CPU) and need for hardware needs could lead to a fewer number of processes executed at the same time, which would impede the 549 parallelization capabilities of the pipeline and result in a waste of computing resources. We are aware that 550 this behaviour could lead to execution problems, more so when processing high resolution images that 551 spans multiple gigabytes on disk. To mitigate the effects of that limitation, we opted for a retry strategy 552 for processes failing for that reason, combined with a restriction of the number of parallel executions of 553 memory hungry processes.

CONCLUSION
In this work, we built a reliable and reproducible pipeline, able to tackle the task of preprocessing and 555 computing diffusion models on multi-resolution diffusion MRI data. We also presented our modular library 556 of Nextflow processes and workflows using state-of-the-art and cutting edge DWI processing technologies, 557 designed to be easily upgradable and adaptable. We used our pipeline to analyze the variability of 3 sites 558 (32 subjects) included in the PRIME-DE database presenting good quality DWI data. We showed that 559 even if promising, that data exhibited a great level of variability and its usage should be done with care to 560 prevent instilling uncertainty in statistical analyses. This is a good example of the effects of the lack of 561 consensus in the diffusion MRI field when it comes to specifying guidelines for acquisition procedures.

562
The hurdle of defining a gold standard in dMRI makes this problem even more difficult to overpass, but 563 a level of standardization and harmonization must be attained in order to reduce variability in imaging 564 between subjects, scanner and sites as much as possible. Without it, conjugate efforts like the PRIME-DE 565 become less relevant for the study of the brain and its properties. It is through conjoint endeavours such 566 as the consensus effort from the ISMRM diffusion study group and others, and by the development of 567 reproducible, robust and maintainable software, that we as a community will make of diffusion MRI a 568 reliable modality.

DATA AVAILABILITY STATEMENT
The pipeline developed for this study and the nextflow modules library are available online at