ACID: A comprehensive toolbox for image processing and modeling of brain, spinal cord, and ex vivo diffusion MRI data

Abstract Diffusion MRI (dMRI) has become a crucial imaging technique in the field of neuroscience, with a growing number of clinical applications. Although most studies still focus on the brain, there is a growing interest in utilizing dMRI to investigate the healthy or injured spinal cord. The past decade has also seen the development of biophysical models that link MR-based diffusion measures to underlying microscopic tissue characteristics, which necessitates validation through ex vivo dMRI measurements. Building upon 13 years of research and development, we present an open-source, MATLAB-based academic software toolkit dubbed ACID: A Comprehensive Toolbox for Image Processing and Modeling of Brain, Spinal Cord, and Ex Vivo Diffusion MRI Data. ACID is an extension to the Statistical Parametric Mapping (SPM) software, designed to process and model dMRI data of the brain, spinal cord, and ex vivo specimens by incorporating state-of-the-art artifact correction tools, diffusion and kurtosis tensor imaging, and biophysical models that enable the estimation of microstructural properties in white matter. Additionally, the software includes an array of linear and nonlinear fitting algorithms for accurate diffusion parameter estimation. By adhering to the Brain Imaging Data Structure (BIDS) data organization principles, ACID facilitates standardized analysis, ensures compatibility with other BIDS-compliant software, and aligns with the growing availability of large databases utilizing the BIDS format. Furthermore, being integrated into the popular SPM framework, ACID benefits from a wide range of segmentation, spatial processing, and statistical analysis tools as well as a large and growing number of SPM extensions. As such, this comprehensive toolbox covers the entire processing chain from raw DICOM data to group-level statistics, all within a single software package.


INTRODUCTION
Diffusion MRI (dMRI) exploits the self-diffusion of water molecules to produce images that are sensitive to tissue microstructure by measuring the diffusion along various spatial directions ( Callaghan et al., 1988;Le Bihan et al., 1988;Stejskal & Tanner, 1965).dMRI has been applied to study a number of phenomena including normal brain development ( Dubois et al., 2014;Miller et al., 2002), aging ( Draganski et al., 2011;Sullivan et al., 2010), training-induced plasticity ( Scholz et al., 2009), and monitoring progression of and recovery from neurological diseases ( Farbota et al., 2012;Meinzer et al., 2010).Clinical applications of dMRI include the diagnosis of ischemic stroke ( Urbach et al., 2000), multiple sclerosis ( Horsfield et al., 1996), cancer and metastases ( Gerstner & Sorensen, 2011), and surgical planning of brain tumors ( Chun et al., 2005).Although the vast majority of dMRI applications has focused on the brain, there is a growing interest in spinal cord dMRI, as researchers seek sensitive and predictive markers of spinal cord white matter damage ( Cohen et al., 2017;Martin et al., 2016).Furthermore, an increasing number of studies utilize dMRI on ex vivo specimens for comparative analysis with other imaging modalities, such as electron microscopy ( Barazany et al., 2009;Kelm et al., 2016;Papazoglou et al., 2024).
While the majority of toolboxes have been designed for brain dMRI, ACID has introduced several features and utilities that make it particularly suitable for spinal cord and ex vivo dMRI as well.Specifically, ACID addresses the higher level and different nature of artifacts in spinal cord dMRI ( Barker, 2001;Stroman et al., 2014), and the highly variable geometry and diffusion properties in ex vivo dMRI (see Sébille et al., 2019 for a list of ex vivo/ postmortem dMRI studies).Although there are some software options available for processing spinal cord images, most notably the Spinal Cord Toolbox ( De Leener et al., 2017), these tools lack comprehensive artifact correction and biophysical modeling capabilities for estimation of dMRI-based metrics related to microscopic tissue properties.Biophysical modeling estimates microstructural properties, such as axonal water fraction and orientation dispersion, as aggregated measures on the voxel level, providing greater specificity than standard diffusion tensor (DTI) or diffusion kurtosis imaging (DKI).Toolboxes dedicated for biophysical modeling of the dMRI signal, such as the NODDI ( Zhang et al., 2012) or SMI toolbox ( Coelho et al., 2022), typically do not include a comprehensive processing pipeline to correct for artifacts in dMRI data.In addition, to date, only a few of the dMRI toolboxes support the Brain Imaging Data Structure (BIDS; Gorgolewski et al., 2016) standard for organizing and annotating raw and processed dMRI data.The lack of standardization not only complicates the sharing and aggregation of processed dMRI data but also the application of automated image analysis tools designed for big data, such as machine learning techniques.Over the past two decades, tens of thousands of dMRI datasets have been made openly available in large neuroimaging databases (e.g., HCP ( Van Essen et al., 2013) and the UK Biobank ( Littlejohns et al., 2020)), underscoring the importance of consistent data storage practices.
Building upon 13 years of research and development, we introduce an open-source MATLAB-based extension to the Statistical Parametric Mapping (SPM) software, the ACID toolbox: A Comprehensive Toolbox for Image Processing and Modeling of Brain, Spinal Cord, and Ex Vivo Diffusion MRI Data.ACID was initially developed as a collection of artifact correction tools but has now been extended to a comprehensive toolbox for processing and modeling of dMRI data.In particular, ACID offers (i) stateof-the-art image processing tools as well as (ii) DTI, DKI, and white matter biophysical model parameter estimation methods optimized for brain, spinal cord, and ex vivo dMRI data.Additionally, (iii) ACID adheres to the BIDS standard for organizing the output, making the processed images compliant with other BIDS software and facilitating data sharing.Finally, (iv) ACID is embedded in the SPM framework to benefit from its established functions including spatial processing tools and statistical inference schemes.ACID tools can be combined with other SPM functions to create pipelines in the SPM batch system, which offers an all-in-one software solution from conversion of DICOM data to statistical group analysis.ACID also benefits from a large and growing number of SPM extensions.For example, ACID can be combined with the SPM12-based hMRI toolbox ( Tabelow et al., 2019) to perform multicontrast analysis of dMRI and other quantitative MRI data, such as relaxation rates, acquired from the same subject, all within a single pipeline.Many of the methods used in the ACID toolbox have already been published in the scientific dMRI literature (Table 1).In this paper, we detail the design and function of the ACID modules and provide guidance on their optimal combination for brain, spinal cord, and ex vivo applications.

Overview
The ACID toolbox is a comprehensive toolbox for processing and analyzing dMRI data, built upon the following four pillars: (1) preprocessing of dMRI data (Pre-processing module), (2) physical models of the diffusion signal (Diffusion tensor/kurtosis imaging module), (3) white matter biophysical models of the diffusion signal (Biophysical models module), and (4) additional features referred to as Utilities.The Pre-processing module consists of state-ofthe-art methods for retrospective correction of the dMRI data.The Diffusion tensor/kurtosis imaging module contains tensor and kurtosis models that can be applied to dMRI data from various tissues or samples, including gray and white matter, as well as diffusion phantoms ( Woletz et al., 2024).In contrast, the Biophysical models module can only be applied to samples that fall within their validity ranges (see Section 4.2.2).The Utilities module contains various useful tools, including masking and noise estimation.The ACID toolbox follows the BIDS convention and enables the seamless integration of external tools into its processing pipeline in a modular fashion (External tools module).More details about the implementation and organization of ACID are provided in Appendix A.

Preprocessing
In this section, we provide brief descriptions of each artifact correction tool currently implemented in ACID.For detailed recommendations on various dMRI datasets (in vivo brain, in vivo spinal cord, ex vivo/postmortem), refer to Sections 3.2 and 4.1, as well as Table 5.

Eddy-current and motion correction (ECMOCO)
ACID uses the eddy-current and motion correction (ECMOCO) algorithm ( Mohammadi et al., 2010) to correct for spatial misalignments that may occur between dMRI volumes.These misalignments can be caused by motion and eddy currents induced by the rapidly varying field of the diffusion-sensitizing gradients ( Jezzard et al., 1998), which may lead to biased diffusion estimates ( Mohammadi, Freund, et al., 2013).ECMOCO aligns all source volumes to a target volume using a coregistration algorithm with an affine transformation ( Friston & Ashburner, 1997) implemented in the SPM function spm_coreg.It was previously shown that the robustness of registration can be increased by separately registering diffusion-weighted (DW) and nondiffusion-weighted (b0) volumes to their corresponding target volumes ( Mohammadi, Carey, et al., 2015).ECMOCO features the multitarget registration mode, where source volumes from each diffusion shell (b-value) are coregistered to their shell-specific target volume (Appendix Fig. B1).ECMOCO rotates the b-vectors by the obtained rotational parameters; these rotated b-vectors can be passed on to subsequent processing steps.Of note, the affine transformation of ECMOCO can only correct for first-order eddy-current displacements.The advantages and disadvantages of ECMOCO compared with other established tools, such as FSL eddy, are discussed in Section 4.1.
In spinal cord dMRI, eddy-current and motion correction is more challenging than in brain dMRI due to the considerably lower number of voxels and lower signalto-noise ratio (SNR), particularly in volumes with high b-values (>1000 s/mm 2 ) or with diffusion-sensitizing gradients parallel to the spinal cord.While movement of the brain can be considered approximately rigid, the spinal cord may experience varying degrees of displacement along the rostrocaudal axis caused by factors such as breathing, pulsation of the cerebrospinal fluid, or swallowing ( Yiannakas et al., 2012).To address this, we introduced slice-wise (2D) registration, which independently aligns each slice of the source volume to the corresponding slice of the target volume, thereby correcting for nonrigid, slice-dependent displacements ( Mohammadi, Freund, et al., 2013).For more details on ECMOCO, including other recently introduced features (initialized registration and exclusion mode), refer to Appendix B.

Adaptive denoising (msPOAS)
The Multi-shell Position-Orientation Adaptive Smoothing (msPOAS) is an iterative adaptive denoising algorithm designed to adaptively reduce noise-induced variance in dMRI data while preserving tissue boundaries, as illustrated in Figure 3 ( Becker et al., 2012( Becker et al., , 2014;;Tabelow et al., 2015).The algorithm adapts to the intensity values and their distance in both voxel space and the spherical space of diffusion directions, allowing smoothing only within spatially homogeneous areas of the DW images.One of the key advantages of msPOAS is its compatibility with all diffusion models as it operates on the raw dMRI data.Adjustable parameters include kstar (number of iterations that define the image smoothness), lambda (adaptation parameter that defines the strength of edge detection), kappa (initial ratio of the amount of smoothing between the local space of neighboring voxels and the spherical space of diffusion gradients), and ncoils (the effective number of receiver coils that contributed to the measured signal).To distinguish random fluctuations from structural differences, msPOAS requires an estimate of SNR, or equivalently the noise standard deviation (sigma).A higher kstar leads to greater smoothness within homogeneous image regions, while a larger lambda results in weaker adaptation and more blurring at tissue edges.The optimal kappa depends on the number of directions per shell, while ncoils should be the same as the value used for noise estimation.When using msPOAS, we recommend starting with the default parameters and the sigma estimated with the Noise estimation utility function (Table 2).In case of insufficient noise reduction, parameters should be adjusted according to Appendix D.

Rician bias correction
The voxel intensities of MRI magnitude images exhibit a Rician distribution in case of a single receiver coil ( Gudbjartsson & Patz, 1995) and a noncentral χ distribution in case of multiple receiver coils ( Aja-Fernández et al., 2014).When fitting diffusion signal models (Section 2.3), this distribution leads to a bias, known as the Rician bias, in the estimated tensor ( Basser & Pajevic, 2000;Gudbjartsson & Patz, 1995;Jones & Basser, 2004) and kurtosis parameters ( Veraart et al., 2011;Veraart, Rajan, et al., 2013), as well as in biophysical parameter estimates ( M. Andersson et al., 2022;Fan et al., 2020;Howard et al., 2022).This Rician bias is particularly relevant in low SNR situations ( Polzehl & Tabelow, 2016).Two approaches of Rician bias correction (RBC) are implemented in ACID.The M2 approach, introduced in Miller & Joseph (1993), and later extended to multichannel receiver coil ( André et al., 2014), operates on the dMRI data and uses the second moment of the noncentral χ distribution of the measured intensities and noise estimates to estimate the true voxel intensities.The second approach modifies the parameter estimation by considering the noncentral χ distribution to account for the Rician bias during model fitting ( Oeschger et al., 2023a).Note that the latter approach assumes uncorrected data, therefore, it must not be combined with the first method and is currently only available for nonlinear least squares fitting.Both methods require an estimate of the noise standard deviation, which can be obtained using either the standard or the repeated measures method within the Noise estimation utility function (Table 2).Details on noise estimation are available in Appendix C. In addition, ACID offers the Rician bias simulation utility function to determine the optimal RBC method for the dMRI dataset and SNR at hand (Table 2).An example of how RBC influences the estimation of biophysical parameters is illustrated in Appendix Figure F1.

Susceptibility artifact correction (HySCO)
Hyperelastic Susceptibility Artifact Correction (HySCO) is a technique used to correct for geometric distortions caused by susceptibility artifacts ( Ruthotto et al., 2012( Ruthotto et al., , 2013)).These artifacts can occur at interfaces between tissues with different magnetic susceptibilities, such as those found near paranasal sinuses, temporal bone, and vertebral bodies.To correct for these artifacts, HySCO estimates the bias field based on a reversed-gradient spin-echo echo planar imaging (EPI) acquisition scheme.This requires the acquisition of at least one image with identical acquisition parameters as the dMRI data but with reversed phase-encoding direction, also referred to as "blip-up" or "blip-down" acquisitions.The bias field map, estimated from the blip-up and blip-down images, is applied to the entire dMRI data to unwarp the geometric distortions (see Fig. 3 for examples).For datasets that include full blip-reversed acquisition, that is, each image was acquired with two phase-encoding directions (blip-up and blip-down), the reverse phase-encoded images can Table 2. List of ACID utility functions.

Function Description
Cropping Crops images to a smaller size for less storage space and faster processing.
Input: image(s) to crop, new matrix size, and voxel coordinates of the center of cropping.The center of cropping can also be selected manually through a pop-up window.
Output: cropped image(s) and the cropping parameters.Application: typically in spinal cord dMRI, where the spinal cord occupies a small portion of the image.

Resampling
Resamples images to the desired resolution.
Input: image(s) to be resampled, desired resolution, and type of interpolation (as defined in spm_slice_vol).
Application: for example, when performing voxel-wise arithmetic between two or more images with different resolutions (e.g., g-ratio mapping).

Slice-wise realignment
Enables manual translation and scaling of images along the x and y dimensions on a slice-by-slice basis, facilitated by intensity contour lines of the source image superimposed on the target image.
Input: image to be realigned, target image, and other images to which the realignment parameters are applied.
Output: realigned image(s) and the realignment parameters.Application: useful for realigning spinal cord images, where residual misalignments are often slice dependent.

Fusion
Merges two images with different field of views (FOV), such as a brain and a spinal cord image, into a single combined image (Fig. 5).
Input: two images to be merged and a target image (typically a structural image with a larger FOV).
Output: a combined image, resampled onto to the target image.The voxel intensity values in overlapping regions are the average of the intensity values in both images.Note that before merging the images, they must be in the correct spatial position; if necessary, image realignment can be performed using the SPM Realign or the Slice-wise realignment utility function.
Application: useful for merging a brain and a spinal cord image into a single image before applying a warping field obtained from a large-FOV structural image.

Create brain mask
Creates a binary brain mask by (i) segmenting the brain image into gray matter, white matter, and cerebrospinal fluid using SPM12's unified segmentation tool ( Ashburner & Friston, 2005), (ii) summing up the resulting probability maps, and (iii) thresholding it at a certain value (accessible through the script acid_local_defaults.m; default: 0.8).
Input: a single brain image or tissue probability maps for gray matter, white matter, and cerebrospinal fluid, and optionally a dMRI dataset to be masked.
Output: binary brain mask and optionally a masked dMRI dataset.Application: to restrict the estimation of DTI, DKI, and biophysical parameters to the brain for increased speed and efficiency.

Reliability masking
Aims to identify "unreliable" voxels, i.e., voxels irreversibly corrupted by artifacts.Reliability masks are generated by thresholding the root-mean-square model-fit error (rms(ε)) map ( David et al., 2017).Input: rms(ε) maps (output by tensor fitting methods with label: RMS-ERROR) and the desired threshold value.The optimal threshold can be determined using the Determine threshold submodule.
Output: a binary reliability mask.Application: Reliability masks can serve as binary masks in region-of-interest-based analyses.In principle, reliability masking as an outlier rejection technique is applicable after each model fitting method.It is particularly useful in situations where many data points are affected by outliers (often the case in spinal cord dMRI), which could otherwise lead to unstable tensor fits and inaccurate tensor estimates (see David et al., 2017, for examples).

DWI series browser
Enables browsing through the slices of the dMRI data for quality assessment.Slices with low SNR and/or artifacts can be identified and labeled.
Output: list of labeled slices.
Application: The saved labels can be used to inform ECMOCO about unreliable slices (see Exclusion mode in Appendix B).

DWI series movie
Enables simultaneous streaming of images from multiple dMRI datasets in video mode for quality assessment.
Input: a reference image and up to three dMRI datasets.
Output: a video file containing the image streams.
Application: useful for visual assessment of a single dMRI dataset or for comparing images before and after a specific processing step (e.g., ECMOCO).
be combined using the submodule HySCO: combine blip-up and blip-down images.

Diffusion signal models
The dependence of dMRI signal on the direction and strength of diffusion weighting is commonly described by mathematical models.Two of the most widely used models are DTI ( Basser et al., 1994) and DKI ( Hansen et al., 2016;Jensen et al., 2005).

Diffusion tensor imaging (DTI)
DTI describes the anisotropic water diffusion in the white matter by a diffusion tensor with six independent diffusion parameters.The eigenvalues of the tensor can be used to compute rotationally invariant DTI scalar metrics including fractional anisotropy (FA) and mean (MD), axial (AD), and radial diffusivities (RD).The interpretation of DTI assumes that the direction of axial diffusivity is aligned with the white matter tracts, which may not be the case in complex fiber geometry such as crossing or fanning fibers.ACID provides four algorithms to obtain the diffusion tensor (see Appendix E for details).Ordinary least squares (OLS) fits the tensor model by minimizing the sum of squared model-fit errors, while weighted least squares (WLS) minimizes the weighted sum of squared model-fit errors, accounting for the distortion of noise distribution in the linearized (logarithmic) data.Robust fitting is similar to WLS but factorizes the weights into three components to account for local and slice-specific artifacts as well, while also featuring Tikhonov regularization to handle illconditioned weighting matrices resulting from a high occurrence of outliers.Robust fitting is designed to down-weight outliers in the model fit, which can otherwise introduce a bias in the fitted model parameters Function Description

Noise estimation
Estimates the noise standard deviation (σ) in the dMRI data using either the standard or the repeated measures method.The standard method uses the formula σ ≈ S i 2 / 2Ln ( ) , where S i is the voxel intensity within a background mask defined outside the body, L is the number of voxels within the background mask, and n is the effective number of coil elements that contributed to the measured signal ( Constantinides et al., 1997).The repeated measures method uses the formula , where S(i, k ) is the voxel intensity at voxel i in the kth repeated image ( Dietrich et al., 2007).The standard deviation and mean operators are performed across the repetitions and voxels, respectively.The set of repeated images can be either the nondiffusion-weighted (b ≈ 0) or strongly diffusion-weighted (the highest b-value) images (see Appendix C for recommendations).
Input: the raw (unprocessed) dMRI dataset, a mask (standard method: background mask; repeated measures method: see Appendix C), n (for the standard method only), and b-values (for the repeated measures method only).
Output: a single σ (assuming a homogeneous variance).
Application: σ serves as input for msPOAS, Rician bias correction, and diffusion tensor imaging (for fitting methods WLS and robust fitting).

Rician bias simulation
Simulates diffusion-weighted MRI signals at specified SNR values in voxels within the brain white and gray matter.The simulated signals are corrected using the specified Rician bias correction (RBC) methods (for details, see Oeschger et al., 2023a).
Input: a voxel from a list of 27 predefined voxels, each with different diffusion and kurtosis tensor metrics 1 (for details, see Oeschger et al., 2023a), a list of SNR values, and the number of noise samples.
Output: a figure showing the distance between the estimated metric and the ground truth value for each RBC method.
Application: useful for computing the required SNR for DTI, DKI, and biophysical parameter estimation.

ROI analysis
Calculates the mean value within a specified region of interest (ROI).
Input: list of images and various types of ROIs including (i) global ROIs, applied to all images in the list, (ii) subject-specific ROIs, applied only to the corresponding image, and (iii) subject-specific reliability masks, again applied only to the corresponding image (see Reliability masking).
Output: an array containing the mean values within the specified ROIs per subject, ROI, and (optionally) slice.When multiple types of ROIs are specified, their intersection is applied.
Application: the function offers flexibility for a range of ROI-based analyses; for example, ROI-based analysis in the native space requires a set of subject-specific ROIs, while a single global mask is sufficient in the template space (with optional reliability masks in both cases).An example application including reliability masks can be found in David et al. (2017).

Diffusion kurtosis imaging (DKI)
DKI expands the diffusion tensor model by the kurtosis tensor, a fourth-order tensor with 15 independent parameters, which captures the effects of non-Gaussian water diffusion.From the 15 kurtosis parameters, several kurtosis metrics can be estimated including the mean (MK), axial (AK), and radial kurtosis (RK), as well as the mean (MW), axial (AW), and radial (RW) kurtosis tensor ( Tabesh et al., 2011) (Fig. 1).These metrics provide additional information about tissue complexity beyond what can be captured by diffusion tensor metrics alone.DKI requires the acquisition of a second diffusion shell with higher b-value (typically between 2000 and 2500 s/mm 2 ).ACID also includes the axisymmetric DKI model, a recent modification of DKI which reduces the parameter space to eight independent parameters by imposing the assumption of axisymmetrically distributed axons ( Hansen et al., 2016).Currently, ACID offers the OLS and NLLS algorithms for fitting the kurtosis tensor, and the NLLS algorithm for fitting the axisymmetric kurtosis tensor.Note that the diffusion tensor parameters from DKI might differ from standard DTI parameters.In particular, diffusivities (AD, MD, and RD) derived from the DTI model are often underestimated compared with those derived from the DKI model (referred to as kurtosis bias) ( Edwards et al., 2017).By incorporating higher-order moments of the diffusion signal, DKI can address kurtosis bias, resulting in more accurate diffusivity estimates (see Supplementary Fig. S3 in the Supplementary Material for a comparison of MD derived from DTI and DKI).

Biophysical models
Biophysical models separate the dMRI signal into distinct signal components from various tissue compartments, each with their own underlying assumptions.Biophysical models provide more specific and biologically interpretable metrics that are linked to tissue microstructure ( Jelescu et al., 2020).The application of biophysical models is often referred to as dMRI-based in vivo histology ( Mohammadi & Callaghan, 2021;Weiskopf et al., 2021) or microstructural dMRI ( Jelescu et al., 2020;Novikov, 2021;Novikov et al., 2019).In the following, we briefly describe the two white matter biophysical models currently implemented in ACID (WMTI-Watson and NODDI-DTI), while recommendations on their usage are provided in Section 4.    ( Novikov et al., 2018).We recommend using the plus branch (default in the toolbox), as in our experience, and also reported by others ( Jelescu et al., 2020;Jespersen et al., 2018), the minus branch is the biologically invalid solution.

NODDI-DTI
NODDI-DTI ( Edwards et al., 2017) is based on the neurite orientation dispersion and density imaging (NODDI) model ( Zhang et al., 2012).While NODDI is a three-compartment biophysical model with intra-and extra-axonal space, and cerebrospinal fluid compartments, NODDI-DTI assumes that the latter compartment can be neglected in normal appearing white matter.NODDI-DTI further For the spinal cord and ex vivo specimen, we refrained from masking for the white matter due to the difficulty of obtaining an accurate white matter mask.
assumes a fixed diffusivity of the intraneurite compartment (D a ).In our implementation, users can either choose from two fixed values tailored for in vivo (D a = 1.7•10 -3 mm2 /s) and ex vivo (D a = 0.6•10 -3 mm 2 /s) datasets, or select their own value.NODDI-DTI estimates the intraneurite (here: f) and extraneurite (1− f ) signal fraction, as well as the Watson concentration parameter κ and the orientation dispersion index (ODI), from the FA and MD maps.While WMTI-Watson requires specific multishell dMRI data for robust parameter estimation, NODDI-DTI parameters can also be obtained from single-shell DTI acquisitions.

Utilities
ACID utilizes SPM's utility functions, available under SPM -> Util in the SPM12 Batch Editor, for handling and manipulating NIfTI images.These functions include mathematical operations on single or multiple images, reorienting images, and concatenating 3D volumes and separating 4D volumes.Additionally, ACID provides its own set of utility functions for image manipulation, mask generation, quality assessment, and other related tasks (refer to Table 2 for more details).

External tools
ACID provides the option to integrate external tools from other packages, which can be accessed directly from the ACID graphical user interface (GUI) (External tools module), ensuring a seamless integration into ACID pipelines.
We included the following external tools in the current release: (i) FSL eddy 2 ( J. L. R.  ( Kellner et al., 2016); and (vii) the WMTI model (part of the DESIGNER toolbox) ( Fieremans et al., 2011).ACID also allows expert users to incorporate their own external tools into the toolbox, using the aforementioned examples as a template.

Output structure and BIDS naming convention
ACID supports the BIDS standard, while also being compatible with non-BIDS data.Following BIDS recommendations, ACID appends a label to the output filename's desc field, which indicates the applied processing step (refer to Table 3 for a list of labels used in the modules Preprocessing, Diffusion tensor/kurtosis imaging, and Biophysical models).For instance, after applying ECMOCO to sub01_dwi.nii,the output file becomes sub01_  Run).ACID retains the same folder structure and naming convention even when non-BIDS input is provided.

Pipelines
ACID is fully integrated into the SPM12 batch system, allowing users to execute its functions individually or combined into linear pipelines with multiple steps.Each step can receive the output of any of the previous steps via flexible and easy-to-use dependencies.While pipelines are typically set up in the SPM batch system, they can also be converted into MATLAB code (SPM batch script) for automation and further customization.In addition to its own functions, ACID integrates seamlessly with a range of standard SPM features, including segmentation, coregistration (based on affine transformation), spatial normalization (including nonlinear registration), and voxel-based statistical analyses, as well as a growing number of SPM extensions. 8For example, combining ACID with the hMRI toolbox enables multicontrast analysis of dMRI and other quantitative MRI data, such as relaxation rates ( Tabelow et al., 2019).

Example applications
To demonstrate the application of ACID toolbox on different types of dMRI data, here we provide three example pipelines for in vivo brain, in vivo spinal cord, and ex vivo dMRI (Fig. 3).Details of these three datasets are summarized in Table 4.The gradient schemes used for all datasets were based on the configurations proposed by Caruyer et al. (2013), available online. 9The design of the sampling schemes followed a uniform coverage on a sphere.Note that data with reverse phase-encoding direction were available for all three datasets, which refers to the acquisition of either a single b0 volume or all volumes with identical geometry and sequence parameters but opposite phase-encoding direction.4 for details on the datasets and Table 5 for details on the pipeline settings).Example batches for each type of dMRI data are stored in the Example_Batches folder of the toolbox.The positions of the displayed slices of the dMRI data are indicated in purple on the corresponding structural images.For the ex vivo specimen (C), the brain region from which the sample was extracted is highlighted in an orange box.Although not explicitly shown here, noise estimation should be performed on the unprocessed data (see Appendix C), which serves as input for msPOAS, Rician bias correction, and diffusion tensor fitting (for fitting methods WLS and robust fitting).However, in case of substantial misalignments across volumes, and when using the repeated measures noise estimation method, it might be beneficial to perform this step after ECMOCO to prevent an overestimation of noise.For msPOAS, a zoomed-in visual comparison is shown between a diffusion-weighted (DW) image before (middle row) and after applying msPOAS (bottom row); the msPOAS-corrected image appears less noisy while preserving tissue edges.For HySCO, contour lines of the corresponding structural image (displayed as red lines) are overlaid on a zoomed-in DW image both before (middle row) and after applying HySCO (bottom row).HySCO improves the alignment between the DW and the structural image.For the in vivo brain dMRI dataset (A), an inferior slice is shown that presents high susceptibility-related distortions, making the effect of HySCO more visible.For the ex vivo dMRI dataset (C), the effect of HySCO is shown in a slice (illustrated in yellow) orthogonal to the original one (illustrated in purple) to better visualize susceptibility-related distortions and their correction.Note that HySCO is applied as the final preprocessing step, that is, after applying msPOAS; however, the HySCO field map used for "unwarping" the diffusion-weighted images is estimated on the ECMOCO-corrected datasets, that is, before applying msPOAS.Rician bias correction (not explicitly shown here) should be applied either before (recommended: between msPOAS and HySCO, using the RBC module) or during model fitting (using the Rician bias correction option in NLLS).Diffusion signal models are fitted on the processed dataset; here, we display the maps of fractional anisotropy (FA) and mean kurtosis tensor (MW) from diffusion kurtosis imaging (DKI).The output from DKI can be used to compute biophysical parameters of the white matter; shown here is the map of Watson concentration parameter (κ) from the WMTI-Watson biophysical model.Note that for the in vivo brain dMRI dataset, the inferior slice displayed contains relatively little white matter; hence, we refrained from using a white matter mask.The less smooth appearance of the κ map is due to the low values in the gray matter.
et al., 2019) and is also available in ACID as an external tool, we refrained from including it in the example pipelines because the interaction between denoising and the interpolation associated with Gibbs ringing removal is not well characterized yet.We emphasize that these example pipelines might not be optimal for all cases; users might find that another combination of preprocessing steps, which might also include Gibbs ringing removal, works even better for their data.
While the pipelines for in vivo brain, in vivo spinal cord, and ex vivo dMRI follow similar concepts, recommended settings for each region may differ (Table 5).It is important to note that the settings listed in Table 5 serve as initial values for typical datasets.The optimal settings   In the "degrees of freedom" settings (ECMOCO), x, y, and z represent the frequency-, phase-, and slice-encoding directions, respectively.
for a particular dataset depend on the sequence parameters, the subject, and the imaged region.Model fitting may be followed by spatial processing, such as coregistration to the structural image or spatial normalization to a template in a standard space (e.g., MNI152 space), and statistical analysis (e.g., ROI-or voxel-based analysis).

DISCUSSION
We have developed the ACID toolbox, which extends the capabilities of the SPM framework by providing comprehensive preprocessing and model fitting techniques for in vivo brain, spinal cord, and ex vivo dMRI data.Besides commonly used diffusion signal models such as DTI and DKI, ACID also offers biophysical models that provide parameters of white matter tissue microstructure such as axonal water fraction and axon orientation dispersion.
Being seamlessly integrated into the SPM batch system, ACID allows for user-friendly access to SPM's powerful spatial processing tools and statistical framework.In addition to offering recommended pipelines for in vivo brain, spinal cord, and ex vivo dMRI, ACID provides the flexibility for users to create customized pipelines tailored to their specific data.Adhering to the BIDS conventions facilitates data sharing, enhances data comprehension for investigators, and makes ACID compliant with software requiring BIDS input (https://bids -apps .neuroimaging .io).

Preprocessing dMRI data
ACID offers artifact correction steps typically applied to dMRI data, including image realignment (ECMOCO), adaptive denoising (msPOAS), Rician bias correction (RBC), and correction for susceptibility-induced geometric distortions (HySCO).Here, we discuss specific considerations regarding their use for various applications.
Correcting for displacements within the dMRI data through image realignment is one of the most important but also challenging tasks.ECMOCO provides users with the flexibility to choose the degrees of freedom for image realignment based on the anticipated type of displacement, but also offers a selection of predefined degrees of freedom that are optimized for brain, spinal cord, and ex vivo dMRI.
In brain dMRI, motion can be approximated as a rigid body displacement with 6 degrees of freedom (DOF).Eddy-current spatial displacements, to a first-order approximation, result in translation and scaling along the phase-encoding direction (typically, the y-axis), and inplane and through-plane shearing ( Mohammadi et al., 2010).Since these displacements affect the entire brain, we recommend employing a 9-DOF volume-wise (volume to volume) registration with translation and rotation along x, y, and z, scaling along y, and shearing in the x-y and y-z planes.First-order approximation of eddy-current displacements might not always be sufficient, as dMRI data can also be affected by higher-order eddy-current field inhomogeneities causing nonlinear distortions ( J. L. R. Andersson & Sotiropoulos, 2016;Rohde et al., 2004).For example, in our observations, ECMOCO was not effective in removing pronounced eddy-current displacements present in the dMRI data of the Human Connectome Project ( Van Essen et al., 2012).In such cases, we recommend using FSL eddy, which incorporates higherorder eddy-current correction terms ( J. L. R. Andersson & Sotiropoulos, 2016) and can be called directly from ACID as an external tool (Section 2.6).In cases where ECMOCO is sufficient, an advantage of ECMOCO is that its performance is largely independent of the number of diffusion directions, whereas FSL eddy requires a minimum number of diffusion directions for good performance (see FSL website10 for recommendations).
In spinal cord dMRI, volume-wise registration has been found to be less effective ( Cohen-Adad et al., 2009;Mohammadi, Freund, et al., 2013) due to displacements that vary along the rostrocaudal axis of the spinal cord.These displacements appear mostly in the phaseencoding direction and are caused by physiological factors such as respiration and cardiac pulsation ( Kharbanda et al., 2006;Summers et al., 2006).We recommend applying volume-wise registration for rough alignment and correction of through-slice displacements, followed by slice-wise (slice to slice) registration for correcting any remaining slice-dependent displacement.This combined approach has demonstrated effectiveness in realigning not only volumes but also individual slices (Appendix Fig. B2), as well as improving the contrast-to-noise ratio between gray and white matter and reducing test-retest variability in DTI maps of the spinal cord ( Mohammadi, Freund, et al., 2013).Eddy-current distortions are typically less severe in the spinal cord compared with the brain, because the in-plane field of view is smaller and located near the scanner isocenter.This makes the firstorder approximation of eddy-current displacements, as supported by ECMOCO, generally adequate.We recommend employing a 4-DOF volume-wise registration (translation along x, y, z; scaling along y) followed by a 3-DOF slice-wise registration (translation along x, y; scaling along y).In datasets with low SNR, slice-wise correction along x can be omitted, given the smaller range of movement which makes reliable estimation difficult.We advise against correcting for in-plane rotation and shearing, as their expected range is very small.Correction for these DOFs might introduce spurious displacements during realignment, a risk we consider greater than not applying correction at all.
Structures surrounding the spinal cord (bones, ligaments, etc.) may move independently from the spinal cord, potentially leading to inaccuracies in transformation parameters.Moreover, as these structures typically occupy a larger portion of the image, they can dominate the estimation of transformation parameters.To address this challenge, ECMOCO provides the option of specifying a spinal cord mask to restrict the estimation of transformation parameters to the spinal cord and its immediate surroundings (Fig. 3).Any residual misalignments can be manually corrected using the Slice-wise realignment utility function (Table 2).
In ex vivo dMRI, specimen motion is not anticipated if the specimen is appropriately fixed, for instance, by using a sample holder or embedding it in agarose.Thus, we recommend correcting only for the four first-order eddy-current displacements (y-translation, y-scaling, x-y shearing, y-z shearing).The first-order approximation is typically adequate for small specimens where eddycurrent displacements are not severe.
In general, the performance of msPOAS and HySCO is largely independent of the anatomical features present in the image; therefore, default parameters are expected to work well for in vivo brain, spinal cord, and ex vivo dMRI data.Nevertheless, the default regularization parameters for HySCO (alpha "diffusion" and beta "Jacobian" regulator), accessible through the script config/local/acid_ local_defaults.m,are optimized for the brain and may require adjustment for the spinal cord if performance is inadequate.
Applying HySCO is particularly important for acquisitions with severe susceptibility-related distortions, such as multiband EPI without parallel imaging, and for multicontrast analyses where dMRI data or other quantitative maps are combined with structural reference images, for example, the dMRI-based axonal water fraction and magnetization transfer saturation maps in g-ratio mapping ( Mohammadi & Callaghan, 2021) or multicontrast MRI in the spinal cord ( David et al., 2019).In these cases, HySCO improves the overlap between the undistorted structural image and the dMRI data, improving the performance of subsequent coregistration and spatial normalization algorithms.HySCO has also been shown to improve the accuracy of g-ratio mapping ( Clark et al., 2021;Mohammadi, Tabelow, et al., 2015).While HySCO is far more efficient than FSL topup in terms of computation time ( Macdonald & Ruthotto, 2018), it does not integrate movement and susceptibility artifact correction into a single model.To mitigate the effects of subject movement, we propose acquiring images with reversed phase-encoding direction (the blip-up and blip-down images) in close succession.
The application of adaptive denoising (msPOAS) is important as it reduces the variance and, therefore, improves the precision of the tensor and kurtosis parameter estimates (see Supplementary Fig. S4 for an example illustrating the effect of msPOAS on DKI parameters, and refer to Becker et al. (2014) for more examples and details).For high-SNR data, denoising might not be advantageous; instead, denoising methods could even introduce additional error (see analysis in Appendix G).For low-SNR data, Rician bias correction (RBC), either applied to the dMRI data or during model fitting, must be performed in addition to msPOAS to mitigate the Rician bias in parameter estimates (see Appendix F for an example).An in-depth analysis of the impact of Rician bias correction on DKI and axisymmetric DKI can be found in Oeschger et al. (2023a).

Physical diffusion models
At a given b-value, the SNR in spinal cord dMRI is typically lower than in brain dMRI due to (i) the smaller crosssectional area that requires higher in-plane resolution (see Fig. 4A for a size comparison), (ii) the high signal attenuation for diffusion-gradient directions parallel to the highly aligned fibers in the head-feet direction (Fig. 4B), and (iii) the suboptimal coil configuration in the thoracic and lumbar regions, which are not covered by the head and neck coil.Lower SNR increases the variance of parameter estimates and makes spinal cord dMRI more susceptible to Rician bias.Consequently, SNR is often prohibitively low at higher b-values necessary for fitting the kurtosis tensor, making the application of DKI in the spinal cord very challenging.
The bias in parameters estimates induced by signal outliers from cardiac, respiratory, and other physiological artifacts ( Mohammadi, Hutton, et al., 2013) can be mitigated by applying robust fitting as a tensor fitting method (Appendix E.3).Given the higher occurrence of signal outliers in the spinal cord, robust fitting holds particular relevance for spinal cord dMRI.In a previous study, we demonstrated that robust fitting leads to higher FA values within the white matter and lower FA values within the gray matter in spinal cord dMRI data, resulting in an approximately 8% enhancement in contrast-tonoise ratio ( Mohammadi, Freund, et al., 2013).While robust fitting demonstrated high resistance to contamination (presence of outliers) compared with OLS and NLLS estimations, it is important to note that robust fitting requires a sufficiently large number of artifact-free data points.Simulations suggested that robust tensor estimates begin to break down when the frequencies of moderately intense cardiac pulsation artifacts exceed 27-30% ( Zwiers, 2010, fig. 5).
One potential limitation of linearized fitting methods is their operation on logarithmically transformed signals, where the assumption of Gaussian (or Rician) error distribution may not hold.The presence of logarithmically distorted Rician noise distribution not only restricts validity but can also impact the accuracy of the parameter estimates ( J. L. R. Andersson, 2008;Chang et al., 2005;Koay et al., 2006), particularly in the low-SNR regime such as in spinal cord dMRI.The WLS and robust fitting algorithms incorporate the signal intensity into the weights of the estimator function (Appendices E.2 and E.3), which was shown to reduce the effect of log-Rician distortion ( Salvador et al., 2005).Alternatively, the NLLS algorithm (Appendix E.4) can be used, which circumvents the distortion of the Rician distribution by operating on the original (nonlogarithmic) signals, and is, therefore, expected to yield more accurate parameter estimates, provided that the numerical fitting problem is sufficiently well conditioned.
In summary, for data with relatively high SNR and a frequent occurrence of outliers, we recommend using robust fitting to mitigate the influence of outliers.NLLS, particularly when combined with Rician bias correction, may be more suitable for dMRI data with lower SNR, which is often encountered in acquisitions for DKI (refer to Oeschger et al., 2023a, for recommended minimum SNR values and the Rician bias simulation utility function in Table 2 for simulating the Rician bias on dMRI data with a given SNR).Low-SNR data with a frequent occurrence of outliers pose challenges for model fitting, where a combination of msPOAS with RBC might reach their limits.In such cases, reliability masking can assist in identifying and excluding corrupted, thus unreliable, voxels from the parameter maps ( David et al., 2017).

Biophysical diffusion models
Of the biophysical models implemented in ACID, WMTI-Watson relies on DKI metrics (requiring at least two diffusion shells), while NODDI-DTI relies on DTI metrics (requiring a single diffusion shell only).This implies that the challenges associated with the estimation of DTI and DKI metrics, as discussed earlier, also apply to derived biophysical models.Accurate and precise estimation of DKI and DTI metrics is essential for the successful application of WMTI-Watson and NODDI-DTI, respectively.
In general, we recommend the DKI-based WMTI-Watson model over NODDI-DTI due to the fewer model assumptions, allowing it to better capture diffusion patterns in complex axonal configurations within the brain white matter.This aligns with the results from our example multishell brain dMRI dataset, where WMTI-Watson yielded more accurate estimates of κ and AWF compared with NODDI-DTI (Supplementary Fig. S5).For a more indepth comparison of biophysically derived values with histological values, refer to Papazoglou et al. (2024).
On the other hand, complex models are more "datahungry" and more susceptible to noise due to the higher number of fitted parameters, which can lead to poorly conditioned optimization problems when the amount and/or the quality of input data are insufficient.Therefore, for low-SNR data, as is often the case in spinal  4 for details on the datasets).(B) Schematic visualization of the spinal cord, highlighting the "butterfly-shaped" gray matter, which is located in the middle of the spinal cord and contains neuronal cell bodies and loosely aligned fibers, and the surrounding white matter, which contains highly aligned fibers.cord dMRI, the less complex but better-conditioned NODDI-DTI model might be the preferred choice.The low b-values often used in spinal cord dMRI could lead to inadequate parameter estimation when using the WMTI-Watson model.Indeed, NODDI-DTI yielded a more accurate estimation of κ in the example spinal cord dMRI dataset, whereas WMTI-Watson highly overestimated it (Supplementary Fig. S5).A drawback of the NODDI-DTI model in the spinal cord is its assumption of fixed intra-and extracellular diffusivities, which were optimized for the brain and might not be valid for the spinal cord.Both low SNR ( Veraart et al., 2011) and kurtosis bias ( Edwards et al., 2017) can lead to an underestimation of MD (Supplementary Fig. S3), impacting the model parameter estimation when MD falls outside the range where the NODDI-DTI model provides a valid representation (refer to Equation (4) in Edwards et al., 2017).This was evident in the estimation of AWF, which proved unfeasible in the spinal cord dataset (see Appendix Fig. F1; Supplementary Fig. S5).We anticipate that future improvements in acquisition methods will enhance the SNR in spinal cord dMRI, enabling the acquisition of higher b-values, which would alleviate many of the above-mentioned drawbacks.
A compromise between these two models could be the white matter tract integrity (WMTI) model, which is included as an external tool in ACID (Section 2.6).WMTI assumes highly aligned fibers, which holds true in white matter regions with high fiber alignment, such as the corpus callosum or the spinal cord white matter, but is less appropriate in regions with more complex axonal configurations, such as parts of the superior longitudinal fasciculus.
Ex vivo neuronal tissues exhibit different diffusivities compared with in vivo tissues due to various factors, including the effect of fixation, changes in chemical properties, and differences in temperature and composition of the embedding fluid.For example, white matter diffusivity was reported to reduce by approximately 85% from in vivo to ex vivo conditions, while the ratio between gray and white matter diffusivities remains similar, around 2-3 ( Roebroeck et al., 2019).To accommodate the reduced diffusivities under ex vivo conditions, ACID offers the option to utilize compartmental diffusivities tailored for ex vivo datasets within the NODDI-DTI model.Such an adjustment is not necessary for WMTI and WMTI-Watson, as their compartmental diffusivities are fitted rather than fixed.
We emphasize that NODDI-DTI, WMTI, and WMTI-Watson have been developed to characterize diffusion in the white matter.Recently, several efforts have been made to extend biophysical models to the gray matter ( Jelescu et al., 2020).Notable examples include the SANDI ( Palombo et al., 2020) and NEXI ( Jelescu et al., 2022) biophysical models.However, these models thus far, no study using these protocols on a clinical MRI system has been published.

Studies quantitatively evaluating the performance of ACID pipelines
Here, we briefly summarize and discuss the studies that quantitatively evaluated the performance of ACID tools individually or in comparison with other tools.

Evaluating preprocessing pipelines
In a previous study, we assessed the performance of ECMOCO as well as the combination of ECMOCO and msPOAS in simulated high-and low-SNR multishell brain dMRI datasets with added motion and eddy-current artifacts (i.e., perturbed data) (Mohammadi, Tabelow, et al., 2015).We found that the performance of ECMOCO in correcting the perturbed volumes was dependent on the SNR, with the number of incorrectly registered volumes increasing at lower SNR (SNR < 16).However, the combined application of msPOAS and ECMOCO effectively reduced the number of incorrectly registered volumes even at low SNR (fig. 3 in Mohammadi, Tabelow, et al., 2015).Additionally, correcting the perturbed volumes with ECMOCO and msPOAS yielded FA maps closer to the "ground truth," that is, the FA map computed on the unperturbed data (fig.5 in Mohammadi, Tabelow, et al., 2015).In another study utilizing clinical spinal cord dMRI data, we evaluated the impact of ECMOCO on the group differences observed in FA between patients with degenerative cervical myelopathy and healthy controls (fig.7 in David et al., 2017).Our analysis revealed that ECMOCO had only a minimal effect on the two-sample t-score computed between the FA values of the two groups.
We also tested the effects of different denoising methods (msPOAS, LPCA, and MP-PCA) on the accuracy of DKI metrics, with the details and results described in Appendix G.In short, we found that denoising (using any of the three methods) is beneficial only in the low-SNR domain (below an SNR of approximately 30).In high-SNR data, denoising did not lead to further improvements with MP-PCA and even introduced additional errors with msPOAS and LPCA.In terms of susceptibility artifacts, we previously found in a brain dMRI dataset that FSL topup was more efficient in correcting susceptibilityrelated distortions than HySCO, even when including a motion correction step between the reverse phaseencoded (blip-up and blip-down) images before running HySCO (fig. 3 in Clark et al., 2021).This is potentially because the HySCO pipeline involved multiple interpolation steps, introducing additional blurring effects, while FSL topup incorporates motion and susceptibility distortion correction within the same model.The same study found that combining reverse phase-encoded images using the "weighted average" method (HySCO: combine blip-up and blip-down images module), as opposed to the "arithmetic average" method, reduces image blurring in the corrected brain dMRI data and achieves greater overlap between the dMRI data and the corresponding structural image.In fact, when using the "weighted average" method, HySCO performed comparably to FSL topup and even outperformed it in regions suffering from high levels of distortion (fig.5 in Clark et al., 2021).In spinal cord dMRI, a previous study found that HySCO is comparable with other distortion correction tools such as FSL topup ( Schilling et al., 2024).

Evaluating diffusion signal models
In brain dMRI datasets, we found that robust tensor fitting can reduce the effect of signal outliers due to motion, eddy-current artifacts, incorrectly registered volumes (fig.5C, D in Mohammadi, Tabelow, et al., 2015), or physiological noise (fig.9 in Mohammadi, Hutton, et al., 2013).In spinal cord dMRI, we quantified the performance of robust fitting and showed that it can reduce the bias in FA, especially at tissue boundaries (fig.7 in Mohammadi, Freund, et al., 2013).On the other hand, robust fitting had only a minor effect on group differences in FA between patients with degenerative cervical myelopathy and healthy controls, regardless of whether using the ACID implementation of robust fitting or using RESTORE (part of the CAMINO toolbox; Chang et al., 2012) (fig.7 in David et al., 2017).However, within the same study, we also found that supplementing the pipeline with reliability masking to exclude outlier voxels (Table 2) considerably increased the statistical differences between patients and controls (fig.7 in David et al., 2017).

Applications
For all applications, it is highly recommended to assess the data quality before and after each processing step.In addition to the quality assessment utility functions DWI series browser and DWI series movie (Table 2), multiple ACID modules generate diagnostic plots to identify the presence and type of artifacts in the dMRI data.Example diagnostic plots are provided in Supplementary Figures S1 and S2.

Integration with SPM modules
ACID can be readily combined with SPM tools for segmentation, spatial processing, and voxel-based analysis of parametric maps.Segmenting the brain or spinal cord is often necessary for coregistration, spatial normalization, or tissue-specific analyses.In the brain, tissue probability maps of white matter, gray matter, and cerebrospinal fluid can be created by unified segmentation, the default segmentation routine in SPM12 ( Ashburner & Friston, 2005).These tissue probability maps can also be used to create a binary brain mask using the Create brain mask utility function (Table 2).To enable SPM's unified segmentation in the spinal cord, the brain tissue priors need to be substituted with the joint brain and spinal cord tissue priors from the probabilistic brain and spinal cord atlas ( Blaiotta et al., 2017).However, this atlas only covers the upper cervical cord down to C3; for other spinal levels, the user is referred to automatic (e.g., deepseg ( Perone et al., 2018)) or semiautomatic (e.g., active surface method ( Horsfield et al., 2010)) segmentation techniques.
Brain dMRI data can be coregistered to the corresponding structural image using spm_coreg.For nonlinear spatial registration to the MNI space, we recommend SPM DAR-TEL ( Ashburner, 2007) or Geodesic Shooting ( Ashburner & Friston, 2011).As SPM registration tools often rely on brain tissue priors, they cannot be applied directly on spinal cord dMRI.For the spinal cord, we recommend utilizing the PAM50 template ( De Leener et al., 2018) and the corresponding normalization tools integrated into the Spinal Cord Toolbox ( De Leener et al., 2017).
While brain and spinal cord images are typically analyzed separately, there are scenarios where combining them into a single image can be beneficial.For example, when registering the brain and spinal cord image to a joint brain-spinal cord template, such as the probabilistic atlas of the brain and spinal cord ( Blaiotta et al., 2017), the warping field is often obtained using a structural image with a large field of view (FOV) covering both regions (Fig. 5).To apply this warping field to the brain and spinal cord images, they need to be fused into a single image.ACID provides the Fusion utility function (Table 2) which merges two distinct images, acquired with different FOV and geometric properties, into a unified large FOV image (Fig. 5).
ACID benefits from SPM's rich statistical framework for voxel-based analysis.SPM's second-level analysis tool (SPM -> Specify 2nd-level) performs voxelbased statistical tests on the parametric maps using t-test, ANOVA, or general linear model.In the SPM -> Results module, the framework also offers (i) multiple comparison correction in the form of family-wise error rate and false discovery rate, (ii) thresholding the test statistics at cluster level and voxel level and providing a list of significant clusters/voxels, and (iii) various visualization tools for displaying and saving the significant clusters.Furthermore, ACID's ROI analysis utility function (Table 2) can be used to extract mean metrics within subject-specific ROIs in the native space or perform atlas-based analysis in the template space.For atlas-based analysis in the spinal cord, the user is referred to the PAM50 white and gray matter atlas ( De Leener et al., 2018).
Although ACID does not provide tractography or tractbased analysis tools, the output of its model fitting methods can be input into tractography tools such as FSL or the SPM12-based Fibertools toolbox (see Wiki11 on the git repository for more details).

Computation time
To speed up the processing and analysis of dMRI data, parallel computing is implemented wherever applicable.This technique can substantially accelerate the most time-consuming ACID modules, including ECMOCO and DTI/DKI fit.Note that parallel computing requires the Parallel Computing Toolbox in MATLAB.Table 6 provides the computation times for selected ACID modules on a typical brain and spinal cord dMRI dataset.

Research applications
ACID has been used in a variety of clinical and neuroscience research, for example, in dMRI studies assessing cerebral changes in patients with multiple sclerosis ( Deppe, Krämer, et al., 2016;Deppe, Tabelow, et al., 2016;Dossi et al., 2018;Kugler & Deppe, 2018) and Parkinson's disease ( Szturm et al., 2021), and to assess gliomas ( Paschoal et al., 2022;Raja et al., 2016).We have also used ACID to investigate spinal cord white matter following spinal cord injury ( Büeler et al., 2024;David et al., 2019David et al., , 2021David et al., , 2022;;Grabher et al., 2016;Huber et al., 2018;Seif et al., 2020;Vallotton et al., 2021).A noncomprehensive list of studies using the ACID toolbox can be found on the project website.12Note that certain ACID functions can be applied to MRI data beyond dMRI as well; for instance, HySCO has been used to correct brain fMRI data for susceptibility artifacts ( De Groote et al., 2020).It is important to note that ACID has not been approved for clinical applications by any health agency and it comes with no warranty.Therefore, it should not be used for diagnosis in clinical settings.2).The two images should ideally share an overlapping region, but they may have different geometric properties such as resolution and number of slices.In the overlapping region, the voxel intensity values are computed as the average of the intensity values from the two images.The merging process requires a structural image as the registration target.The combined FA map is resampled onto the higher resolution structural image, resulting in a smoother appearance.Comparing the tools within the ACID toolbox with alternative implementations in other software presents challenges because their performance depends on the specific dMRI data and the chosen parameter settings from a potentially large parameter space, which necessitates a systematic exploration of the parameter space.In addition, the evaluation of entire processing pipelines would drastically increase the number of parameters to test.While we have outlined the comparisons conducted so far in Section 4.3, we assert that a thorough quantitative comparison between toolboxes warrants a dedicated future study.In general, we encourage users to undertake such comparisons on their own datasets.
The ACID toolbox is the result of a collaborative effort to extend the SPM ecosystem with state-of-the-art processing and modeling tools for dMRI data.Our aim is to make the toolbox widely accessible, leveraging SPM's large and vibrant community.Users can submit their questions, bug reports, and suggestions via the dedicated mailing list or by opening an issue on the git website.This paper offers an overview of the current state of the toolbox, with several ongoing developments not covered here.The modularity of the toolbox allows for integration of newly developed methods, even when used concurrently with old ones.Biophysical modeling is an emerging field, and we expect many methodological advancements to occur in the coming years.To align with this ongoing development, our goal is to consistently integrate state-of-the art biophysical models into ACID.We also plan to add the Rician maximum likelihood estimator ( Sijbers et al., 1998) as an alternative to the existing quasi-likelihood estimators ( Polzehl & Tabelow, 2016).

CONCLUSION
ACID is an open-source extension to SPM12 that provides a comprehensive framework for processing and analyzing in vivo brain, spinal cord, and ex vivo dMRI data.The toolbox was developed to meet the increasing demand for studies involving spinal cord dMRI, research employing biophysical models, and validation studies utilizing ex vivo dMRI.ACID leverages the core SPM tools and other SPM extensions, which can be easily integrated into the ACID pipeline.

DATA AND CODE AVAILABILITY
The source code of ACID is freely available at https:// bitbucket .org / siawoosh / acid -artefact -correction -in -diffusion -mri / src / master/.The authors will make the raw data used for the visualizations in this article available in an associate publication.( Tabelow et al., 2015).It is important to note that msPOAS assumes a single global value of sigma, which may not always hold.If sigma is correctly estimated, the default lambda value will ensure optimal adaptation.Incorrect estimation of sigma can be compensated by the choice of lambda, which makes msPOAS robust against misspecification of sigma ( Becker et al., 2014).We recommend determining kappa automatically based on the number of diffusion directions ( Tabelow et al., 2015).However, manual adjustment of kappa may be necessary in cases where the SNR is low (e.g., for spinal cord dMRI) or if the dataset has more images with high b-values than with low b-values.The effective number of coils (ncoils) is 1 when using SENSE1 reconstructions ( Polzehl & Tabelow, 2016;Sotiropoulos et al., 2013), but the correct value is more difficult to determine when using multiple receiver channels ( Aja-Fernández et al., 2014).It is important to use the same ncoils for the estimation of sigma and in msPOAS to ensure the same number of degrees of freedom.
where Ŝi represents the unknown true signal (without noise) and σ i is the background noise for acquisition i.The estimator function now becomes ρ(ε i ) = ( ω i ε i ) 2 , yielding the solution α α wls = W T B T WB ( ) W T B T Wy, with W being the diagonal matrix of ω i .Note that OLS is a special case of WLS, where ω i = 1 for all i.A practical consideration in obtaining α α wls is related to estimating Ŝi .One approach is to use the measured noisy signal S i as an estimate of Ŝi .Another approach is to start with the OLS solution and use the fitted signal as an estimate of Ŝi , which was shown to be more accurate ( Veraart, Sijbers, et al., 2013).

Appendix E.3. Robust fitting
The concept behind robust fitting is to assign lower weights to data points with higher model-fit errors during the fitting process ( Mangin et al., 2002).The robust fitting method implemented in ACID is based on the "Patching ArTefacts from Cardiac and Head motion" (PATCH) technique introduced by Zwiers (2010).While the form of the estimator function is similar to that of WLS, PATCH factorizes the weighting function ω i into three components as ω i = ω i1 ω i2 ω i3 , where each component is designed to address different types of artifacts: ω i1 and ω i2 account for regional and slice-wise artifacts, respectively, while ω i3 is identical to the weight term in WLS.ω i1 and ω i2 are exponentially decaying functions of ε i :  4 for details on the datasets).Being derived from white matter biophysical models, the parameter maps were masked for the white matter in the brain dataset.For the spinal cord, we refrained from masking due to the difficulty of obtaining an accurate white matter mask.These maps were computed both without (left column) and with (middle column) RBC; their voxel-wise difference, referred to as the Rician bias, is shown in the right column.RBC slightly decreased the mean of the kurtosis tensor in both the brain and spinal cord, which resulted in an increase in κ.The estimation of AWF using the NODDI-DTI model was not feasible in the spinal cord, as the mean diffusivity (MD) values derived from DTI fell below the range where the NODDI-DTI model provides a valid representation (refer to Equation (4) in Edwards et al., 2017).This discrepancy could be attributed to either the underestimation of MD due to kurtosis bias (Supplementary Fig. S3) or the invalidity of fixed compartmental diffusivities in the NODDI-DTI model.
as C 1 = 1.4826•median ε i ( ) and C 2 = 1.4826•median ε i,sl ( ) ( Hampel, 1974;Rousseeuw & Croux, 1993).Note that accurate estimation of C 1 and C 2 is crucial for effectively down-weighting outliers.This holds true as long as outliers are sparsely distributed and the median of the model-fit errors remains unaffected.However, a frequent occurrence of outliers can increase C, leading to a less effective down-weighting of outliers.While OLS and WLS independently fit the tensor in each voxel, PATCH makes use of the observation that physiological noise represents a structured, spatially correlated noise.To accommodate the anticipated smoothness of C 1 , the median operator is spatially smoothed using a 2D Gaussian kernel before computing C 1 ( Zwiers, 2010).
As a modification to PATCH, the robust fitting met hod incorporates Tikhonov regularization to handle ill-conditioned weighting matrices resulting from a high occurrence of outliers.This leads to the solution α α λ λ = W T B T WB + λ B T B ⎡ ⎣ ⎤ ⎦ −1 W T B T Wy , where W represents the diagonal matrix of factorized weights, and λ is the Tikhonov regularization factor.Notably, in the two extreme cases, the Tikhonov solution either becomes α wls (albeit with a different W) (λ = 0) or converges to α α ols (λ → ∞).The above equation cannot be solved readily, as W is a function of ε, which is only available after obtaining the solution.This is addressed by using an iteratively reweighted least squares (IRLS) algorithm.In the first iteration, ω i is set to 1 for all i to obtain the OLS solution α α ols and the initial ε ε.In the second iteration, an updated W is computed based on the initial ε, which is then used to compute α α λ λ .In each further iteration, ε ε from the preceding iteration is used to update W, which is in turn used to compute the updated α α λ λ .This iterative process is repeated until convergence or until the predefined number of iterations is exceeded.
2.2.Example maps are shown in Figure 2, and specific values obtained from the brain and spinal cord are presented and discussed in Supplementary Figure S5 (Supplementary Material).

Fig. 2 .
Fig. 2. Maps of biophysical parameters derived from the WMTI-Watson model using an in vivo brain, in vivo spinal cord, and ex vivo dMRI dataset (refer to Table4for details on the dataset).Shown are maps of Watson concentration parameter (κ), axonal water fraction (f ), parallel and perpendicular extra-axonal diffusivities (D e,! and D e,⊥ ), and intra-axonal diffusivity (D a ).Note that for the in vivo spinal cord dataset, the maximum b-value (b = 1500 s/mm 2 ) was probably too low for an accurate estimation of D e,! , resulting in voxels with negative (hence unphysical) values within the spinal cord.Since WMTI-Watson is a white matter biophysical model, the parameter maps were masked for the white matter in the brain dataset.For the spinal cord and ex vivo specimen, we refrained from masking for the white matter due to the difficulty of obtaining an accurate white matter mask.

Fig. 3 .
Fig. 3. Standard processing pipelines for typical (A) in vivo brain, (B) in vivo spinal cord, and (C) ex vivo dMRI datasets (refer to Table4for details on the datasets and Table5for details on the pipeline settings).Example batches for each type of dMRI data are stored in the Example_Batches folder of the toolbox.The positions of the displayed slices of the dMRI data are indicated in purple on the corresponding structural images.For the ex vivo specimen (C), the brain region from which the sample was extracted is highlighted in an orange box.Although not explicitly shown here, noise estimation should be performed on the unprocessed data (see Appendix C), which serves as input for msPOAS, Rician bias correction, and diffusion tensor fitting (for fitting methods WLS and robust fitting).However, in case of substantial misalignments across volumes, and when using the repeated measures noise estimation method, it might be beneficial to perform this step after ECMOCO to prevent an overestimation of noise.For msPOAS, a zoomed-in visual comparison is shown between a diffusion-weighted (DW) image before (middle row) and after applying msPOAS (bottom row); the msPOAS-corrected image appears less noisy while preserving tissue edges.For HySCO, contour lines of the corresponding structural image (displayed as red lines) are overlaid on a zoomed-in DW image both before (middle row) and after applying HySCO (bottom row).HySCO improves the alignment between the DW and the structural image.For the in vivo brain dMRI dataset (A), an inferior slice is shown that presents high susceptibility-related distortions, making the effect of HySCO more visible.For the ex vivo dMRI dataset (C), the effect of HySCO is shown in a slice (illustrated in yellow) orthogonal to the original one (illustrated in purple) to better visualize susceptibility-related distortions and their correction.Note that HySCO is applied as the final preprocessing step, that is, after applying msPOAS; however, the HySCO field map used for "unwarping" the diffusion-weighted images is estimated on the ECMOCO-corrected datasets, that is, before applying msPOAS.Rician bias correction (not explicitly shown here) should be applied either before (recommended: between msPOAS and HySCO, using the RBC module) or during model fitting (using the Rician bias correction option in NLLS).Diffusion signal models are fitted on the processed dataset; here, we display the maps of fractional anisotropy (FA) and mean kurtosis tensor (MW) from diffusion kurtosis imaging (DKI).The output from DKI can be used to compute biophysical parameters of the white matter; shown here is the map of Watson concentration parameter (κ) from the WMTI-Watson biophysical model.Note that for the in vivo brain dMRI dataset, the inferior slice displayed contains relatively little white matter; hence, we refrained from using a white matter mask.The less smooth appearance of the κ map is due to the low values in the gray matter.

Fig. 4 .
Fig. 4. (A) Illustration of differences in the cross-sectional area between the brain and spinal cord, displaying a single axial slice of the mean T2-weighted (b0) image (refer to Table4for details on the datasets).(B) Schematic visualization of the spinal cord, highlighting the "butterfly-shaped" gray matter, which is located in the middle of the spinal cord and contains neuronal cell bodies and loosely aligned fibers, and the surrounding white matter, which contains highly aligned fibers.

Fig. 5 .
Fig.5.Merging of two fractional anisotropy (FA) maps, covering the brain and cervical cord, respectively, into a unified FA map using the Fusion utility function (Table2).The two images should ideally share an overlapping region, but they may have different geometric properties such as resolution and number of slices.In the overlapping region, the voxel intensity values are computed as the average of the intensity values from the two images.The merging process requires a structural image as the registration target.The combined FA map is resampled onto the higher resolution structural image, resulting in a smoother appearance.
model-fit error, with n being the number of voxels within the slice.A 1 and A 2 are model parameters, by default set to 0.3 and 0.1, respectively, with higher values resulting in a faster exponential decay.C 1 and C 2 are estimates of the standard deviation of ε i and ε i,sl , respectively, in the absence of outliers, and are computed Appendix Fig.F1.The impact of Rician bias correction (RBC) on maps of biophysical parameter estimates, derived from the NODDI-DTI and WMTI-Watson models, including Watson concentration parameter (κ) and axonal water fraction (AWF), in an in vivo brain and spinal cord dataset (refer to Table

Table 1 .
Peer-reviewed methods used in the ACID toolbox.
DKI, diffusion kurtosis imaging; DTI, diffusion tensor imaging; NLLS, nonlinear least squares; NODDI, neurite orientation dispersion and density imaging; WMTI, white matter tract integrity.*The ACID implementation is based on the method introduced by Jespersen et al. (2018).

Table 3 .
List of labels in the output filename's desc field (not comprehensive).When multiple processing steps are involved, the labels are concatenated, as in sub01_desc-ECMOCO-msPOAS_dwi.nii.Model fitting appends three labels indicating the type of diffusion model, algorithm, and parametric map, such as sub01_ desc-ECMOCO-msPOAS-DKI-OLS-FA_dwi.nii.For BIDS-compliant input, ACID generates a bval and bvec file after each processing step.ACID stores all output in the derivatives folder, with separate subfolders for each module's output (e.g., derivatives/msPOAS-

Table 4 .
Scan parameters of the in vivo brain, in vivo spinal cord, and ex vivo dMRI datasets used in this paper.

Table 5 .
Settings of selected modules for in vivo brain, in vivo spinal cord, and ex vivo dMRI datasets.

Table 6 .
Computation times of selected ACID modules on an example in vivo brain and in vivo spinal cord dMRI dataset (refer to Table4for details on the datasets), when run on a MacBook M1 laptop (4 cores, 16GB RAM).