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

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).

To fully utilize the sensitivity of dMRI to tissue microstructure, expert knowledge is required to minimize artifacts both during acquisition, for example, by cardiac gating or twice-refocused spin-echo sequences, and through dedicated retrospective correction methods. Commonly used retrospective correction techniques include motion and eddy-current correction (J. L. R. Andersson & Sotiropoulos, 2016; Mohammadi et al., 2010), susceptibility distortion correction (Gu & Eklund, 2019; Ruthotto et al., 2012), denoising (Becker et al., 2014; Veraart et al., 2016), Rician bias correction (Oeschger et al., 2023a; Sijbers et al., 1998), and robust tensor fitting techniques (Chang et al., 2005; Mohammadi, Freund, et al., 2013). Retrospective artifact correction techniques, along with diffusion signal modeling capabilities, are widely available in open-source toolboxes such as FSL-FDT (Smith et al., 2004), DiPY (Garyfallidis et al., 2014), DESIGNER (Ades-Aron et al., 2018), ExploreDTI (Leemans et al., 2009), MRtrix3 (Tournier et al., 2019), TORTOISE (Pierpaoli et al., 2010), AFNI-FATCAT (Taylor & Saad, 2013), and others.

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) state-of-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.

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-of-the-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.

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.

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 signal-to-noise ratio (SNR), particularly in volumes with high b-values (>1000 s/mm²) 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.

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, 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.

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.

Hyperelastic Susceptibility Artifact Correction (HySCO) is a technique used to correct for geometric distortions caused by susceptibility artifacts (Ruthotto et al., 2012, 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 be combined using the submodule HySCO: combine blip-up and blip-down images.

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).

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 ill-conditioned 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 (Mohammadi, Freund, et al., 2013) (Appendix Fig. E1). Unlike the linearized models, the nonlinear least squared (NLLS) method is based on an implementation (Modersitzki, 2009) of the Gauss-Newton algorithm and operates on the nonlogarithmic data, avoiding the distortion of the noise distribution.

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²). 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 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.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).

The white matter tract integrity (WMTI)-Watson model as an implementation of the Standard Model assumes two nonexchanging water compartments (intra- and extra-axonal tissue water) (Alexander et al., 2019; Novikov et al., 2019). The model characterizes the intra-axonal compartment as “sticks” of zero radius, with an intra-axonal diffusivity $Da$ and axonal water fraction $f$⁠. Axonal alignment is characterized by the Watson concentration parameter $κ$⁠, where higher values indicate higher axonal alignment, and the orientation dispersion index (ODI), where higher values indicate lower alignment. While $κ$ and ODI are mathematically related (Mollink et al., 2017), ACID outputs both for convenience. The extra-axonal space is modeled as a homogeneous medium, described by an axisymmetric diffusion tensor with parallel (⁠$De,∥$⁠) and perpendicular (⁠$De,⊥$⁠) extra-axonal diffusivities. The five biophysical parameters (⁠$Da$⁠, $f$⁠, $κ$⁠, $De,∥$⁠, $De,⊥$⁠) are derived from the axisymmetric DKI tensor metrics (AD, RD, MW, AW, RW) according to the formulas described in Jespersen et al. (2018) and Novikov et al. (2018). Being derived from the biophysical Standard Model, the estimation of WMTI-Watson biophysical parameters is generally degenerate, which leads to two solutions: the plus branch, which assumes $Da$ > $De,∥$⁠, and the minus branch, which assumes $Da$ < $De,∥$ (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 (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 assumes a fixed diffusivity of the intraneurite compartment (⁠$Da$⁠). In our implementation, users can either choose from two fixed values tailored for in vivo (⁠$Da$ = 1.7∙10^–3 mm²/s) and ex vivo (⁠$Da$ = 0.6∙10^–3 mm²/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.

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).

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² (J. L. R. Andersson & Sotiropoulos, 2016); (ii) FSL topup³ (Smith et al., 2004); (iii) dwidenoise⁴ (based on the Marchenko-Pastur principal component analysis (MP-PCA), part of the MRtrix toolbox) (Veraart et al., 2016); (iv) denoising⁵ (based on the local principal component analysis (LPCA), part of the DWI Denoising Software) (Manjón et al., 2013); (v) Koay’s noise estimation⁶; (vi) mrdegibbs⁷ for Gibbs ringing removal, part of the MRtrix toolbox (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.

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 Pre-processing, Diffusion tensor/kurtosis imaging, and Biophysical models). For instance, after applying ECMOCO to sub01_dwi.nii, the output file becomes sub01_desc-ECMOCO_dwi.nii. 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-Run). ACID retains the same folder structure and naming convention even when non-BIDS input is provided.

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.⁸ For 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).

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.⁹ The 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. All example pipelines consist of artifact correction (ECMOCO, msPOAS, RBC, HySCO) and model fitting steps. While Gibbs ringing removal is often part of dMRI processing pipelines (Ades-Aron et al., 2018; Kellner et al., 2016; Tournier 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 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).

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).

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 in-plane 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 higher-order 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 website¹⁰ 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 phase-encoding 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 first-order 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 eddy-current 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).

At a given b-value, the SNR in spinal cord dMRI is typically lower than in brain dMRI due to (i) the smaller cross-sectional 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-to-noise 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).

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 in-depth comparison of biophysically derived values with histological values, refer to Papazoglou et al. (2024).

On the other hand, complex models are more “data-hungry” 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 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.

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

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 susceptibility-related distortions than HySCO, even when including a motion correction step between the reverse phase-encoded (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).

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).

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.

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 DARTEL (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 voxel-based 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 tract-based 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 Wiki¹¹ on the git repository for more details).

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.

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., 2019, 2021, 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.¹² Note 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.

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).

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.

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.

Gergely David: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Resources, Software, Visualization, Writing—original draft, Writing—review & editing. Björn Fricke: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing—original draft, Writing—review & editing. Jan Malte Oeschger: Formal analysis, Methodology, Software, Writing—original draft, Writing—review & editing. Lars Ruthotto: Methodology, Software, Writing—review & editing. Francisco J. Fritz: Data curation, Resources. Ora Ohana: Data curation, Resources. Laurin Mordhorst: Software. Thomas Sauvigny: Data curation, Resources. Patrick Freund: Conceptualization, Project administration, Writing—review & editing. Karsten Tabelow: Conceptualization, Investigation, Methodology, Project administration, Software, Writing—review & editing. Siawoosh Mohammadi: Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Writing—original draft, Writing—review & editing.

The authors declare no competing interests.

Three dMRI datasets from previous studies were reused in this paper. These studies complied with the principles of the Declaration of Helsinki and were approved by the local ethics committee (Ärztekammer Hamburg). The whole-brain dataset was measured in vivo (ethics approval number: PV5600). The dataset of the temporal lobe specimen was acquired ex vivo (PV5034). The spinal cord dataset was measured in vivo (PV5141).

This work was supported by the German Research Foundation (DFG Priority Program 2041 “Computational Connectomics” (MO 2397/5-1, MO 2397/5-2)), the Emmy Noether Stipend (MO 2397/4-1 and 2397/4-2), the BMBF (01EW1711A and B) in the framework of ERA-NET NEURON, and the ERC (Acronym: MRStain, Grant agreement ID: 101089218, DOI: 10.3030/101089218). Views and opinions expressed are, however, those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. Lars Ruthotto is supported in part by NSF awards (DMS 1751636 and DMS 2038118). Patrick Freund is funded by an SNF Eccellenza Professorial Fellowship grant (PCEFP3_181362/1).

Supplementary material for this article is available with the online version here: https://doi.org/10.1162/imag_a_00288