Patryk Filipiak
email: patryk.filipiak@inria.fr

Rutger Fick

Alexandra Petiet

| Mathieu Santin

Anne-Charlotte Philippe

Stephane Lehericy

Philippe Ciuciu

Rachid Deriche

| Demian Wassermann

Reducing the Number of Samples in Spatio-Temporal dMRI Acquisition Design

Keywords: diffusion mri, acquisition design, stochastic optimization

Purpose: Acquisition time is a major limitation in recovering brain white matter microstructure with diffusion magnetic resonance imaging. The aim of this paper is to bridge the gap between growing demands on spatio-temporal resolution of diffusion signal and the real-world time limitations. The authors introduce an acquisition scheme that reduces the number of samples under adjustable quality loss.

Methods: Finding a sampling scheme that maximizes signal quality and satisfies given time constraints is NP-hard.

Therefore, a heuristic method based on genetic algorithm is proposed in order to find sub-optimal solutions in acceptable time. The analyzed diffusion signal representation is defined in the qτ space, so that it captures both spacial and temporal phenomena.

Results:

The experiments on synthetic data and in vivo diffusion images of the C57Bl6 wild-type mouse corpus callosum reveal superiority of the proposed approach over random sampling and even distribution in the qτ space.

Conclusion:

The use of genetic algorithm allows to find acquisition parameters that guarantee high signal reconstruc-1 2 PATRYK FILIPIAK ET AL. tion accuracy under given time constraints. In practice, the proposed approach helps to accelerate the acquisition for the use of qτ -dMRI signal representation.

momentum and propose, for the first time, the qτ -based acquisition design solutions.

The main goal of our study is to find a sampling scheme that maximizes accuracy of a signal representation and satisfies given time constraints. Similar studies were performed for Diffusion Kurtosis Imaging (DKI) [START_REF] Hansen | Experimentally and computationally fast method for estimation of a mean kurtosis[END_REF][START_REF] Poot | Optimal experimental design for diffusion kurtosis imaging[END_REF] or the composite hindered and restricted model of diffusion (CHARMED) [START_REF] Alexander | A general framework for experiment design in diffusion MRI and its application in measuring direct tissuemicrostructure features[END_REF][START_REF] Alexander | Orientationally invariant indices of axon diameter and density from diffusion MRI[END_REF][START_REF] Santis | Improved precision in CHARMED assessment of white matter through sampling scheme optimization and model parsimony testing[END_REF] including gradient waveform optimization [START_REF] Drobnjak | Optimising time-varying gradient orientation for microstructure sensitivity in diffusionweighted MR[END_REF][START_REF] Drobnjak | Optimizing gradient waveforms for microstructure sensitivity in diffusion-weighted MR[END_REF]. In this paper, we optimize for the accuracy of qτ -dMRI signal reconstruction. Additionally, we want our approach to be usable on real-world applications, considering that different acquisition protocols are recommended for different WM microstructure models [START_REF] Sepehrband | A time-efficient acquisition protocol for multipurpose diffusion-weighted microstructural imaging at 7 Tesla[END_REF]. For this, we discretize the spatio-temporal search space by performing a dense pre-acquisition of dMRI signal with uniform coverage of the unit sphere, as suggested by Caruyer et al. [START_REF] Caruyer | Design of multishell sampling schemes with uniform coverage in diffusion MRI[END_REF].

Despite discretization, the problem is computationally difficult. It requires selecting an optimal subset of Diffusion Weighted Images (DWIs), which is NP-hard, as we will show in Section 3.

Taking into account that the time complexity of our problem grows exponentially with the domain size, such that global optima cannot be found deterministically within few hours or even few days, we apply a stochastic search engine instead. We use Standard Genetic Algorithm (SGA) [START_REF] Holland | Adaptation in natural and artificial systems. An introductory analysis with application to biology, control, and artificial intelligence[END_REF][START_REF] Pena-Reyes | Evolutionary computation in medicine: an overview[END_REF][START_REF] Sabat | Three dimensional k-space trajectory design using genetic algorithms[END_REF] for this purpose, which allows us to find approximate solutions in acceptable time.

Our experiments comprise of two granularity levels with respect to the domain size. The coarse-grained level allows us to present the effectiveness of the proposed approach by comparing our results with the global optimum found with lengthy exhaustive search. The fine-grained level outcomes provide us with crucial information about the structure of the optimized acquisition schemes.

We validate our approach on both synthetic diffusion model and real data comprising in vivo diffusion images of the C57Bl6 wild-type mouse corpus callosum.

| THEORY

In a nutshell, our spatio-temporal signal representation is a cross-product of a 3D space of diffusion gradient directions and a 1D range of diffusion times, hence 4D. This section briefly introduces the mathematical formulation of the suggested representation. For more details we refer to Fick et al. [START_REF] Fick | Multi-Spherical Diffusion MRI: Exploring Diffusion Time Using Signal Sparsity[END_REF][START_REF] Fick | Non-parametric graphnet-regularized representation of dMRI in space and time[END_REF][START_REF] Fick | MAPL: Tissue microstructure estimation using Laplacian-regularized MAP-MRI and its application to HCP data[END_REF].

| Four-dimensional Ensemble Average Propagator

Let us first define the qτ -diffusion signal space and its relationship to the 4D Ensemble Average Propagator (EAP) [START_REF] Callaghan | Pulsed-gradient spin-echo NMR for planar, cylindrical, and spherical pores under conditions of wall relaxation[END_REF][START_REF] Stejskal | Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient[END_REF]. In dMRI, the EAP P (R; τ ) describes the probability density that a particle undergoes a displacement R ∈ R 3 after diffusion time τ > 0. The EAP is estimated from a set of DWIs, which are obtained by applying two sensitizing diffusion gradients G ∈ R 3 of pulse length δ > 0, separated by separation time ∆ > 0. Assuming that no diffusion takes place during the pulses (δ → 0), the EAP is related to the dMRI signal through the Inverse Fourier Transform (IFT) [START_REF] Stejskal | Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient[END_REF] defined as

P (R; τ ) = R 3 E(q, τ )e i2πq•R dq, (1)
where the signal attenuation E(q, τ ) = S(q, τ )/S 0 , and S(q, τ ) is the signal measured at the diffusion encoding position q ∈ R 3 and diffusion time τ = ∆ -δ/3, whereas S 0 is the baseline image acquired without diffusion sensitization, i.e. q = 0. We denote q = |q| and G = |G|, such that q = qu and G = Gv for some 3D unit vectors u, v ∈ S 2 . The wave vector q in Equation ( 1) is defined as q = γδG/2π, where γ is the nuclear gyromagnetic ratio and G is the applied diffusion gradient vector.

| qτ -space signal representation with GraphNet regularization

We reconstruct the continuous EAP from a finite set of DWIs by representing the discretely measured attenuation E(q, τ ) in terms of the basis coefficients c of a "Multi-Spherical" 4D qτ -Fourier basis [START_REF] Fick | Multi-Spherical Diffusion MRI: Exploring Diffusion Time Using Signal Sparsity[END_REF]. The qτ -basis is formed by the cross-product of a 3D q-space basis Φ i (q) [START_REF] Özarslan | Mean apparent propagator (MAP) MRI: A novel diffusion imaging method for mapping tissue microstructure[END_REF] and 1D diffusion time basis T j (τ ) [START_REF] Fick | Non-parametric graphnet-regularized representation of dMRI in space and time[END_REF]. The approximated signal attenuation Ê(q, τ, c) is given as

Ê(q, τ, c) = Nq i=1 Nτ j=1 c ij Φ i (q)T j (τ ) with c = [c ij ] ∈ R Nq×Nτ , (2)
where Nq and Nτ are the maximum expansion orders of spatial and temporal bases, respectively, and c ij are the weights of the contribution of the ij th basis function to Ê(q, τ, c). As Φ is a Fourier basis over q, the EAP can be recovered for each voxel as P (R; τ, c) = IFTq Ê(q, τ ; c) .

To estimate c from a noisy and sparsely sampled E(q, τ ), we use GraphNet regularization [START_REF] Grosenick | Interpretable whole-brain prediction analysis with Graph-Net[END_REF], as advised by Fick et al. [START_REF] Fick | Non-parametric graphnet-regularized representation of dMRI in space and time[END_REF]. Note that the commonly assumed Rician noise perturbing the dMRI signal can be approximated with Gaussian distribution for SNR > 5 [START_REF] Gudbjartsson | The Rician distribution of noisy MRI data[END_REF]. In our GraphNet approach, we use the Laplacian regularization term to compensate for the Gaussian noise and the l1-norm regularization to impose sparsity over the basis coefficients. Mathematically, the above are formulated as follows

argmin c E(q, τ ) -Ê(q, τ, c) 2 dqdτ + λ ∇ 2 Ê(q, τ, c) dqdτ + α c 1 subject to Ê(0, τ, c) = 1 and Ê(q, 0, c) = 1. (3)
The parameters λ, α > 0 stand for the smoothness and sparsity regularization weights, respectively. We optimize them using five-fold cross-validation.

Having introduced our 4D dMRI signal representation, let us now focus on the acquisition design, which is the main contribution of this paper. Note that our optimization methodology is not strictly related to qτ -dMRI. In fact, one can replace it with any other approach as long as it allows to recover the original dMRI signal from the subsampled DWIs.

| METHODS

We aim to find a sampling scheme that satisfies a given time constraint under adjustable quality loss. For this, we first perform a dense acquisition of dMRI signal. Then, among those densely acquired samples, we seek a fixed-sized subset for which our qτ -dMRI signal representation reaches the highest accuracy taking the dense acquisition as reference.

In the following subsections, we first formulate this goal mathematically as an optimization problem and study its complexity. Later on, we address the problem with a Genetic Algorithm (GA). Then, we define the accuracy measures and describe the experiments.

| Optimal Acquisition Design

Let N > 0 be the number of DWIs in a dense pre-acquisition. As a matter of fact, each DWI corresponds with a certain (q, τ ) pair, although for the moment it is more convenient to think of them as an enumerated set of samples. Among those, we want to select a subset of up to nmax < N samples, such that the reconstruction accuracy of the dMRI signal is maximized. To this end, we define the objective function F : {0, 1} N → R in the space of binary vectors x = (x 1 , ..., x N ) ∈ {0, 1} N . The assignment x i = 1 for a given i = 1, ..., N indicates that the i-th sample from the pool of DWIs is included in the subset of interest, whereas x i = 0 determines its exclusion. Formally, our goal is to solve the following optimization problem arg min

x

F (x) = 1 M M j=1 E (j) - Ê(j) x 2 2 subject to N i=1 x i ≤ nmax with 1 ≤ nmax ≤ N -1, (4)
where M > 0 is the number of voxels in each DWI, E (j) is the normalized signal attenuation captured in the j-th voxel with the qτ measurements, and Ê(j)

x is the corresponding signal obtained by fitting the qτ -dMRI representation (defined in Equation 3) to the subsamples indicated by x. Note that the normalization of E (j) is necessary to equalize signal intensities obtained with different (q, τ ) parameters.

From now on, we will omit the voxel indexing (j) while referring to E and Êx for simplicity. Our goal thus comes down in one sentence: minimize the residuals between E and Êx using up to nmax samples.

The problem posed in Equation ( 4), despite its simple formulation, is in fact very difficult to solve analytically. In Combinatorics, it is known under the name Knapsack Problem (KP) and is proved to be NP-hard [START_REF] Hochbaum | Approximation algorithms for NP-hard problems[END_REF]. Let us remind that the objective of KP is to pick a finite set of items that maximize the total value of the knapsack, while respecting its capacity limitation. In our case, the goal is to select a subset of qτ -indexed samples that maximize accuracy of dMRI signal recovery, while satisfying given time constraints. Assuming that acquisition time of each DWI is constant, we express the time budget as the maximum number of qτ -indexed samples in a subset, nmax.

The only known deterministic solver of KP is the exhaustive search procedure comprising N !/(nmax!•(N -nmax)!) evaluations of the objective function F , which is unacceptable for large N . Thus, we address the problem with a stochastic search engine. For this purpose, we choose SGA [START_REF] Goldberg | Genetic algorithms in search, optimization, and machine learning[END_REF][START_REF] Holland | Adaptation in natural and artificial systems. An introductory analysis with application to biology, control, and artificial intelligence[END_REF] due to its reported high performance in large KP instances [START_REF] Chu | A genetic algorithm for the multidimensional knapsack problem[END_REF][START_REF] Khuri | The zero/one multiple knapsack problem and genetic algorithms[END_REF][START_REF] Michalewicz | Genetic algorithms for the 0/1 knapsack problem[END_REF][START_REF] Olsen | Penalty functions and the knapsack problem[END_REF][START_REF] Spillman | Solving large knapsack problems with a genetic algorithm[END_REF]. Formally, the time complexity of SGA is polynomial, although practically it is difficult to estimate, since certain input parameters, including the number of iterations, need to be chosen experimentally [START_REF] Holland | Adaptation in natural and artificial systems. An introductory analysis with application to biology, control, and artificial intelligence[END_REF].

| Objective function & performance measures

As defined in Equation ( 4), the objective of our optimization mechanism is to minimize the residuals between E and ÊX .

We use this quantity as a primary measure of the dMRI signal reconstruction accuracy. Additionally, we are interested in verifying how well the spatio-temporal indices are preserved while using our scheme. To this end, we study a set of commonly used metrics [START_REF] Fick | MAPL: Tissue microstructure estimation using Laplacian-regularized MAP-MRI and its application to HCP data[END_REF][START_REF] Özarslan | Mean apparent propagator (MAP) MRI: A novel diffusion imaging method for mapping tissue microstructure[END_REF], namely

• Return To Origin Probability:

RTOP(τ ) = P (0; τ ) [1/m 3 ]
• Return To Axis Probability:

RTAP(τ ) = R P (Rr ; τ )dR [1/m 2 ]
• Return To Plane Probability:

RTPP(τ ) = R {r∈S 2 :r•r =0} P (Rr ⊥ ; τ )dr ⊥ dR [1/m]
• Mean Squared Displacement:

MSD(τ ) = R S 2 P (Rr; τ )R 2 drdR [m 2 ]
for a given displacement R = Rr, as defined in Section 2. Let us mention that two of the above metrics, i.e. RTAP and RTPP, assume that WM is modeled by parallel cylinders with the vectors r parallel and r ⊥ perpendicular to the cylinder axis.

Granularity levels:

The acquisition schemes that we study are 800-dimensional binary feature vectors. Dealing with such high-dimensional data involves the risk of bias towards a dense pre-acquisition scheme or over-fitting to dMRI signal noise. In order to avoid these pitfalls, we consider the following two granularity levels:

• Optimization by shells -where each of the shells in our multishell dense acquisition is either taken as a whole or discarded completely. In this scenario, the task is to subsample among 5 available τ values and 8 available G values, which mimics the experiment described by Alexander et al. [START_REF] Alexander | A general framework for experiment design in diffusion MRI and its application in measuring direct tissuemicrostructure features[END_REF][START_REF] Alexander | Orientationally invariant indices of axon diameter and density from diffusion MRI[END_REF]. Particularly, for the budget size nmax = 100, the search space is reduced to 658,008 possible solutions, hence we are able to find the global optimum of the problem by performing exhaustive search and use it as a reference (cf. Figure 1). Our proposed optimizer uses a set of 100 candidate solutions per iteration in this scenario, which is a commonly used setting for SGA [START_REF] Baluja | Removing the genetics from the standard genetic algorithm[END_REF].

• Optimization by measures -where solutions are picked from all the 800 variants of densely acquired samples.

Considering relatively high complexity of this scenario, we chose experimentally the population size of our SGA optimizer to be 1000 candidate solutions per iteration.

| Experiments setup

In the experiments, we analyze both synthetic and real diffusion data using the protocol described below.

Dense pre-acqusition:

Time constraints:

We considered four variants of time limits expressed as budget sizes nmax ∈ {100, 200, 300, 400} out of 800 DWIs. The remaining 40 b0-images were excluded from the optimization domain, as they were used by default in every acquisition scheme.

Other approaches:

We compared our method with two alternative sampling schemes. One of them, called random, used the uniform random distribution of qτ samples in the index space {1, ..., N }. In the second one, referred to as even, we picked each i-th sample for i = kN/nmax and k = 1, ..., nmax. Considering that the samples were ordered by G and ∆, the even subsampling scheme ensured nearly equal number of samples per shell, i.e. nmax/40 or nmax/40 .

| Diffusion data

We used the following two data sets:

Synthetic data:

As first data set, we generated time-dependent diffusion data using our Python-based [START_REF] Fick | An Open-source Framework for Reproducible dMRI-Based Microstructure Research (Version 0.1)[END_REF][START_REF] Fick | Mipy: An Open-Source Framework to improve reproducibility in Brain Microstructure Imaging[END_REF] implementation of the two-compartment model with the intra-cellular fraction modeled with a set of Watson dispersed cylinders [START_REF] Callaghan | Pulsed-gradient spin-echo NMR for planar, cylindrical, and spherical pores under conditions of wall relaxation[END_REF][START_REF] Zhang | NODDI: practical in vivo neurite orientation dispersion and density imaging of the human brain[END_REF] and the extra-cellular one modeled with a temporal zeppelin [START_REF] Burcaw | Mesoscopic structure of neuronal tracts from time-dependent diffusion[END_REF]. The Watson distribution is an antipodally symmetric distribution, centered around the principal orientation µ, describing a density controlled by concentration parameter κ > 0, which is inversely related to the dispersion of cylinders [START_REF] Callaghan | Pulsed-gradient spin-echo NMR for planar, cylindrical, and spherical pores under conditions of wall relaxation[END_REF][START_REF] Zhang | NODDI: practical in vivo neurite orientation dispersion and density imaging of the human brain[END_REF]. Ronen et al. [START_REF] Ronen | Microstructural organization of axons in the human corpus callosum quantified by diffusion-weighted magnetic resonance spectroscopy of N-acetylaspartate and postmortem histology[END_REF] modeled the axon dispersion in the corpus callosum with a spherical Gaussian distribution and found out that typical values of its standard deviation vary between 10 and 25 degrees, which translates to the interval 3 < κ < 16. In another study, Zhang et al. [START_REF] Zhang | NODDI: practical in vivo neurite orientation dispersion and density imaging of the human brain[END_REF] used the values of κ ∈ {0, 0.25, 1, 4, 16}. In our experiment, we chose κ = 4 from the intersection of both sets.

Table 1 summarizes all the model's parameters that we took from Burcaw et al. [START_REF] Burcaw | Mesoscopic structure of neuronal tracts from time-dependent diffusion[END_REF], Alexander et al. [START_REF] Alexander | Orientationally invariant indices of axon diameter and density from diffusion MRI[END_REF], and Zhang et al. [START_REF] Zhang | NODDI: practical in vivo neurite orientation dispersion and density imaging of the human brain[END_REF].

[Table 1 about here.]

Apart from the original noiseless signal, we also studied the two variants of the signal with incorporated Rician noise, having respective Signal-to-Noise Ratio (SNR) set to 20 (13dB) and 10 (10dB).

Real data:

In the second data set, we used in vivo diffusion images of the corpus callosum of C57Bl6 wild-type mouse.

| RESULTS

Exhaustive search:

As we mentioned earlier, we defined the coarse-grained optimization by shells scenario in such a way that it could be solved with exhaustive search and serve as a reference. The running time of a parallel code for this approach was about one week. Let us mention that using the same exhaustive search implementation for the optimization by measures scenario would take more then 10 200 centuries. [Figure 1 about here.]

In silico experiments:

As shown in Figure 2, convergence of SGA in the optimization by shells scenario is apparently reached before the 30th iteration. Based on that observation, we cautiously chose the termination condition for all experiments to be the fixed threshold of 40 iterations, thus leaving a safety margin of 10 iterations.

[Figure 2 about here.] Figure 3 (left-column plots) and Table 2 summarize the results of the tested subsampling schemesours, random, and even -for the in silico experiment, averaged over 30 runs. In the rows of Table 2, we present respectively the noiseless case and the two variants with incorporated Rician noise having SNR = 20 (13dB) and 10 (10dB). Within each case, we consider four time budget limits nmax ∈ {100, 200, 300, 400}. The scores are expressed with normalized root mean squared errors (NRMSEs) and the corresponding standard deviations (STDs). The STD is omitted in the even scheme which is deterministic and thus its STD=0.

[Figure 3 about here.]

As Table 2 shows, our approach outperforms the other two in all studied cases, in both granularity levels, reaching lowest error values and standard deviations. It is also worth noticing that the addition of Rician noise induces overfitting of the qτ -dMRI signal representation for nmax > 200. As a consequence, our approach handles better the noisy data in the optimization by shells scenario, having 20 times less free parameters, than the optimization by measures.

We compared all the pairs of results, i.e. ours vs. random and ours vs. even, in both granularity levels, using paired two-sample Student's t-tests with the Bonferroni adjusted significance level α = 10 -5 and the number of degrees of freedom 2n -2 = 58. In each case, our approach was statistically significantly better than the other two schemes.

[Table 2 about here.]

In order to study the stability of our proposed method, we distilled a single best acquisition scheme found in each of the 30 runs of SGA. Next, we arranged these solutions by MSEs and selected the top 10% out of them. Figure 4 illustrates the averages of those schemes. Note the perfect similarity between the best schemes found by SGA, presented in Figure 4(a), and the global optimum depicted in Figure 1(a). Furthermore, the concentrations of high likelihood (red squares) in Figure 4 tend to form shapes that are especially visible in the optimization by shells scenario. The latter suggests that SGA repeatedly converges to the same or highly similar solutions in each run, thus giving stable outputs.

[Figure 4 about here.]

In vivo experiments: 3 present the results for the C57Bl6 wild-type mouse corpus callosum (CC).

Table 3 details the three regions of interest, i.e. genu, body, and splenium, and the time budgets nmax ∈ {100, 200, 300, 400}. Our method outperforms the other two significantly, which we verified using the same type of Student's t-test as above. However, unlike in silico experiment, there is no superiority of optimization by shells over optimization by measures here. In fact, the residuals in both scenarios are comparable, with slightly lower averages in optimization by measures for small nmax, whereas optimization by shells gives lower averages for large nmax.

[Table 3 about here.]

Analogously to the in silico experiment, we illustrate the averaged top 10% acquisition schemes obtained with SGA in the form of heat maps depicted in Figure 5. Here also the gradually smoothing patterns are present as in Figure 4.

[Figure 5 about here.]

Spatio-temporal indices: As it is seen in Figure 6, all the indices are well preserved under subsampling by the three tested schemes in the noiseless case of in silico experiment. However, the discrepancies between the reconstructed and the reference curves readily increase with the addition of noise. Also, index dispersion grows as SNR decreases and diffusion time increases.

The differences between the reconstructed and the reference curves are similar in the in vivo experiment presented in Figure 7. Generally, all the tested schemes give comparable outputs, although the results of our approach are less dispersed in this case. The highest accuracy is reached at the splenium region, while there are considerable errors in the body.

[Figure 6 about here.]

[Figure 7 about here.]

| DISCUSSION

It is beyond doubt that acquisition time matters, especially in long-lasting protocols like dMRI. In this paper, we propose a method to shorten acquisition time with reduced signal loss by optimizing the sampling scheme in the (q, τ ) parameter space. Additionally, we want our scheme to preserve temporal changes in diffusion signal, hence we study the qτ -indices presented in Section 3.

In practice, our approach helps to accelerate the acquisition for the use of qτ -dMRI signal representation. It is crucial, since the latter allows to extrapolate the signal outside qτ -samples.

| Our approach maximizes signal reconstruction accuracy

The accuracy of WM microstructure recovery grows monotonically with an increase of dMRI acquisition density, although the contributions of particular DWIs to the measured signal attenuation are not equal. This indicates the need for identifying a variable-density for acquiring the dMRI signal. Our method, described in Section 3, seeks a fixed-sized set of samples that provide the most accurate signal reconstruction with respect to the reference data obtained with the dense acquisition scheme. In other words, our approach allows to reach high accuracy within a controllable time budget.

As we can see, by comparing Figures 1 and4, our approach is able to find nearly-optimal solutions in less then 30 iterations of the proposed algorithm (cf. Figure 2). The accuracy of dMRI signal reconstruction obtained with our technique is significantly better than random or even subsampling of the densely acquired signal, as we can observe in Tables 2 and3, as well as Figure 3 5

.2 | Scheme optimization is most profitable in short acquisition times

In our experiments, we studied the four time budgets nmax ∈ {100, 200, 300, 400} out of 800 densely acquired samples.

Both, in silico (cf. Table 2 and Figure 3) and in vivo experiments (cf. Table 3 and Figure 3), showed that superiority of our technique was most evident in the cases with tightest time constraints nmax = 100, which is 1/8 of the original time span. For nmax ≥ 200, the differences in residuals of the dMRI signal reconstruction among the three tested approaches decreased gradually. Also, overfitting of the qτ -dMRI signal emerged as nmax grew. This means that the less imaging time we had, the more we benefited from the acquisition scheme optimization. On the other hand, when small speed-ups were concerned (e.g. 1/2 of the original time span), our optimization seemed unnecessary. The way to go in such cases was as simple as to remove gradient directions from the dMRI scanning protocol evenly, i.e. every second one for 1/2 time save, every third one for 1/3, and so on.

| The optimized schemes form repetitive patterns

In agreement with Caruyer et al. [START_REF] Caruyer | Design of multishell sampling schemes with uniform coverage in diffusion MRI[END_REF], we believe there can be no single all-purpose optimal acquisition scheme. The imaging parameters depend on numerous factors including the needs of a given WM microstructure recovery method, the time constraints, and the physical limitations of a scanner. That is why we propose a data-driven mechanism which is able to adapt to varying imaging conditions by maximizing the accuracy of the signal reconstructed with qτ -dMRI.

Alexander [START_REF] Alexander | A general framework for experiment design in diffusion MRI and its application in measuring direct tissuemicrostructure features[END_REF] and Caruyer et al. [START_REF] Caruyer | Design of multishell sampling schemes with uniform coverage in diffusion MRI[END_REF] postulated that an angular coverage of gradient directions across the sphere should be "as uniform as possible" in order to ensure orientational invariance of a scheme. The results of our study partly comply with those remarks. Note that the additional parameters which controlled the distribution of gradient directions introduced in the optimization by measures scenario improved the reconstruction accuracy in the noiseless cases only, compared to the coarse-grained optimization by shells. Apparently, the ad hoc alterations that we made possible after freeing the space of gradient directions were of little importance for the signal reconstruction accuracy.

We thus conclude that studying our problem at such a fine-grained level might be an overhead in real-world applications.

On the other hand, there is a possibility that choosing a different (i.e. non-uniform) scheme for the dense pre-acquisition would lead to slightly different conclusions in this regard, which leaves a room for further improvement of the proposed technique.

Apparently, setting up the parameters G and ∆ is more difficult than handling gradient directions. An analogue of the simplistic "as uniform as possible" strategy, that we had called the even scheme in our experiments, turned out insufficient in this case. Tables 2 and3 show that the above approach gave mediocre results, roughly comparable to the naive random subsampling scheme.

The heat maps presented in Figures 4 and5 shed some light on the optimized schemes obtained with our approach.

The distribution of parameters forms clusters of points in the (G, ∆) space, especially in the optimization by shells scenario. In the optimization by measures, both in silico and in vivo, the dominance of low G-values emerges as nmax increases.

In the similar experiment described by Alexander [START_REF] Alexander | A general framework for experiment design in diffusion MRI and its application in measuring direct tissuemicrostructure features[END_REF], the distribution of pulse lengths ∆ was close to uniform, although it came with the dominance of largest G-values which is not the case in our study. The above discrepancy is most probably caused by the differences in objective functions. Alexander's work targeted particular biomarkers like axon density and radius in WM, whereas we maximize reconstruction accuracy of a dMRI signal using our qτ signal representation. This in turn brings us again to the earlier remark that the notion of acquisition scheme optimality is highly problem-dependent.

| Globally optimal acquisition schemes may be beyond reach

It is proved that the SGA optimization engine, that we use here, will always converge to global optimum when given enough time [START_REF] Davis | A Markov chain framework for the simple genetic algorithm[END_REF]. Nevertheless, the question "how much time is enough?" is generally unanswerable, as it requires thorough analysis of each problem in hand.

For better understanding of our parameter space, we illustrate the full spectrum of feasible acquisition schemes at the top row of Figure 1. The solutions are arranged from best to worst, i.e. from the lowest MSE of signal reconstruction on the left-hand side boundary to the highest one on the right. Such arrangements is convenient to visualize the accuracy of potential solutions, although we must point out that it disregards the space topology. On the top-left plot in Figure 1, there is a narrow valley in the proximity of global optimum and the high peek on the other side of the plot. Between them, there is a large interval of low steepness. Such a curve shape suggests that the globally best solution can be difficult to find using randomized approaches, since it is not surrounded by a considerable neighborhood of nearly-optimal ones. Indeed, the curve of the normalized Hamming distance to the global optimum (in the feature space), depicted at the top-right plot in Figure 1, grows rapidly with the quantiles pc. To make matters worse, the low steepness in the middle of the top-left plot indicates a possibly large number of local minima or plateau areas which in turn would hinder the numerical optimization techniques like GAs. Both those observations suggest that with an increase of problem complexity, chances of finding global optima will decrease rapidly.

Our experiments showed that SGA was able to find the global optimum in the optimization by shells scenario with synthetic diffusion data and nmax = 100 (cf. Figures 1 and4). However, it is very likely that in difficult problems like optimization by measures, we found sub-optimal solutions instead.

| Our approach preserves spatio-temporal indices when reducing acquisition time

As it is seen in Figures 6 and7, the four analyzed spatio-temporal indices are preserved when using our approach. The even subsampling scheme produces comparably good averages, however the curves are much more dispersed. On the other hand, the random subsampling scheme often gives closer approximations, albeit disregarding the shape of the curve producing a linear fit. Our approach better reproduces the shape of the reference curves although with an increased bias. Additionally, our method performs better on the diffusion tensor indices, particularly Fractional Anisotropy, presented in Supporting Information Figures S1 andS2.

The observed agreement with the reference qτ indices, which were not included in the objective function introduced in Section 3, suggests that our method is able to generalize well. With nmax = 100, we avoided overfitting to noise that was present in the densely pre-acquired signal. On the other hand, it might be interesting to extend the objective function, by adding multiple criteria covering, for instance, mean squared displacement or fractional anisotropy.

We also find it promising that the spatio-temporal signal representation itself turned out to harmonize efficiently with our approach. Note that all the plots of qτ indices presented here illustrate the results for the tightest time budget nmax = 100. This implies again that a great deal of acquisition time can be saved without much compromise on dMRI signal accuracy.

The studied diffusion times, 10-20ms, cover the range of ∆ values that were feasible for the pulsed-gradient spin echo mode due to the T2 decay time. Note that the qτ indices computed with our approach smoothly interpolate between sampled diffusion times. Also, our qτ -dMRI signal representation is able to extrapolate the index values outside the measured time interval, although the accuracy in such cases requires further investigation.

| CONCLUSIONS

We proposed the spatio-temporal dMRI acquisition design that greatly reduces the number of qτ samples under the adjustable quality loss. Despite the fact that selecting a sampling scheme that maximizes brain white matter reconstruction accuracy and satisfies given time constraints is NP-hard, our stochastic optimization mechanism based on genetic algorithm found sub-optimal solutions efficiently.

The experiments on both synthetic diffusion data and in vivo images of the C57Bl6 wild-type mouse corpus callosum revealed superiority of our technique over random subsampling and even distribution in the qτ space. Our approach performed best under the tightest among all the considered time constraints, leading to reduction of acquisition time to 1/8 of the original time span. Additionally, we observed repetitive patterns in our optimized schemes that gave us crucial insight about the space of acquisition parameters. The latter is really promising for future identification of optimal sampling density.

In this study, we assumed availability of a densely acquired dMRI signal for reference, although it is not often the case. Nonetheless, we believe that our preliminary work will allow to uncover the optimal sampling density in the (q, τ )

space for optimal recovery of the EAP, and thus the optimal range of acquisition parameters G and ∆.

Future work should target the reproducibility of our approach among different subjects and scanners. For instance, the risk of peripheral nerve stimulation while using high gradient magnitudes should be taken into account. Also, different areas of WM tissue need to be studied, including for instance crossing fibers. Finally, the optimizer itself might be improved to ensure faster convergence and adaptability, and thus achieve lower average quality loss of solutions.

AC K N O W L E D G E M E N T S

This Convergence of SGA

F I G U R E 2
The convergence of our method is reached before the 30th iteration of the algorithm. The plots present the convergence curves averaged over 30 runs of SGA in the optimization by shells scenario for the in silico experiment with the Rician noise, SNR=20db (the upper plot), and the C57Bl6 wild-type mouse corpus callosum body (the lower plot). The time budgets nmax ∈ {100, 200, 300, 400} are color-coded. TA B L E 3 Our approach significantly outperforms the other two subsampling schemes (with p-value < 10 -5 ) in both granularity levels (by shells and by measures), reaching lower means and standard deviations. The tables present the residuals of the dMRI signal reconstruction for the three regions of C57Bl6 wild-type mouse corpus callosum (CC) and the time budgets nmax ∈ {100, 200, 300, 400}. The results are expressed as normalized root mean squared errors (NRMSEs) of signal reconstruction with standard deviations (STDs) aggregated over 30 runs. The STD is omitted in the even scheme which is deterministic and thus its STD=0.

Figure 1

 1 Figure 1 illustrates all the 658,008 feasible acquisition schemes in the optimization by shells scenario for the in silico experiment with nmax = 100. The schemes are arranged from best to worst. The top-left plot shows the mean squared errors (MSEs) of signal reconstruction, whereas the top-right one depicts the normalized Hamming distances from the global optimum ± 1 standard deviation. In order to visualize the analyzed (G, ∆) parameter space, the percentiles pc = 0%, 1%, 10%, 50%, 90% are annotated on both plots, showing respectively the global optimum, the top 1% solutions, the top 10% solutions, etc. The cumulative averages of acquisition schemes, corresponding with these percentiles, are shown in the heat maps at the bottom. The colors reflect the likelihood of a given (G, ∆) pair in the scheme. The heat maps (a) and (b) represent, respectively, the global optimum and its proximity. The interval of percentiles between pc = 10% and pc = 90%, as shown in the heat maps (c)-(e), contains a spectrum of feasible acquisition schemes with similar MSEs and almost equally large Hamming distances from the global optimum. Note the arrangements of high likelihood (red-colored squares) in Figure 1(a). There is one red field in the area of highest G-values located in the middle of the ∆ range. Then, there are four other red fields spread evenly across G-and ∆ parameter spaces. The pattern gradually disappears as we move away from the global optimum. Starting from about 10th percentile, the differences in intensities on the heat map slowly fade.

Figure 3 (

 3 Figure 3 (right-column plots) and Table3present the results for the C57Bl6 wild-type mouse corpus callosum (CC).

Figures 6 and 7

 7 Figures6 and 7present the reconstruction of spatio-temporal indices RTOP, RTAP, RTPP, MSD ± 1 standard deviation with nmax = 100, obtained for in silico and in vivo experiments respectively. The black lines show the reference curves obtained from the densely acquired signal, whereas the color lines represent even (blue), random (green), and ours (red) subsampling schemes. In addition, Supporting Information FiguresS1 and S2present the diffusion tensor related curves of mean, axial and radial diffusivity, and fractional anisotropy for in silico and in vivo experiments respectively.

  work has received funding from the ANR/NSF award NeuroRef; the European Research Council (ERC) under the Horizon 2020 research and innovation program (ERC Advanced Grant agreement No 694665 : CoBCoM); the MAXIMS grant funded by ICM's The Big Brain Theory Program and ANR-10-IAIHU-06; the programs "Institut des neurosciences translationnelle" ANR-10-IAIHU-06 and "Infrastructure d'avenir en Biologie Santé" ANR-11-INBS-0006. the least dispersed outputs. The plots present the reconstruction of spatiotemporal indices RTOP, RTAP, RTPP, MSD ± 1 standard deviation with nmax = 100, obtained for the three regions of C57Bl6 wild-type mouse corpus callosum (CC). The black plots show reference curves obtained from the densely acquired signal, the color plots represent even (blue), random (green), and ours (red) subsampling schemes. As expected, the RTOP, RTAP, and RTPP indices are decreasing, whereas MSD is increasing with diffusion time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 S1 Our approach ensures the least dispersed outputs. The plots present the reconstruction of Mean Diffusivity, Axial Diffusivity, Radial Diffusivity, and Fractional Anisotropy indices ± 1 standard deviation with nmax = 100 in the optimization by measures, obtained for the in silico experiment either with or without Rician noise. The black plots show reference curves obtained from the densely acquired signal, the color plots represent even (blue), random (green), and ours (red) subsampling schemes. . . . . . . . . . 25 S2 Our approach ensures the least dispersed outputs. The plots present the reconstruction of Mean Diffusivity, Axial Diffusivity, Radial Diffusivity, and Fractional Anisotropy indices ± 1 standard deviation with nmax = 100, obtained for the three regions of C57Bl6 wild-type mouse corpus callosum (CC). The black plots show reference curves obtained from the densely acquired signal, the color plots represent even (blue), random (green), and ours (red) subsampling schemes. . . . . . . . . . . . . . . . . . . . . . . . 26 %n )cheme) arranged fr%m be)t t% w%r)t [% rank] pc = 50% pc = 90% Exhaustive search results F I G U R E 1 Exhaustive search results of the optimization by shells for the in silico experiment with nmax = 100. The plots at the top present all the 658,008 feasible acquisition schemes arranged from best to worst, illustrating the mean squared errors (MSEs) of signal reconstruction (top-left plot) and the normalized Hamming distances from the global optimum ± 1 standard deviation (top-right). In order to visualize the analyzed (G, ∆) parameter space, the percentiles pc = 0%, 1%, 10%, 50%, 90% are annotated on both plots, showing respectively the global optimum, the top 1% solutions, the top 10% solutions, etc. The corresponding cumulative averages of acquisition schemes are depicted in the heat maps at the bottom. The colors reflect the likelihood of a given (G, ∆) pair in the scheme. The heat maps for pc ≤ 0% and pc ≤ 1% represent, respectively, the global optimum and its proximity. The interval between pc = 10% and pc = 90% contains a huge spectrum of schemes with similar MSEs and almost equally large distances from the global optimum.

3

 3 Our approach significantly outperforms the other two subsampling schemes (with p-value < 10 -5 ) in both granularity levels: optimization by shells (top-row plots) and optimization by measures (bottom-row plots), reaching lower means and standard deviations. The plots present the residuals of the dMRI signal reconstruction for the in silico experiment with the dispersion controlled by the concentration parameter κ = 4 and SNR=20 (13dB

  s(0(on sche)es ((n 2(2o e3pe.()en0)F I G U R E 5The concentrations of high likelihood (red squares) tend to form consistent shapes that are especially visible in the optimization by shells scenario. The plots present the averages of the top 10% acquisition schemes found by SGA in the optimization by shells scenario (top row) and the optimization by measures scenario (bottom row) for the body region of the C57Bl6 wild-type mouse corpus callosum, and the time budgets nmax ∈ {100, 200, 300, 400}. The colors reflect the likelihood of a given (∆, G) pair in the scheme.F I G U R E 6 All the tested schemes give comparable index curves for the in silico experiment. The plots present the reconstruction of spatio-temporal indices RTOP, RTAP, RTPP, MSD ± 1 standard deviation with nmax = 100. The black plots show reference curves obtained from the densely acquired signal, the color plots represent even (blue), random (green), and ours (red) subsampling schemes. As expected, the RTOP, RTAP, and RTPP indices are decreasing, whereas MSD is increasing with diffusion time.F I G U R E 7 Our approach ensures the least dispersed outputs. The plots present the reconstruction of spatio-temporal indices RTOP, RTAP, RTPP, MSD ± 1 standard deviation with nmax = 100, obtained for the three regions of C57Bl6 wild-type mouse corpus callosum (CC). The black plots show reference curves obtained from the densely acquired signal, the color plots represent even (blue), random (green), and ours (red) subsampling schemes. As expected, the RTOP, RTAP, and RTPP indices are decreasing, whereas MSD is increasing with diffusion time.S U P P O R T I N G I N F O R M AT I O N F I G U R E S 1Our approach ensures the least dispersed outputs. The plots present the reconstruction of Mean Diffusivity, Axial Diffusivity, Radial Diffusivity, and Fractional Anisotropy indices ± 1 standard deviation with nmax = 100 in the optimization by measures, obtained for the in silico experiment either with or without Rician noise. The black plots show reference curves obtained from the densely acquired signal, the color plots represent even (blue), random (green), and ours (red) subsampling schemes.S U P P O R T I N G I N F O R M AT I O N F I G U R E S 2Our approach ensures the least dispersed outputs. The plots present the reconstruction of Mean Diffusivity, Axial Diffusivity, Radial Diffusivity, and Fractional Anisotropy indices ± 1 standard deviation with nmax = 100, obtained for the three regions of C57Bl6 wild-type mouse corpus callosum (CC). The black plots show reference curves obtained from the densely acquired signal, the color plots represent even (blue), random (green), and ours (red) subsampling schemes.L I S T O F TA B L E S 1Parameters used for generating the in silico diffusion data with the values taken from the referenced literature. Only the principal orientation was chosen arbitrary. The Watson distribution is an antipodally symmetric distribution, centered around the principal orientation µ, describing a density controlled by concentration parameter κ > 0, which is inversely related to the dispersion of cylinders. . . .

  Our initial dense pre-acquisition covered 40 shells, each of which comprised 20 directions and one b 0 -image, i.e. to parametrize our (q, τ ) acquisition space. The corresponding b-values ranged from 41 to 5248 s/mm 2 . The gradient duration δ = 5 ms remained constant throughout the experiments. In each of the 8 × 5 = 40 shells, we followed the uniform distribution of directions suggested by Caruyer et al.[START_REF] Caruyer | Design of multishell sampling schemes with uniform coverage in diffusion MRI[END_REF].

	40 × 20 = 800 DWIs plus 40 non-weighted images. We used combinations of 8 gradient strengths G ∈ {50, 100,
	150, 200, 250, 300, 350, 400} [mT/m] (magnitudes available in pre-clinical scanners) and 5 separation times ∆ ∈ {10.8,
	13.1, 15.4, 17.7, 20.0}[ms]

  ) (left-column plots) and the body part of C57B16 wild-type mouse corpous callosum (right-column plots), for the time budgets nmax ∈ {100, 200, 300, 400}. The results are expressed as normalized root mean squared errors (NRMSEs) of signal reconstruction with standard deviations (STDs) aggregated over 30 runs. The concentrations of high likelihood (red squares) tend to form consistent shapes that are especially visible in the optimization by shells scenario. The plots present the averages of the top 10% acquisition schemes found by SGA in the optimization by shells scenario (top row) and the optimization by measures scenario (bottom row) for the in silico experiment with the time budgets nmax ∈ {100, 200, 300, 400}. The colors reflect the likelihood of a given (∆, G) pair in the scheme.

			Be01 c.2(0(1(on 0cheme0 ((n 0()(co e3pe/(men1)
	400 350 300 250 200 150 100 50 400 350 300 250 200 150 100 50 350 300 250 200 150 100 50 400 350 300 250 200 150 100 50 F I G U R E 4 400 OPTIMIZATION BY SHELLS gradient strengths G [mT/m] OPTIMIZATION BY MEASURES gr d(en1 01reng1h0 G [mT/m] OPTIMIZATION BY SHELLS gradient strengths G [mT/m] OPTIMIZATION BY MEASURES g. d(en0 s0.eng0hs G [)T/m]	(a) n max = 100 10.8 13.1 15.4 17.7 20.0 0ep r 1(on 1(me0 4 [m0] (e) n max = 100 (a) n max = 100 10.8 13.1 15.4 17.7 20.0 (e) n max = 100 sep . 0(on 0()es 4 [)s]	0.00 0.01 0.02 0.03 0.04 0.05 0.00 0.01 0.02 0.03 0.04 0.05 0.00 0.01 0.02 0.03 0.04 0.05 0.00 0.01 0.02 0.03 0.04 0.05	(b) n max = 200 10.8 13.1 15.4 17.7 20.0 0ep / 1(on 1(me0 4 [m0] (f) n max = 200 (b) n max = 200 10.8 13.1 15.4 17.7 20.0 (f) n max = 200 sep . 0ion 0i)es 4 [)s]	0.00 0.01 0.02 0.03 0.04 0.05 0.00 0.01 0.02 0.03 0.04 0.05 0.00 0.01 0.02 0.03 0.04 0.05 0.00 0.01 0.02 0.03 0.04 0.05	(c) n max = 300 10.8 13.1 15.4 17.7 20.0 0ep / 1(on 1(me0 4 [m0] (g) n max = 300 (c) n max = 300 10.8 13.1 15.4 17.7 20.0 (g) n max = 300 sep . 0(on 0()es 4 [)s]	0.00 0.01 0.02 0.03 0.04 0.05 0.00 0.01 0.02 0.03 0.04 0.05 0.00 0.01 0.02 0.03 0.04 0.05 0.00 0.01 0.02 0.03 0.04 0.05	(d) n max = 400 10.8 13.1 15.4 17.7 20.0 0ep / 1(on 1(me0 4 [m0] (h) n max = 400 (d) n max = 400 10.8 13.1 15.4 17.7 20.0 sep . 0(on 0()es 4 [)s]	0.00 0.01 0.02 0.03 0.04 0.05 0.00 0.01 0.02 0.03 0.04 0.05 0.00 0.01 0.02 0.03 0.04 0.05

  . . . . . .28 2 Our approach significantly outperforms the other two subsampling schemes (with p-value < 10 -5 ) in both granularity levels (by shells and by measures), reaching lower means and standard deviations. The tables present the residuals of the dMRI signal reconstruction for the in silico experiment and the time budgets nmax ∈ {100, 200, 300, 400}. The results are expressed as normalized root mean squared errors (NRMSEs) of signal reconstruction with standard deviations (STDs) aggregated over 30 runs. The STD is omitted in the even scheme which is deterministic and thus its STD=0. . . . . . . . . . . . . . . . . . . . . 29 3 Our approach significantly outperforms the other two subsampling schemes (with p-value < 10 -5 ) in both granularity levels (by shells and by measures), reaching lower means and standard deviations. The tables present the residuals of the dMRI signal reconstruction for the three regions of C57Bl6 wildtype mouse corpus callosum (CC) and the time budgets nmax ∈ {100, 200, 300, 400}. The results are expressed as normalized root mean squared errors (NRMSEs) of signal reconstruction with standard deviations (STDs) aggregated over 30 runs. The STD is omitted in the even scheme which is deterministic and thus its STD=0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 Parameters used for generating the in silico diffusion data with the values taken from the referenced literature.Only the principal orientation was chosen arbitrary. The Watson distribution is an antipodally symmetric distribution, centered around the principal orientation µ, describing a density controlled by concentration parameter κ > 0, which is inversely related to the dispersion of cylinders. Our approach significantly outperforms the other two subsampling schemes (with p-value < 10 -5 ) in both granularity levels (by shells and by measures), reaching lower means and standard deviations. The tables present the residuals of the dMRI signal reconstruction for the in silico experiment and the time budgets nmax ∈ {100, 200, 300, 400}. The results are expressed as normalized root mean squared errors (NRMSEs) of signal reconstruction with standard deviations (STDs) aggregated over 30 runs. The STD is omitted in the even scheme which is deterministic and thus its STD=0.

	intra-cellular fraction parameter's name value fraction weight 0.3 principal orientation µ azim.: π/2, elev.: 0 parallel diffusivity λ 1.7 × 10 -9 m 2 /s cylinder diameter d 1.0 × 10 -6 m concentration parameter κ 4 TA B L E 1 time noise budget optimization by shells ref. [11] [2] [2] [60] NRMSE ± STD ×10 -3 extra-cellular fraction parameter's name value fraction weight 0.7 principal orientation µ azim.: π/2, elev.: 0 parallel diffusivity λ 1.7 × 10 -9 m 2 /s bulk diffusivity constant λ∞ 6.5 × 10 -10 m 2 /s [11] ref. [11] [2] characteristic coefficient A 7.41 × 10 -12 m 2 [11] optimization by measures nmax ours random even ours random even 100 17.3 ± 0.72 43.7 ± 22.70 23.5 14.9 ± 0.14 23.3 ± 3.23 28.1 noiseless 200 13.2 ± 0.26 21.8 ± 4.85 15.7 12.7 ± 0.13 15.8 ± 0.71 21.9 signal 300 12.1 ± 0.15 18.5 ± 4.35 16.4 12.2 ± 0.07 14.5 ± 0.53 14.0 400 11.8 ± 0.06 16.2 ± 2.55 13.1 11.7 ± 0.08 14.0 ± 0.46 14.3 100 165.6 ± 1.06 225.4 ± 40.15 192.1 166.6 ± 0.41 184.4 ± 2.49 186.6 SNR=20 200 161.5 ± 0.39 196.8 ± 29.45 190.4 162.7 ± 0.21 177.9 ± 3.03 183.8 (13dB) 300 158.6 ± 0.15 184.3 ± 18.75 171.1 160.5 ± 0.18 173.2 ± 1.75 172.4 400 157.3 ± 0.18 174.5 ± 9.16 169.6 159.0 ± 0.22 170.2 ± 2.12 171.9 100 275.6 ± 1.51 329.0 ± 28.08 304.5 273.6 ± 0.64 298.5 ± 4.64 304.7 SNR=10 200 266.8 ± 0.48 301.1 ± 20.03 299.5 268.4 ± 0.43 288.8 ± 3.41 297.8 (10dB) 300 263.7 ± 0.28 291.2 ± 16.55 281.3 266.5 ± 0.97 283.6 ± 3.47 283.0 400 262.1 ± 0.51 283.7 ± 13.34 279.7 264.6 ± 0.18 279.4 ± 2.24 280.3 time NRMSE ± STD ×10 -3 of budget optimization by shells optimization by measures interest nmax ours random even ours random even 100 110.2 ± 0.45 131.6 ± 14.98 120.1 106.9 ± 0.33 118.9 ± 1.95 120.8 CC 200 103.7 ± 0.17 112.7 ± 4.18 112.1 102.8 ± 0.14 109.8 ± 1.10 114.1 genu 300 101.7 ± 0.14 108.3 ± 3.11 106.9 101.5 ± 0.09 106.4 ± 0.78 106.5 400 100.6 ± 0.05 104.9 ± 1.75 105.4 100.8 ± 0.13 104.1 ± 0.52 104.1 100 127.6 ± 0.85 172.2 ± 44.30 144.2 122.9 ± 0.54 142.1 ± 3.28 155.3 CC 200 116.1 ± 0.45 136.2 ± 9.83 132.7 116.6 ± 0.19 128.4 ± 1.76 149.7 body 300 111.9 ± 0.36 124.8 ± 6.41 122.0 113.4 ± 0.20 121.4 ± 1.73 124.1 400 109.9 ± 0.22 118.9 ± 3.13 119.1 111.5 ± 0.10 116.8 ± 1.16 119.9 100 107.0 ± 1.09 133.3 ± 19.75 123.1 104.4 ± 0.37 117.9 ± 2.80 131.1 CC 200 101.4 ± 0.41 111.2 ± 4.82 109.6 100.9 ± 0.11 108.4 ± 1.32 141.0 splenium 300 99.9 ± 0.42 105.7 ± 2.23 106.8 99.6 ± 0.09 104.2 ± 0.75 104.4 TA B L E 2 region 400 99.0 ± 0.42 103.2 ± 1.46 103.7 98.6 ± 0.12 101.9 ± 0.53 103.8

Funding information the ANR/NSF award NeuroRef; the European Research Council (ERC) under the Horizon 2020 research and innovation program (ERC Advanced Grant agreement No 694665 : CoBCoM); the MAXIMS grant funded by ICM's The Big Brain Theory Program and ANR-10-IAIHU-06; the programs "Institut des neurosciences translationnelle" ANR-10-IAIHU-06 and "Infrastructure d'avenir en Biologie Santé" ANR-11-INBS-0006.

L I S T O F F I G U R E S

Exhaustive search results of the optimization by shells for the in silico experiment with nmax = 100.

The plots at the top present all the 658,008 feasible acquisition schemes arranged from best to worst, illustrating the mean squared errors (MSEs) of signal reconstruction (top-left plot) and the normalized Hamming distances from the global optimum ± 1 standard deviation (top-right). In order to visualize the analyzed (G, ∆) parameter space, the percentiles pc = 0%, 1%, 10%, 50%, 90% are annotated on both plots, showing respectively the global optimum, the top 1% solutions, the top 10% solutions, etc. The corresponding cumulative averages of acquisition schemes are depicted in the heat maps at the bottom.

The colors reflect the likelihood of a given (G, ∆) pair in the scheme. The heat maps for pc ≤ 0% and