{'fraction_non_alphanumeric': 0.038021298001176894, 'fraction_numerical': 0.01833678176192555, 'mean_word_length': 4.68, 'pattern_counts': {'":': 1, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Neural networks for parameter estimation in microstructural MRI: Application to a diffusion-relaxation model of white matter

Specific features of white matter microstructure can be investigated by using biophysical models to interpret relaxation-diffusion MRI brain data. Although more intricate models have the potential to reveal more details of the tissue, they also incur time-consuming parameter estimation that may converge to inaccurate solutions due to a prevalence of local minima in a degenerate fitting landscape. Machine-learning fitting algorithms have been proposed to accelerate the parameter estimation and increase the robustness of the attained estimates. So far, learning-based fitting approaches have been restricted to microstructural models with a reduced number of independent model parameters where dense sets of training data are easy to generate. Moreover, the degree to which machine learning can alleviate the degeneracy problem is poorly understood. For conventional least-squares solvers, it has been shown that degeneracy can be avoided by acquisition with optimized relaxation-diffusion-correlation protocols that include tensor-valued diffusion encoding. Whether machine-learning techniques can offset these acquisition requirements remains to be tested. In this work, we employ artificial neural networks to vastly accelerate the parameter estimation for a recently introduced relaxation-diffusion model of white matter microstructure. We also develop strategies for assessing the accuracy and sensitivity of function fitting networks and use those strategies to explore the impact of the acquisition protocol. The developed learning-based fitting pipelines were tested on relaxation-diffusion data acquired with optimal and sub-optimal acquisition protocols. Networks trained with an optimized protocol were observed to provide accurate parameter estimates within short computational times. Comparing neural networks and least-squares solvers, we found the performance of the former to be less affected by sub-optimal protocols; however, model fitting networks were still susceptible to degeneracy issues and their use could not fully replace a careful design of the acquisition protocol.

far, learning-based fitting approaches have been restricted to microstructural models with a reduced number of independent model parameters where dense sets of training data are easy to generate. Moreover, the degree to which machine learning can alleviate the degeneracy problem is poorly understood. For conventional least-squares solvers, it has been shown that degeneracy can be avoided by acquisition with optimized relaxation-diffusion-correlation protocols that include tensor-valued diffusion encoding. Whether machine-learning techniques can offset these acquisition requirements remains to be tested. In this work, we employ artificial neural networks to vastly accelerate the parameter estimation for a recently introduced relaxation-diffusion model of white matter microstructure. We also develop strategies for assessing the accuracy and sensitivity of function fitting networks and use those strategies to explore the impact of the acquisition protocol. The developed learning-based fitting pipelines were tested on relaxationdiffusion data acquired with optimal and sub-optimal acquisition protocols. Networks trained with an optimized protocol were observed to provide accurate parameter estimates within short computational times. Comparing neural networks and least-squares solvers, we found the performance of the former to be less affected by sub-optimal protocols; however, model fitting networks were still susceptible to degeneracy issues and their use could not fully replace a careful design of the acquisition protocol.

Introduction
Microstructure imaging uses compartment modelling of diffusion MRI (dMRI) data with the aim to map specific tissue quantities (Alexander et al., 2019;Nilsson et al., 2013;Novikov et al., 2019). A central goal in microstructure imaging has been to estimate the volume fractions of different microstructural components such as axons (Lampinen et al., 2020(Lampinen et al., , 2019Veraart et al., 2018). Estimating volume fractions rather than signal fractions is challenging, however, because it requires the simultaneous estimation of both diffusion and relaxation properties of the different model compartments. This kind of inverse problem is sensitive to degeneracy issues (Jelescu et al., 2016;Lampinen et al., 2019), in which multiple sets of model parameters can describe the acquired data equally well. Parameter estimation can also be computationally slow, preventing real-time mapping. A potential solution is to employ machine learning to accelerate the parameter estimation process (Golkov et al., 2016). However, the current literature lacks systematic assessments of the advantages and drawbacks of this approach, which is surprising considering the exponential increase in interest for such methods.
Artificial neural networks (ANNs) and other machine learning approaches have been applied previously to accelerate the estimation of microstructure parameters from dMRI data (Barbieri et al., 2020;Bertleff et al., 2017;Golkov et al., 2016;Grussu et al., 2020;Gyori et al., 2019;Hill et al., 2021;Kaandorp et al., 2021;Nedjati-Gilani et al., 2017;Palombo et al., 2020;Reisert et al., 2017). For example, a random forest regressor has been used to fit a compartment model for white matter (WM) microstructure in the presence of water exchange (Nedjati-Gilani et al., 2017) and to fit the SANDI model for grey matter properties . Reisert et al. (2017) applied machine learning to a Bayesian estimation approach which dramatically accelerated the fitting of two-and threecompartment models. Barbieri et al. (2020) applied ANNs to the intra-voxel incoherent motion model. An important open question, however, is what impact the training strategy has on the fitting performance. This is particularly relevant when applied to non-linear multi-compartment models with many independent model parameters, which we here refer to loosely as 'high-dimensional models'. The generation of training data scales poorly with the number of model parameters, as sampling each combination of p model parameters in m steps requires m p samples. As p increases, it is unavoidable that a finite set of samples becomes sparse in the p-dimensional model parameter space, risking selection bias.
Here, we investigate the impact of different sampling patterns within this space on the performance of the neural network. Apart from accelerating model fitting, neural networks may in principle also reduce the requirements on the imaging protocol by learning priors from training examples (Golkov et al., 2016). For example, neural networks have been used to learn a mapping between fully-sampled and sub-sampled datasets, which can in turn be used to stabilise model fitting performance against substantial degrees of data down-sampling Tian et al., 2020). However, we do not expect machine learning approaches to completely alleviate degeneracy issues. Indeed, for cases where the acquisition protocol does not provide sufficient information to resolve between different parameter values, the learning-based estimates will simply equal the mean of the model parameter distribution used for training (Reisert et al., 2017).
The aims of this study were to compare training strategies, to pro-pose tools to evaluate the performance of model fitting neural networks, and to test to what degree neural networks can solve problems with degeneracy. As a testbed, we use a high-dimensional relaxationdiffusion microstructure model of WM (Lampinen et al., 2020(Lampinen et al., , 2019Veraart et al., 2018). For this model, parameter estimation is enabled by state-of-the-art imaging protocols featuring so-called b-tensor encoding (Topgaard, 2017;Westin et al., 2016) combined with diffusion-relaxation correlations (de Almeida Martins et al., 2020;de Almeida Martins and Topgaard, 2018;Lampinen et al., 2019). We investigated the ability of neural networks to speed up model fitting, and explored the extent to which they can offset the requirements on the acquisition protocol.

Theory
White matter microstructure can be modelled by multiple compartments with different microstructural properties but a common orientation distribution (Alexander et al., 2019;Novikov et al., 2019). In this description, the measured signal is the convolution between an orientation distribution function (ODF) P (n) and a microstructural kernel K(u ⋅ n) where n and u are unit vectors defining the symmetry axes of the ODF and of the diffusion encoding process, respectively. Note that the microstructural kernel depends on the relative angle between n and u, cos β = u ⋅ n. In this work, we assign an effective transverse relaxation time T 2 and an apparent microscopic diffusion tensor D to each component, and use exponentially decaying functions to model the effect of these microstructural properties on the relaxation-diffusion-weighted signal . Under these assumptions, the microstructure kernel is written as a weighted sum of exponentials corresponding to a mixture of J components with signal fraction f j , transverse relaxation time T 2; j , and diffusion tensor D j . The colon ":" denotes the Frobenius inner product, B: D = ∑ i ∑ j B ij D ij . Information about T 2; j and D j is encoded into the signal by the echo time τ E and diffusion encoding tensor B(u), respectively, both of which are experimental variables. To simplify the model, we only consider axisymmetric B(u) and additionally assume that the component-wise D j are axisymmetric.
The convolution expressed in Eq.
(1) can be simplified by factorizing both P (n) and K(u ⋅ n) in their spherical harmonic coefficients p lm and k lm , respectively: and where Y lm are the spherical harmonics basis functions with the L l m (x) term denoting the associated Legendre polynomials. The summations in Eqs.
(4), we have taken the axial symmetry of the microstructural kernel K(u ⋅ n) into account (Lampinen et al., 2020;Novikov et al., 2018). Symmetry around the polar axis implies k l'm' = 0 for either m' ≠ 0 or odd l'. Taken together, this means that the k l'm' coefficients are reduced to their 0th degree terms k l'0 (typically written as k l' ) and only even-ordered spherical harmonic terms (l' = 0, 2, …) provide non-trivial contributions. Using the spherical harmonics addition theorem, Eq. (4) can be rewritten as Inserting Eqs.
The spherical harmonic coefficients of the microstructure kernel (k l0 ) and the ODF (p lm ) are estimated as the inner products between a given spherical harmonics basis function Y lm and either K(u ⋅ n) or P (n). Due to the orthonormality of the spherical harmonics basis, the inner products are given by multiplication with the complex conjugates of the Y lm , followed by integrations over the unit sphere. For the microstructural kernel, this procedure results in (Lampinen et al., 2020) k l0 ≡ k l = S 0 ∑ j = 1 J f j 4π(2l + 1)I lj exp −bD I; j 1 − b Δ D Δ; j exp − τ E T 2; j , where b is the conventional b-value and b Δ denotes the normalized anisotropy of the diffusion encoding tensor B (Eriksson et al., 2015). The isotropic diffusivity and the normalized diffusion anisotropy (D I and D Δ ) are related to the axial and radial diffusivities (D ∥ and D ⊥ ) of the diffusion tensor according to D I = D ∥ + 2D ⊥ /3 and D Δ = D ∥ − D ⊥ /3D I (Conturo et al., 1996); in its principal axis, a given D can thus be represented by a diagonal matrix parametrized as diag diag D I 1 − D Δ , D I 1 − D Δ , D I 1 + 2D Δ . The I lj factors are a function of the regular Legendre polynomials, L l , and defined as with α j = 3bD I; j b Δ D Δ; j . Different diffusion MRI models feature different numbers of components and impose different constraints on the component properties. Here we consider a two-compartment model (J = 2) comprising a "stick " component (S) with D Δ;S = 1 and a "zeppelin " (Z) component with D Δ;Z ≠ 1. Truncating the spherical harmonic summation at the second order (l max = 2) then yields the signal according to where m ∈ − 2, − 1, 0, 1, 2 . The derivation of Eq. (10) uses the p 00 = Y 00 = 1/ 4π ODF normalization (Lampinen et al., 2020;Novikov et al., 2018). The vectors e and m capture the experiment-related parameters, e = τ E , b, b Δ , θ, ϕ , and scalar model parameters, m = f S , D I;S , D I;z , D Δ;z , T 2; s , T 2; z , p 20 , Re p 21 , Im p 21 , Re p 22 , Im p 22 , where Re p lm = p lm + ( − 1) m p l − m /2 and Im p lm = p lm − ( − 1) m p l − m /2i denote the real and imaginary parts of the p lm coefficients, respectively. We refer to the model expressed by Eq. (10) as the Standard Model with Relaxation (SMR). This name is chosen to mark its descendance from the "standard model" of WM microstructure (Novikov et al., 2019) and to emphasize the fact that it accounts for compartment-specific T 2 times.
The SMR model parameters can be determined by fitting Eq. (10) directly to the acquired signals (Lampinen et al., 2020). An alternative strategy is to fit to some representation of the signal, such as the spherical harmonics coefficients. Veraart et al. (2018) used a model fitting framework that effectively reduces the dimensionality of the parameter space by means of performing a rotationally invariant factorization of the voxel-wise ODFs Reisert et al., 2017). The initial step of such framework consists in projecting the measured signal onto a spherical harmonic basis The S lm coefficients are subsequently converted to rotational invariants S l , and fitted to the corresponding rotationally invariant terms of the P (u) ⊗ K(g ⋅ u) convolution S l = p l k l , where k l is the 0th degree term of the microstructural kernel as defined by Eq. (8). The rotationally invariant coefficients, S l and p l , are computed from  x l = 4π∑ m x lm 2 (2l + 1) , where x lm are the spherical harmonics coefficients, and x l ≡ S l or x l ≡ p l . At sufficiently low b-values, signal projections with l > 2 have small contributions to the measured signal Jespersen et al., 2007) and the sum in Eq. (11) is typically truncated at the second order term (l = 2). The fitting framework summarized by Eqs. (11)-((13) is commonly referred to as the "RotInv " approach due to its use of rotational invariants. The l = 2 RotInv approach condenses the five p 2m , m ∈ − 2, − 1, 0, 1, 2 parameters of the SMR model onto a single p 2 invariant capturing the orientation coherence of the sub-voxel diffusion domains, thus reducing the dimensionality of the fitting problem by four parameters.

Neural network architecture and training
In this work, we constructed feedforward neural networks in MATLAB R2020b (The MathWorks, Inc.), and used them to fit vectors of scalar parameters, m = f S , D I;S , D I;Z , D Δ;Z , T 2; S , T 2; Z , p 20 , Re p 21 , Im p 21 , Re p 22 , Im p 22 to sets of measurements S τ E , B . We explored various network designs with different numbers of hidden nodes and/or layers before deciding on two final network architectures: an artificial neural network featuring 3 fully connected hidden layers with a decreasing number of nodes (180, 80, and 55) and a deeper/wider neural network featuring 4 fully connected hidden layers with 250 nodes each. All hidden layers were activated by hyperbolic tangent (tanh) functions and the deeper/wider network also featured batch normalization layers between the fully connected inner layers and their respective tanh activations. To distinguish the networks, we refer to them as the shallower neural network (SNN) and deeper neural network (DNN), respectively. Both SNN and DNN comprise an output layer with 11 nodes corresponding to the parameters in m. The input comprised a given number (E) of signal samples acquired with a pre-defined relaxation-diffusion encoding protocol. We considered three different acquisition protocols; with E = 164, E = 242, and E = 270 samples τ E , B .
Independent networks were trained for each protocol, meaning that 3 SNNs and 3 DNNs were evaluated. To remove the influence of S 0 from the fitting problem, we normalized the input vector to the median signal acquired at the lowest b-value and shortest echo-time.
Supervised network training was performed using a mean squared error loss where m targ is the ground-truth target vector, m net is the corresponding network output vector, and ∥ ⋅ ∥ 2 denotes the Euclidean norm. The m targ parameters were rescaled between 0 and 1 using a min-max normalization strategy before being supplied to the networks. The networks were trained with sets of voxels with randomly generated model parameters and noisy signal samples S τ E , B , as detailed in Section 3.2. The SNNs were trained with a batch size of 0.5•10 6 and a scaled conjugate gradient optimiser. The DNNs were trained in a mini-batch fashion using a total of 5•10 6 training sets, a mini-batch size of 50•10 3 , and an Adam optimiser. Throughout, training data was divided such that 75% of the original data was used to update the weights and biases and 25% was used for cross-validation. Overfitting was addressed by an early stopping method that terminated training following an increase of the MSE of the validation data for 5 (SNN) or 20 (DNN) consecutive epochs.
Network GPU training took approximately 83 min for the SNNs and 74 min for the DNNs on two parallel NVIDIA GeForce RTX 2080 SUPER, each with 8 GB of memory. Both graphic cards were installed on a high-end consumer-grade desktop computer with 32 GB memory and an 8-core Intel i9-9900k 3.6 GHz CPU with 2 threads per core.

Generating training data
We studied the impact of training data generation strategies on the network performance, including training based on uniformly sampled and real brain data. Training parameter vectors were created by two strategies: • m unif was synthetically constructed by random sampling of uncorrelated uniform distributions within the bounds described in Table 1; • m brain was constructed from in vivo brain data by randomly sampling parameter vectors estimated from a NLLS fit of Eq. (10). This dataset contains parameter correlations found in a typical brain dataset from a healthy adult.
The m brain vectors comprise the solutions of a nonlinear least-squares (NLLS) fit of Eq. (10) to in vivo signal data, referred to as m fit , together with an additional parameter set m mut , consisting of random mutations of the fitted solutions, given by where '∘' denotes the element-wise (Hadamard) product, and X is an 11-dimensional vector of normally distributed numbers. Each element of X is an independent and identically distributed random variable sampled from a normal distribution with mean 1 and standard deviation 0.3. The standard deviation of X was chosen following brief in silico experiments which revealed that virtually indistinguishable training/test results are obtained for standard deviations within the [0.2, 0.5] interval, provided all other training/ network parameters are kept constant. The number of m fit vectors was kept constant n fit ≈ 1.5 ⋅ 10 5 , and the total number of mutated vectors was defined as n mut = n brain − n fit .
The introduction of mutated parameters is a data augmentation technique, designed to simultaneously compensate for the relative low number of m fit vectors and expand the f S , D I;S , D I;Z , D Δ;Z , T 2; S , T 2; Z , p 20 , Re p 21 , Im p 21 , Re p 22 , Im p 22 domain of the m brain parameter targets.
The training vectors, m train , were combinations of m brain and m unif parameter vectors. Using a given total number of vectors n tot and varying number of m brain parameters (n brain ), we modulated the fraction of in vivo brain data, f brain = n brain /n tot , between 0 and 1 in steps of 0.05. The SNN training sets contained a total of n tot = 5 ⋅ 10 5 parameter vectors, while the DNN training sets contained n tot = 5 ⋅ 10 6 . Fig. S1 in the Supporting Information shows the distribution of m fit , m mut , and m unif parameters that compose a typical n tot = 5 ⋅ 10 5 SNN training dataset.
Signal data were generated from m train using Eq. (10) and one of three different τ E , B acquisition protocols:

•
The optimized protocol comprises tensor-valued encoding with full relaxationdiffusion-correlation optimized for minimal SMR parameter variance (Lampinen et al., 2020) • The unoptimized protocol comprises tensor-valued encoding with relaxationdiffusion-correlations restricted to low b-values (Lampinen et al., 2019). This protocol was an early attempt to design a diffusion-relaxation protocol with b-tensor encoding. It preceded the optimized protocol and was configured to fit into an available timeslot by following heuristics without a formal performance optimization, and was later found to yield degenerate results in white matter (Lampinen et al., 2020).

•
The LTE-only protocol comprises diffusion-relaxation optimized for minimal SMR parameter variance but includes only linear b-tensor encoding (b Δ = 1) (Lampinen et al., 2020). Just as the unoptimized protocol it has been found to yield degenerate solutions in white matter.
Additional details on the various protocols can be found in their respective references and in Table S1 of the Supporting Information. We emphasize that all training data used in this study was generated using the SMR forward model, Eq. (10), rather than using raw in vivo brain data.
Noise was sampled from the Rice distribution and added to the ground-truth synthetic signals. Because relaxation-diffusion MRI data comprises voxels with different signal-tonoise ratio (SNR), the amplitude of the SNR at S 0 = S τ E = 0, B = 0 was uniformly varied in the interval SNR ∈ [80,160]. Considering the relaxation-diffusion properties of typical healthy WM T 2 ≈ 70 ms, D I ≈ 0.9 μm 2 /ms , this choice results in SNR ∈ [30,60] at the point of maximal signal of the optimized protocol τ E = 63 ms, b = 0.1 ms/μm 2 , SNR amplitudes that are consistent with tensor-valued dMRI measurements of the in vivo brain (Szczepankiewicz et al., 2019a). Finally, networks were trained using m train vectors as targets and their corresponding in silico noisy signals as inputs.

Network evaluation
To find the optimal fraction of m unif and m brain parameters (adjusted by the f brain parameter), we trained SNNs with varying values of f brain , deployed them on in silico data generated from an unseen subject, and compared the various networks in terms of accuracy of the resulting parameter estimates. Network accuracy was assessed via normalized root-meansquared errors (NRMSE) and linear correlations with ground-truth values in terms of the Pearson correlation coefficient (ρ). The NRMSE captures the absolute agreement between the target ground-truth parameters and their corresponding network estimates, whereas ρ captures the linear target-to-estimate correlation strength. The f brain optimization process is discussed in more detail in section S3 of the Supporting Information. Briefly, the f brain hyper-parameter controls a trade-off between accuracy to WM-relevant parameters and network generalizability, and we found f brain = 0.5 to provide an optimal balance between accuracy and generalizability. From this point onward, we concentrate on networks trained with f brain = 0.5 datasets and evaluate them in further detail using correlation plots.
The accuracy performance of an f brain = 0.5 SNN, an f brain = 0.5 DNN, and a standard NLLS solver were compared on the basis of NRMSEs and Pearson correlation coefficients. The comparison was performed using two distinct in silico datasets: one based on m fit vectors from WM and deep GM data (m fit;WM-like ), and another based on m unif vectors. Each dataset comprised a total of 10•10 3 parameter vectors and their respective in silico signals.
The ground-truth synthetic signals were corrupted with Rician distributed noise and the SNR at the S 0 point was sampled uniformly from the [80, 160] range.
The effects of different acquisition protocols on network performance were evaluated in terms of NRMSE and sensitivity to parameter changes. The latter was gauged by modulating the non-orientational parameters of an SMR solution f S , D I;S , D I;Z , D Δ;Z , T 2; S , T 2; Z one at a time by 10% and measuring the response in all parameters. The original parameter set was based on in vivo data from the corona radiata where f S = 0.45, D I;S = 0.58 μm 2 /ms, D I;Z = 1.36 μm 2 /ms, D Δ;Z = 0.44, T 2;S = 69 ms, T 2;Z = 60 ms (Lampinen et al., 2020).
Subsequently, in silico datasets were generated for each of the 6 modulated datasets, Rice noise was added with SNR = 160 at S 0 , and parameter estimation was performed with protocol-specific networks.
To investigate if the reduced parameter space of RotInv fitting impacts the performance of ANN-based fitting, we trained an SNN using rotationally invariant in silico datasets and the same optimal f brain value found for the SMR networks. RotInv training vectors, m train;RI , were generated from the m train vectors (Section 3.2), using Eq. (13) to convert the full SMR parameters to RotInv parameters f S , D I;S , D I;Z , D Δ;Z , T 2; S , T 2; Z , p 2 . The RotInv in silico signal data was generated in four steps: (1) signals were calculated using m train and Eq. (10); (2) noise was added to the in silico signal data with a SNR ∈ [80, 160] at S 0 ; (3) S lm components were estimated by projecting the noisy S (τ E , B) signals to a spherical harmonics basis; and (4) S l, l = {0,2} signals were calculated from S lm using Eq.
(13). As with the full SMR model, training was performed using m train;RI as targets and their corresponding synthetic noisy signals as ANN inputs.
Trained SMR (RotInv) networks were tested on previously unseen m unif (m unif;RI ) and m brain (m brain;RI ) synthetic datasets at an SNR ∈ [80, 160] at S 0 . Performance was compared in terms of their respective target-estimate correlations. All networks were trained/ tested in a leave-one-out fashion where the training and testing m brain (m brain;RI ) datasets were generated using in vivo data from different subjects (Section 3.5).

In vivo data acquisition
We analysed data from three adult volunteers previously reported in (Lampinen et al., 2020). The study was approved by the regional ethical review board in Lund and written informed consent was obtained from all volunteers prior to scanning. Measurements were performed on a MAGNETOM Prisma 3T system (Siemens Health-care, Erlangen, Germany) using a prototype spin-echo EPI sequence that facilitates user-defined gradient waveforms for diffusion encoding (Szczepankiewicz et al., 2019a). Data were collected using echo times between 63 and 130 ms, repetition time of 3.4 s, voxel size of 2.5 mm 3 , 40 slices, matrixsize of 88 × 88, in-plane and through plane acceleration factor of 2 × 2 (GRAPPA), partial-Fourier of 3/4, band-width = 1775 Hz/pixel, and "strong " fat saturation. Diffusion encoding was performed with gradient waveforms optimized to maximize the encoding strength per unit time and to suppress concomitant field effects (Sjölund et al., 2015;Szczepankiewicz et al., 2019b). A total of 270 combinations of τ E and B were used, according to the optimized protocol in Table S1 of the Supporting Information. Total acquisition time was 15 min.

In vivo data processing and parameter estimation
Prior to analysis, all in vivo data were corrected for eddy-currents and subject motion using ElastiX (Klein et al., 2009) with extrapolated target volumes . Susceptibility-induced geometric distortions were corrected using the TOPUP tool in FMRIB software library (FSL) (Smith et al., 2004). Gibbs ringing artefact correction was performed according to the method described in (Kellner et al., 2016). To suppress the influence of noise, we filtered data with a 3D Gaussian kernel with a standard deviation of 0.45 times the voxel dimensions (Lampinen et al., 2020).
The SMR model parameters were estimated from a voxel-by-voxel NLLS fit of Eq. (10) to the post-processed data. The fitting process was performed with the multidimensional dMRI toolbox (https://github.com/markus-nilsson/md-dmri) (Nilsson et al., 2018), with MATLAB's built-in lsqcurvefit function. To remove outliers, model fitting was performed twice in each voxel and the result with lowest residual was retained (Lampinen et al., 2020). The initial guesses were sampled uniformly from the parameter bounds in Table 1. The resulting estimates were stored and used to compute in silico signal data following the procedure detailed in Section 3.2. NLLS fitting of a single in vivo brain dataset took approximately 8 h (approximately 5.5 s per voxel) on the CPU described in Section 3.1. The computations were carried out using parallel computing and multi-threading.
Finally, previously trained networks were used to estimate the parameters from Eq. (10) from in vivo data, which took approximately 2 and 20 s for the whole brain using the SNN and DNN, respectively. Training was performed on in silico m train data with an optimal f brain fraction. The training process followed a leave-one-out scheme, where the networks were trained on synthetic data generated from two subjects before being deployed/tested on a third, previously unseen, subject. Neural network fitting provided voxel-wise parameter maps that were compared to the ones obtained from a conventional NLLS fitting approach.

Neural network parameter estimates
SNN-based parameter estimation was approximately 10 4 times faster than NLLS fitting on the same computer, and yielded parameters in good agreement with the ground-truth targets and preserved contrast between regions characterized by distinct (T 2 , D) properties ( Fig.  1). For example, the estimated f S and p 2 are high in WM regions generally and highest in orientationally coherent WM regions such as the corpus callosum, similar to the in-silico ground-truth. However, a reduced contrast was observed in the T 2;Z maps, where the distinction between WM (darker) and cortical GM (brighter) regions is more prominent in the ground-truth map. The T 2;Z estimates are also characterized by considerable differences between ground-truth and estimated parameters in the long T 2 regions such as the lateral ventricles. The largest overall discrepancy between estimated and ground-truth parameters was found for D Δ;Z , likely because the signal is insensitive to it when |D Δ;Z | < 0.5 (Eriksson et al., 2015). Using an ANN trained on synthetic data to directly fit in vivo experimental data resulted in noisier maps. Nevertheless, it preserved an anatomically plausible contrast. Given the strong correlations between in silico ground-truth maps and network estimates, the noisier appearance of the in vivo parameter maps is likely because the SMR model cannot accurately represent the underlying in vivo data. In vivo SNN parameter estimates from WM regions of interests are displayed in Table S2 of the Supporting Information, where they are additionally compared to NLLS estimates. Fig. 2 shows that SNN-based estimates correlated well with the ground-truth parameter targets, with most parameters yielding linear correlation coefficients close to or above 0.9. The referenced figure focuses on the performance of a network trained with in silico S τ E , B data generated with the optimized protocol and an even mix of random and WM-like samples (f brain = 0.5), and distinguishes between performance on parameters obtained by uniform random sampling (light blue points) and parameters derived from in vivo non-cortical brain data (dark blue points). Red points correspond to parameter vectors derived from low component-specific signal fractions, as described in the figure caption. Poor performance is observed for low D Δ;Z values, where the network yields D Δ;Z ≈ 0.3 regardless of the underlying ground-truth. This can be attributed to an intrinsic difficulty in distinguishing between the diffusion-weighted signals of components with D Δ;Z < 0.5 components (Eriksson et al., 2015). Moreover, a poor target-to-estimate correspondence was seen for T 2;Z -times where these were longer than the maximal echo time.
The parameter maps estimated from the deeper network are in good agreement with their respective ground-truth targets (Figs. S3 and S4 of the Supporting Information correspond to Figs. 1 and 2). In Fig. S4, we observed that DNN-based fitting resulted in slightly stronger correlations between network estimates and ground-truth parameter targets. Although in vivo maps from DNN and SNN are similar, differences can be found in D Δ;Z and T 2;Z ; the DNN produces a noisier D Δ;Z and the T 2;Z map has a higher contrast between WM and cortical GM. Both of these features are likely artefactual, and suggest that the DNN is more susceptible to differences between the SMR signal predictions and the measured in vivo data.
The errors and prediction-target correlations of the ANN-based estimates are compiled in Table 2, where they are also compared to a conventional NLLS solver. The NLLS, SNN, and DNN approaches all have a comparable accuracy for in silico datasets designed to capture WM (T 2 , D) properties. By contrast, the function-fitting networks are observed to be more accurate than the NLLS approach for synthetic m unif parameter vectors.

Effect of acquisition protocol on network accuracy and sensitivity
In this section, we focus on the relationship between acquisition protocol design and network performance. Fig. 3 shows that ANN-based fitting could partly but not completely eliminate the known fit degeneracy in the unoptimized and LTE-only protocols: ANNs based on the optimized protocol provide lower estimation errors (NRMSE) than the ANNs based on the other two protocols. Comparing these two protocols, we note that the unoptimized protocol yields relatively more accurate estimates of D I;S , and T 2;Z , while the LTE-only protocol yields more accurate estimates of f S , D I;Z , D Δ;Z , T 2;S , and p 2 . Fig. 3 also shows that the performance of both SNN and DNN is less affected by sub-optimal acquisition protocols than the traditional NLLS approach (Fig. 3). For the NLLS approach, the use of the unoptimized or LTE-only protocols leads to a considerable increase of the estimation errors, while only a slight increase of NRMSE is observed for the DNN or SNN approaches. This suggests that ANNs may partly alleviate parameter estimation difficulties caused by a protocol that is inadequate in relation to the microstructure model. Fig. 4 shows the sensitivity of the various protocols to parameter changes. Networks trained on data generated with the optimized protocol are sensitive to all parameters, but slightly underestimate the magnitude of the change, particularly in D Δ;Z . The parameter-specific modulations did not have a major effect on the estimation of the remaining unmodulated parameters. An exception was found when the underlying T 2;Z is increased by 10%, which results in a 3% overestimation of the unchanged D I;S . Compared to the optimized protocol, the unoptimized and LTE-only protocols exhibit a lower sensitivity to the small parameter modulations and appear to be unresponsive to changes in D Δ;Z (both protocols) and D I;S (LTE-only). In addition to lower sensitivity, the unoptimized protocol also resulted in less accurate estimations of the unmodulated parameters, with a 10% modulation of f S leading to an erroneous 7% increase in T 2;S .

Neural network fitting of rotationally invariant microstructural features
Fig. 5 A shows that training a SNN with rotational invariants results in slightly stronger correlation between target and estimated parameters (compare with the scatter plots of Fig.  2). We note a considerable improvement in accuracy at low D Δ;Z values, where the constant D Δ;Z ≈ 0.3 behaviour observed for the full SMR model (see Fig. 2) is no longer present.
Applying the RotInv network to an unseen in vivo S l = {0,2} dataset results in parameter maps with anatomically plausible contrast (see Fig. 5 B). Consistent with the better D Δ;Z accuracy performance of the RotInv approach, we note that the RotInv D Δ;Z in vivo map has a smoother appearance and better demarcates cortical/non-cortical parenchyma than its non-rotationally invariant SMR counterpart (compare the fourth column of Figs. 1 with 5 B).
Interestingly, in vivo maps smoother than the ones displayed in Fig. 5 B can be attained from an ANN that was trained on unreasonably noisy in silico data. Fig. 6 displays the in vivo parameter maps obtained from a RotInv network trained with SNR ∈ [20, 40] at S 0 , which is 4 times lower than that used in Fig. 5. The resulting maps have a smooth appearance and exhibit anatomically plausible contrast. For example, regions with high f S correspond to WM regions, the lateral ventricles are characterized by low f S and high D I;Z values, and darker/brighter D Δ;Z regions demarcate cortical/non-cortical parenchyma. While it is tempting to favour the seductively 'robust' maps of Fig. 6 over the noisier maps of Fig.  5 B, we note that the low-SNR RotInv network results in weak correlations between target and estimated parameters (compare the scatter plots of Fig. 6 with those of Fig. 5 B). For example, SNN-based estimates of D Δ;Z may yield a smooth map that appears robust, but a closer inspection reveals that the D Δ;Z estimates in WM and deep GM regions are equal to the mean of the target D Δ;Z distribution and constitute an exceedingly inaccurate estimate of the underlying ground-truth. The tendency for networks to return the mean of the training parameter distribution has been reported in studies of the RotInv model (Reisert et al., 2017) and the behaviour was explained in detail by Coelho et al. (2021).

Discussion and conclusions
Replacing traditional NLLS solvers with function-fitting neural networks enables vastly faster parameter estimation when using high-dimensional microstructural models. On a consumer-grade desktop computer, the fitting time was reduced from hours (NNLS) to seconds (ANN). Naturally, the NNLS fitting times, based on the relatively slow trust-regionreflective algorithm, can be improved by linearizing the fitting problem (Daducci et al., 2015) or by using GPU-based solvers (Harms et al., 2017). However, while such procedures have enabled whole brain fitting of non-linear models within minutes (Daducci et al., 2015;Harms et al., 2017), we still expect the seconds-long forward pass of an ANN to provide a competitive choice in terms of computation time.
The ANN-based estimates were observed to be in good agreement with synthetic data that mimicked healthy WM as well as data that spanned the entire space of allowed model parameters. When deployed on unseen in vivo brain data, neural networks provide maps that are consistent with known brain anatomy and preserve contrast be-tween regions with different relaxation-diffusion properties. Our findings are encouraging and in line with recent advanced dMRI modelling studies that use machine learning techniques for parameter estimation (Barbieri et al., 2020;Bertleff et al., 2017;Golkov et al., 2016;Grussu et al., 2020;Gyori et al., 2019;Hill et al., 2021;Kaandorp et al., 2021;Nedjati-Gilani et al., 2017;Palombo et al., 2020;Reisert et al., 2017). A combination of error metrics, correlation analysis, and sensitivity matrices was found to provide a useful set of tools for quantitatively assessing parameter-specific accuracy/sensitivity and for identifying the limitations of learning-based approaches. These tools facilitate a survey of the performance across all dimensions of the SMR model, for example, revealing that D Δ;Z was consistently less accurate than other parameters, as expected from previous studies that have emphasized that it is difficult to estimate (Eriksson et al., 2015;Lampinen et al., 2020Lampinen et al., , 2019. By contrast, visual inspection of ANN-based parameter maps was found to provide limited insight on the general performance of the networks. Indeed, smooth and anatomically plausible maps can be achieved even with poor network performance and data with low SNR. This is a common and deceptive pitfall that has strong implications for the evaluation of performance in machine learning approaches (Reisert et al., 2017).
We found no evidence that voxel-wise ANN-based parameter estimation can fully alleviate the degenerate fitting landscape typically present when working with biophysical models in dMRI (Jelescu et al., 2016) or replace an exhaustive sampling of all relevant experimental dimensions (Coelho et al., 2019;Lampinen et al., 2020). Fig. 3 shows worse performance in terms of parameter estimation errors (NRMSE) for the two protocols with known degeneracy problems (Lampinen et al., 2020). Similarly, Fig. 4 shows that only the optimized protocol can faithfully recover parameter-specific changes while the other two cannot. These are both signs of unresolved degeneracies. Indeed, we cannot expect good performance for LTE-only and unoptimized protocols because these protocols can yield virtually identical signal vectors for different model parameters; the inverse problem has many solutions (Lampinen et al., 2020(Lampinen et al., , 2019. Nevertheless, ANN-based fitting showed an advantage compared with the traditional NLLS approach, as it yielded lower estimation errors in the degenerate cases (unoptimized and LTE-only protocols; Fig. 3). Our de Almeida Martins et al. Page 14 interpretation is that the NLLS approach returns one out of the many solutions, whereas the ANN-based estimate tends to an average across the many solutions.
The 11-dimensional parameter space of the SMR model is difficult to sample densely and thus presents a challenge when designing training datasets that are representative of the vast fitting landscape. In this work, we addressed this challenge by constructing training data based on in vivo healthy brain data (m brain ) and more naïve parameter vectors randomly sampled from the entire model parameter space (m unif ). Networks trained exclusively with m brain vectors displayed the best accuracy in terms of expected WM properties, but their domain of validity is restricted to the relatively small space spanned by m brain solutions. This raises questions about their generalizability, i.e., their performance in cases where atypical microscopic tissue structures are present (Alexander et al., 2019). To find a good trade-off between accuracy and generalizability we optimized the fraction of in vivo-based training data (f brain ). However, we expect that more work is needed to define a truly optimal strategy for network training.
In this study, we focused on fully connected networks that follow the design of multilayer perceptrons (MLPs), a traditional ANN class that is well-suited for regression problems (Cybenko, 1989;Hornik et al., 1989). Alternatives or complements to the fully connected ANN architecture should also be explored in future works. Promising avenues include the use of dropout (Gal and Ghahramani, 2016;Tanno et al., 2021) or deep ensemble strategies (Lakshminarayanan et al., 2016;Qin et al., 2021) as a means to derive uncertainty metrics, the use of rolled-out network structures inspired by non-learning-based iterative fitting frame-works (Ye, 2017), the use of auto-encoders (Zucchelli et al., 2021), or the use of denoising networks (Fadnavis et al., 2020;Wang et al., 2019) to minimize the amount of noise present in the data that is supplied to the function-fitting ANN. While fundamentally different network architectures may considerably boost the performance of the ANN-based fitting approach, the modest differences found between the SNN and DNN designs suggests that simply increasing the width and/or depth of the fully connected ANN architecture is not a promising avenue. Despite the potential for improvement, we note that the plots in Fig. 2 constitute an improvement over similar target-estimate correlation plots reported in (Reisert et al., 2017), where supervised learning based on polynomial regressors was used to fit a three-compartment diffusion model, and are equivalent to the correlation plots reported in more recent works on learning-based fitting of diffusion (Gyori et al., 2019;Palombo et al., 2020) and diffusion-relaxation MRI models (Grussu et al., 2020).
The fully connected ANNs we considered here are not invariant to sample rotations. The input to the ANNs is a vector of E signal samples measured at a pre-defined set of both directional (θ, ϕ) and rotationally invariant τ E , b, b Δ experimental points. The ordering in which the E measurements are provided is kept fixed and samples with similar microstructural kernels K(u ⋅ n) but different orientations will result in distinctive network input vectors. This places a burden on the training data generation, which has to span a sufficient set of possible tissue orientations. Using the RotInv formulation renders the network invariant to sample rotations, which considerably reduces the dimensionality of the parameter space that has to be represented in training data. The higher training efficiency likely explains its slightly higher accuracy performance relative to the full SMR networks. de Almeida Martins et al. Page 15 Alternatives to the RotInv formulation presented in this study include a framework based on a different set of rotationally invariant features of the dMRI signal (Zucchelli et al., 2021) or fitting the full SMR model with equivariant network architectures (Cohen et al., 2018;Thomas et al., 2018).
A potential limitation of the present study is the focus on a single multi-compartment model of tissue microstructure whose range of application is mostly limited to WM and deep GM tissues. Applications for cortical GM should therefore consider models tailored to the appropriate microstructure . Our decision to focus on a single model follows from previous dMRI literature which has presented the "Standard Model " of tissue microstructure -from which our SMR model descends -as an overarching signal model that encompasses several other WM models as particular cases (Novikov et al., 2019). Fur-thermore, the "Standard Model" has been used to reveal general degeneracy problems in microstructure parameter estimation . Given the generality of our model and the prevalence of degeneracies in advanced dMRI modelling, we expect the degradation of performance with less optimal protocols to also be found in alternative multi-compartment models or when using different learning-based fitting algorithms (e.g. : polynomial Reisert et al. 2017 or random forest Nedjati-Gilani et al. 2017. However, future work is needed to fully characterize the general relationship be-tween machine learning approaches and degenerate fitting landscapes.
In conclusion, function fitting neural networks can be used to vastly accelerate parameter estimation with high-dimensional microstructural MRI models. The accuracy of ANN-based estimates was observed to degrade less with sub-optimal protocols than traditional NLLS fitting. However, the performance of function fitting networks was still observed to primarily depend on the amount of information sampled by the underlying measurements, and we found no evidence that ANN-based approaches can offset the need for a rich set of data. Therefore, machine learning methodology in MRI microstructure modelling should be matched with comprehensive data acquisition. This work presents a learning-based fitting framework, as well as tools for evaluating combinations of networks and measurement protocols in terms of error metrics, estimate-target correlation plots, and sensitivity matrices.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material. Deploying trained networks on previously unseen in silico and in vivo data provides anatomically plausible parameter maps in under 10 s (including data management times). The first and second columns compare the ground-truth targets and network predictions, respectively, of the in silico dataset. Difference maps are shown in the third column. Parameter maps obtained from applying a trained network on in vivo brain data are displayed in the fourth column. Scatter plots of ground-truth parameters vs. neural network predictions. Light blue points show results when the network is deployed on uniformly distributed random model parameters. The dark blue points correspond to an in silico dataset derived from a nonlinear least-squared fit to measured brain data where voxels within CSF and cortical GM were excluded by masking out regions where microscopic anisotropy (Lasič et al., 2014), μFA, is lower than 0.6. The red points correspond to regions where poor accuracy is expected, i.e., where the signal fraction of the relevant component ("stick " or "zeppelin " depending on the parameter) accounts for less than 15% of the total signal or, for the p 2 map, parameter vectors where the "zeppelin " component accounts for more than 85% of the total signal fraction and | D Δ;Z | < 0.4. The inner legends show the Pearson correlation coefficients (ρ) of the blue points. For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article. Optimized acquisition protocols result in ANN-and NLLS-based parameter estimates with smaller errors. The bar plots indicate the normalized root-mean-squared errors (NRMSE) between ground-truth and predicted parameters, for learning-based (DNN and SNN) and NLLS fitting approaches, and for in silico datasets generated with different acquisition protocols. The leftmost plots correspond to a tensor-valued τ E , B protocol optimized for minimal parameter variance, the optimized protocol (Lampinen et al., 2020); the middle plots correspond to a sub-optimal tensor-valued τ E , B protocol where relaxation-diffusion correlations are exclusively established at low b-values, the unoptimized protocol (Lampinen et al., 2019); the rightmost plots show the results for a τ E , B protocol optimized for minimal parameter variance when limited to linear diffusion encoding b Δ = 1 , the LTEonly protocol (Lampinen et al., 2020). Panel A shows network performance on parameters sampled from a uniform distribution, and panel B shows the performance on in silico data based on least-squares fitting results to in vivo non-cortical brain tissue data. Sensitivity of acquisition protocols to 10% parameter modulations. The matrices display the relation between an induced parameter change and the observed response. When a single parameter on the y-axis is modulated by 10%, the response can be read in all other parameters along the x-axis. An ideal network would report a diagonal matrix with the value 10% on the diagonal, and zero otherwise. The optimized protocol appears sensitive in all parameters, whereas the unoptimized protocol lacks sensitivity D Δ;Z and the LTE-only protocols lacks sensitivity to both D Δ;Z and D I;S .  Neural network fitting of a rotationally invariant (RotInv) model results in strong target-estimate correlations and plausible maps. (A) Correlations between network-based parameter estimates and ground-truth parameter targets. The estimates were obtained from a RotInv network trained using a fraction of f brain = 0.5 between rotationally invariant m brain and m unif training parameter vectors. The colour-coding and legends follow the same convention as Fig. 2. (B) Maps of microstructural diffusion parameters -f S , D I;S , D I;Z and D Δ;Z -obtained from fitting a RotInv network to rotationally invariant in vivo brain data.   SMR parameter bounds. The diffusivity bounds were enforced by limiting D ∥ ; s , D ∥ ; Z and D ⊥ ; Z to the [0.2, 4.0] μm 2 /ms interval. For T 2;S and T 2;Z , the lower bound removes the influence of the assumedly fully-attenuated myelin water, and the large upper bound of T 2;Z enables it to capture effects of increased values in white matter lesions (Lampinen et al., 2019) as well as possible contamination with cerebrospinal fluid which is expected to have a larger influence on the more isotropic zeppelin compartment (Lampinen et al., 2020).  Table 2 Accuracy performance of DNN-, SNN-and NLLS-based fitting approaches. Performance is evaluated on synthetic data simulated from two different sets: uniformly sampled random parameters (m unif ), and parameters derived from least-squared model fitting to in vivo WM and deep GM data (m fit;WM-like ).