Pierre Ablin

Jean-François Cardoso

Alexandre Gramfort

Spectral independent component analysis with noise modeling for M/EEG source separation

Keywords: ICA, EEG, MEG, source separation

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Magnetoencephalography and Electroencephalography (M/EEG) are popular non-invasive techniques to record brain activity [START_REF] Hämäläinen | Magnetoencephalographytheory, instrumentation, and applications to noninvasive studies of the working human brain[END_REF][START_REF] Niedermeyer | Electroencephalography : basic principles, clinical applications, and related fields[END_REF]. They capture respectively the magnetic and electric signals produced by active neurons from the scalp surface or close to it. Each M/EEG sensor captures a combination of the different brain activities. The physics of the mixing is well understood: it is a linear process and can be considered instantaneous.

Independent Component Analysis (ICA) [START_REF] Hyvärinen | Independent component analysis[END_REF] is extensively used in neuroscience for processing M/EEG signals [START_REF] Makeig | Independent component analysis of electroencephalographic data[END_REF]. In its simplest form, it models the observed signals as a linear combination of statistically independent signals called sources. Remarkably, ICA can identify these sources 'blindly', that is, without prior knowledge of the underlying physics of the system (except linearity).

Applied on EEG signals, it separates meaningful brain signals from artifacts

(eye blinks, heartbeats, line noise, muscle, . . . ) [START_REF] Jung | Removing electroencephalographic artifacts by blind source separation[END_REF], making it an algorithm of choice for artifact rejection [START_REF] Urigüen | EEG artifact removalstate-of-theart and guidelines[END_REF]. ICA is widely used for the same purpose in MEG studies [START_REF] Mantini | Improving MEG source localizations: an automated method for complete artifact removal based on independent component analysis[END_REF][START_REF] Vigário | Independent component approach to the analysis of EEG and MEG recordings[END_REF][START_REF] Ikeda | Independent component analysis for noisy dataMEG data analysis[END_REF][START_REF] Dammers | Integration of amplitude and phase statistics for complete artifact removal in independent components of neuromagnetic recordings[END_REF].

Beyond artifact removal, ICA is also used to reveal and study brain activity.

In [START_REF] Makeig | Mining event-related brain dynamics[END_REF], ICA is successfully applied to recover evoked and induced event-related dynamics in EEG signals. In [START_REF] Gómez-Herrero | Measuring directional coupling between EEG sources[END_REF], ICA is used to extract brain sources, on which causal relations are exhibited, uncovering directional coupling. In [START_REF] Subasi | EEG signal classification using PCA, ICA, LDA and support vector machines[END_REF],

independent EEG sources are used in a machine learning pipeline, predicting epileptic seizures. ICA is used on MEG signals to identify links between function and structure in the brain in [START_REF] Stephen | Using joint ICA to link function and structure using MEG and DTI in schizophrenia[END_REF]. It can also be used on MEG data, coupled with Hilbert filtering, to uncover resting state networks [START_REF] Brookes | Investigating the electrophysiological basis of resting state networks using magnetoencephalography[END_REF].

Finally, ICA components can be mapped to certain brain areas via source localization. Indeed, the individual contribution of each source to each sensor can be represented as a topography on the scalp for EEG or on the helmet for MEG.

An equivalent current dipole (ECD) can then be fitted to the topography [START_REF] Scherg | Two bilateral sources of the late AEP as identified by a spatio-temporal dipole model[END_REF], yielding at the same time an estimate of the source location, and its dipolarity (how close it can be explained by a focal activity in the brain modeled with a single dipole).

The hypothesis of independence of the sources is at the heart of ICA. However, independence is a statistical property which is difficult to quantify on real data. In neuroscience, the most widely used algorithms are Infomax [START_REF] Bell | An information-maximization approach to blind separation and blind deconvolution[END_REF] and FastICA [START_REF] Hyvärinen | Fast and robust fixed-point algorithms for independent component analysis[END_REF].

These algorithms perform non-Gaussian ICA: they quantify independence on the marginal (instantaneous) distribution of the data. They ignore any time correlation and focus entirely on the non-Gaussianity of the data. In this case, the sources can be recovered when at most one source has a Gaussian density [START_REF] Comon | Independent component analysis, a new concept?[END_REF].

Brain sources and artifacts are usually heavy-tail signals which depart form

Gaussianity, justifying the use of non-Gaussian algorithms for M/EEG processing.

Another route to ICA is to leverage the time correlation of the sources. In this case, the sources can be recovered when the sources are spectrally diverse, that is, when their power spectrum are non-proportional [START_REF] Pham | Blind separation of mixture of independent sources through a quasi-maximum likelihood approach[END_REF]. Among these algorithms, Second Order Blind Identification (SOBI) [START_REF] Belouchrani | A blind source separation technique using second-order statistics[END_REF] is one of the most widely used. It jointly diagonalizes a set of time correlation matrices. Another approach closely related to our work consists in the joint diagonalization of spectral covariance matrices [START_REF] Pham | Blind separation of mixture of independent sources through a quasi-maximum likelihood approach[END_REF]. ICA methods based on joint-diagonalization of second order statistics might be less popular than non-Gaussian ICA methods, but have encountered some success in M/EEG processing. Congedo et al. [START_REF] Congedo | On the blind source separation of human electroencephalogram by approximate joint diagonalization of second order statistics[END_REF] argues that Pham's approach [START_REF] Pham | Blind separation of mixture of independent sources through a quasi-maximum likelihood approach[END_REF] should be prefered to SOBI for M/EEG preprocessing, one reason being that Pham's approach does not enforce orthogonal constraints.

Finally, in order to leverage both non-Gaussianity and spectral diversity, several methods based on short-time Fourier transform (STFT) have been proposed. For instance, Fourier-ICA [START_REF] Hyvärinen | Independent component analysis of short-time Fourier transforms for spontaneous EEG/MEG analysis[END_REF] leverages both non-Gaussianity and spectral diversity with a hybrid method, consisting of the non-Gaussian ICA of concatenated short-time Fourier transforms. Other approaches consist in joint diagonalization of cospectral matrices (covariance matrices of STFT frames) [START_REF] Congedo | On the blind source separation of human electroencephalogram by approximate joint diagonalization of second order statistics[END_REF] or work in the wavelet domain [START_REF] Pham | Source adaptive blind source separation: Gaussian models and sparsity[END_REF]. While all these algorithms rely on various independence measures, they make the strong assumption that there is no sensor noise: they assume that the signal of each sensor is a linear and noiseless combination of sources. A consequence of the noiseless model is that it enforces that there are as many sources as sensors, while the number of sensors is generally fixed by hardware constraints and not by the actual number of brain or artifactual sources present in the data.

Unfortunately these noiseless models lead to a degenerate likelihood when there are fewer sources than sensors. This is why, when fewer sources than sensors are expected to be present in the data, a dimension reduction technique like Principal Component Analysis (PCA) is often applied before ICA. However, this two-stage approach, consisting of first applying PCA and then ICA, is heuristic as based on the assumption that independent sources have high variance, which is not necessary. As it is argued in [START_REF] Artoni | Applying dimension reduction to EEG data by principal component analysis reduces the quality of its subsequent independent component decomposition[END_REF] applying PCA before ICA can degrade the quality of the recovered sources. To avoid relying on PCA from dimensionality reduction, it is also sometimes suggested to simply discard some channels. Throwing away data without a clear motivation is arguably questionable.

In order to alleviate this problem, some ICA algorithms incorporate a noise model. As explained in [START_REF] Hyvärinen | Independent component analysis in the presence of Gaussian noise by maximizing joint likelihood[END_REF], when the noise statistics are known, maximizing the likelihood of such a model is an optimization problem sharing many similarities with dictionary learning [START_REF] Olshausen | Emergence of simple-cell receptive field properties by learning a sparse code for natural images[END_REF]. Such procedure is typically much more costly than regular ICA, and it is seldom used in M/EEG processing (See [START_REF] Barthélemy | Multivariate temporal dictionary learning for EEG[END_REF][START_REF] Liu | A sparse dictionary learning framework to discover discriminative source activations in eeg brain mapping[END_REF] for instance). In [START_REF] Parra | Convolutive blind separation of non-stationary sources[END_REF], noise is modelled in a non-stationary framework: the sources -and noise-are assumed non-stationary (their instantaneous variance varies over time). The model is estimated by fitting it to the data, via the minimization of a simple quadratic criterion, which deviates from the probabilistic model.

In this article, we study Spectral Matching ICA (SMICA) for M/EEG processing. This ICA model has been first investigated in astronomy for separation of the cosmic microwave background [START_REF] Cardoso | Blind separation of noisy Gaussian stationary sources. application to cosmic microwave background imaging[END_REF][START_REF] Delabrouille | Multidetector multicomponent spectral matching and applications for cosmic microwave background data analysis[END_REF]. SMICA models the observations as a sum of a linear mixture of independent sources and noise. It assumes that the sources and noise are Gaussian, and that the sources have non-proportional spectra. This assumption makes it well suited for brain rhythms and artifacts extraction as they are known to have prototypical spectra. Brain sources tend to exhibit so-called "1/f" power spectral densities, while the spectra of artifacts are often localized in certain frequency bands (e.g. muscle artifacts or line noise). Importantly, the statistics of the noise are parameters of the model, and are estimated along the other parameters of the model. Thanks to its noise model, SMICA can estimate fewer sources than sensors without preprocessing for dimension reduction. The sources can be estimated by Wiener filtering, which takes the noise estimation into account and denoises the sources.

The article is organized as follows. In section 2, the SMICA statistical model is introduced and the estimation strategy based on an Expectation-Maximisation (EM) algorithm is described. In section 3, the usefulness of SMICA is demonstrated on various MEG and EEG datasets.

Notation The trace of a matrix M ∈ R p×p is Tr(M ), and its determinant is |M |. A matrix is invertible when |M | = 0, and we write M ∈ GL p . Given a vector u ∈ R p , the matrix diag(u) ∈ R p×p is the matrix containing the elements of u on its diagonal, and 0 elsewhere. If M is a p × p matrix, then diag(M ) is the diagonal matrix with the same diagonal as M . Given A ∈ R p×q , the vectorization of A is a vector vec(A) ∈ R pq of entries vec(A) i+p(j-1) = A ij . The Moore-Penrose Pseudo-Inverse of a tall matrix A ∈ R p×q is A † = (A A) -1 A .

A maximum likelihood approach to noisy ICA

This section introduces our approach to blind source separation for noisy observations (Sec. 2.1). Its application is then discussed in detail (Sec. 2.3).

The SMICA method in theory

In a noisy ICA model, the outputs of p sensors, e.g. M/EEG recordings, collected in a vector X(t) ∈ R p , are modelled as noisy instantaneous mixtures of q independent sources represented by a vector S(t) of size q with an additive noise term N (t) of size p, this is

X(t) = AS(t) + N (t), ( 1
)
where A is the p × q mixing matrix. The noise is assumed independent from the sources and uncorrelated across sensors.

This model readily translates into the spectral domain. Recall that for a zero-mean p-dimensional stationary time series, {X(t)}, the p × p autocovariance matrix E[X(t)X(t + τ ) ] does not depend on t and that its Fourier transform 1 :

C(f ) = τ E[X(t)X(t + τ ) ] e -2iπf τ (2)
defines p × p spectral covariance matrices C(f ). The diagonal entry C aa (f ) is the power spectrum of {X a (t)} while C ab (f ) contains the cross-spectrum between {X a (t)} and {X b (t)}.

The linear relation between data and sources of Eq. ( 1) translates into the spectral model

C(f ) = AP (f )A + Σ(f ) (3)
where P (f ) and Σ(f ) are the spectral covariance matrices of sources and of the noise. In this work, we assume that the sources and the noise terms are independent, which means that P (f ) and Σ(f ) are diagonal matrices. That property is the key feature enabling the blind estimation of the model since it

1 For simplicity we have set the sampling period to one time unit.

is how statistical independence between sources is expressed. Incidentally, it also entails that the spectral covariance matrices C(f ) are real-valued. Indeed, matrices P (f ) and Σ(f ) are real-valued since they are diagonal matrices and their diagonals are made of power spectra and the mixing matrix A also is real-valued. It follows from Eq. (3) that, unlike generic spectral matrices, a matrix C(f ) only has real entries, as a consequence of our statistical model.

This particular structure of the spectral covariance matrices is preserved when spectra are averaged over frequency bands. Define B frequency intervals

I 1 , . . . , I B by I b = [f b min , f b max ]
and consider frequency averages over those bands:

C b = 1 f b max -f b min f b max f b min C X (f )df (4)
Then, upon averaging, Eq. ( 3) becomes

C b = AP b A + Σ b (5)
where P b and Σ b denote the corresponding averages for P (f ) and Σ(f ). As a consequence, the noisy ICA model in Eq. ( 1) is transformed in the simpler model of Eq. ( 5), where the parameters are the mixing matrix A, the source powers in each band P b , and the noise powers in each band, Σ b .

The noisy ICA model is inferred from by connecting the spectral matrices C b of model ( 5) to samples estimates. If T data samples X(0), . . . , X(T -1) are avialble, spectral matrices are classically estimated from the Fourier coefficients

xk = 1 √ T T -1 t=0 X(t) e -2iπkt/T (6)
by averaging over the relevant frequency bands. These estimates are:

C b = 1 n b k: k T ∈I b Re xk xH k ( 7
)
where n b = #{k : by C b . This spectral matching principle is illustrated by Figure 1. We advocate using a specific spectral matching criterion:

L(θ) = B b=1 2n b KL C b , AP b A + Σ b ( 8
)
where KL is the Kullback-Leibler divergence between two p × p positive matrices:

KL(C 1 , C 2 ) = 1 2 Tr(C 1 C -1 2 ) -log det(C 1 C -1 2 ) -p . (9)
The KL-divergence KL(C 1 , C 2 ) is non-negative and cancels if and only if

C 1 = C 2 .
The particular measure [START_REF] Belouchrani | A blind source separation technique using second-order statistics[END_REF] of spectral adjustment between data and model has been chosen because it is (up to an irrelevant constant) asymptotically (for large T and narrow bands) equal to minus the log likelihood of a Gaussian stationary model. Hence, the SMICA estimates inherits some good properties of maximum likelihood estimates, in particular they enjoy a built-in scale invariance and they can be easily computed using the EM algorithm. Those properties and other considerations are discussed in the remaining of this section. The derivation of the SMICA criterion from a likelihood function is explained in Appendix A.

Statistical properties of SMICA

The novelty in SMICA comes from its noise model: ICA methods that work in the spectral domain are common. In fact, one recovers a popular method by removing the noise in Eq. ( 5), and by assuming a square mixing matrix p = q, that is, finding as many sources as sensors. This was considered by Pham who used the Whittle approximation in his seminal paper [START_REF] Pham | Blind separation of mixture of independent sources through a quasi-maximum likelihood approach[END_REF]. In that case (p = q, N (f ) = 0), the spectral mismatch can be shown to reduce to a joint-diagonality criterion. Estimating A amounts to solving the joint-diagonalization of the spectral covariances: A is such that the matrices A -1 C b A -are as diagonal as possible. The use of such approach for EEG processing is advocated in [START_REF] Congedo | On the blind source separation of human electroencephalogram by approximate joint diagonalization of second order statistics[END_REF]. In the experimental section 3, we compare SMICA to this approach of a plain joint diagonalization of spectral matrices and we refer to it as JDIAG.

The noise model in SMICA makes the estimation of the model parameters harder than joint-diagonalization, but it enables finer source estimation, through

S b = (A Σ -1 b A + P -1 b ) -1 A Σ -1 b X b , (10)
where X b is the matrix of signals filtered in the frequency band [f b , f b+1 ]. This operation is linear in X, and is adaptive to the level of noise: if in frequency band b the estimated noise Σ b for a sensor is large, then its contribution in the source estimate shrinks towards 0. Note that the standard ICA source estimation formula is recovered when the noise is equal on all channels and tends to 0 (i.e.

Σ b = λ b I p with λ b → 0), yielding at the limit S = A † X.

Next, we study the identifiability conditions for SMICA. The issue of identifiability is to find the conditions allowing for a unique recovery of the sources.

Identifiability conditions of noiseless spectral ICA are well established [START_REF] Pham | Blind separation of mixture of independent sources through a quasi-maximum likelihood approach[END_REF]: the sources should not have proportional spectra, that is, for any pair (i, j) of sources the quantities [P b ] ii /[P b ] jj should not be constant with respect to b = 1, . . . , B.

The identifiability conditions for the SMICA model ( 1) are more complicated, and the non-proportionality of the source spectra is only a necessary condition.

Hence, it is well suited to isolate all kinds of artifacts from brain activity: line noise, muscular activity, heartbeats, eye blinks, since their spectra are usually very different from those of brain activity. Among brain sources, some might have similar power spectrum, such as occipital dominant alpha rhythms and therefore cannot be accurately separated by SMICA.

Finally, we discuss the algorithm used to estimate the parameters, that is the specific numerical procedure to minimize the spectral matching criterion (or maximize the likelihood function) of Eq. ( 8). Since this is a likelihood function with latent (or unobserved) variables -namely the source and noise signalsit can be optimized by resorting to the celebrated Expectation-Maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF].

The EM algorithm is appealing because it does not require any hyperparameters like learning rates, and is guaranteed to decrease the loss function at each iteration. Still, this approach is generally slow (it might require many iterations to reach a satisfactory set of parameters), and other optimization techniques could be investigated for the fast minimization of the negative log-likelihood L. The EM algorithm for SMICA is described in the appendix (section Appendix C).

Next, we discuss the practical application of SMICA for M/EEG processing.

SMICA for M/EEG processing

This section goes through various considerations regarding the application of SMICA to M/EEG signals. A first advantage of SMICA over other ICA algorithms for M/EEG processing is its embedded dimension reduction. It is often the case that there are more sensors than sources which can be significantly recovered from the observations. When p > q matrix A cannot be inverted and one usually performs some dimension reduction in order to apply ICA algorithms which require a square matrix (p = q) to operate. In SMICA, thanks to the presence of a noise term, spectral covariance matrices are not degenerate even if q < p and a non-square matrix A can be fitted by SMICA in a statistically sound way even when there are fewer sources than sensors, without resorting to a pre-processing stage of dimension reduction. In some sense, one can say that SMICA has a built-in dimension reduction of the signal part because it can fit a tall (q < p) mixing matrix. It often argued that reduction dimension in EEG processing deteriorates the quality of the subsequent ICA decomposition [START_REF] Artoni | Applying dimension reduction to EEG data by principal component analysis reduces the quality of its subsequent independent component decomposition[END_REF].

SMICA offers a simple way to circumvent this problem, by embedding dimension reduction in the ICA in a statistically sound way.

Signal denoising using SMICA. Like any ICA algorithm, SMICA estimates sources which can be marked as spurious / non-biological by specialists: in addition to brain sources, it usually recovers physiological artifacts (heartbeats, eye blinks) and external electromagnetic perturbations (room and line noise).

The remaining sources can then be projected back in the signal space, giving clean M/EEG signals. Thanks to its noise modeling, SMICA makes these two operations statistically sound: the sources S are estimated by Wiener filtering, some sources can be manually or automatically marked as spurious. If the source i is marked as spurious, we set S i = 0. The cleaned M/EEG signals are computed as X cleaned = A S.

Combining SMICA with another ICA algorithm. In practice, it might happen that some sources recovered by SMICA have similar spectra, which indicates that these sources are not well separated. Another ICA procedure (typically, based on non-Gaussianity) can then be applied on these sources in order to better separate them.

Taking it a step further, SMICA can also be used as a source subspace identifier. Applied on the M/EEG signals X ∈ R p×T , SMICA produces a mixing matrix A 1 ∈ R p×q and a source matrix S 1 ∈ R q×T estimated by Wiener filtering.

The sources in S 1 are maximally independent with respect to the separation criterion of SMICA, but might not be maximally independent with respect to another criterion. Applying another ICA algorithm on S 1 yields a square mixing matrix A 2 ∈ R q×q , and a new source matrix S 2 = A -1 2 S 1 . The overall mixing matrix linking the sources S 2 and the original dataset X is the matrix product A T otal = A 1 A 2 . For instance, using a non-Gaussian ICA algorithm like Infomax or FastICA on the sources found by SMICA may disentangle sources that share similar power spectrum. The practical benefits of such approach are demonstrated in the experiments presented in section 3. In practice, we found that fitting SMICA on a 102 sensors MEG dataset with 40 frequency bins and 100 sources takes about 15 minutes using one CPU of a recent laptop.

Complexity

Frequency selection. SMICA exploits spectral information, but not necesarily over the whole frequency range. Typically, one may exclude the highest frequencies if they are dominated by noise or cut off by a sampling filter. One may also ignore the very lowest frequencies if they are dominated by slow drift artifacts.

In general, there is no counter indication to restricting the sum (8) to the high SNR part of the frequency range.

Experiments

We report some experiments comparing our approach to other (noiseless) ICA algorithms: these algorithms all model the dataset X as X = AS, where A ∈ R p×p is the square mixing matrix, and S is the source matrix. However, they estimate the independent components based on different independence criterion.

We start by briefly describing the ICA algorithms which are compared to SMICA.

Non-Gaussian ICA. Non-Gaussian ICA algorithms model the source time-series as independent and identically distributed, with non-Gaussian probability density functions. Some of the most popular ICA algorithms fall in this category:

FastICA [START_REF] Hyvärinen | Fast and robust fixed-point algorithms for independent component analysis[END_REF], Infomax [START_REF] Bell | An information-maximization approach to blind separation and blind deconvolution[END_REF], its extended version [START_REF] Lee | Independent component analysis using an extended infomax algorithm for mixed subgaussian and supergaussian sources[END_REF], JADE [START_REF] Cardoso | Blind beamforming for non-gaussian signals[END_REF] and more recently AMICA [START_REF] Palmer | AMICA: An adaptive mixture of independent component analyzers with shared components[END_REF].

Second order blind identification. The SOBI algorithm [START_REF] Belouchrani | A blind source separation technique using second-order statistics[END_REF] aims at recovering sources with spectral diversity, just like SMICA and JDIAG. It does so in a heuristic way, by joint-diagonalization of a set of correlation matrices

1 T -τ T -τ
t=1 X(t)X(t -τ ) for a set of time lags τ 1 , . . . , τ B , rather than spectral covariance matrices. Choosing an appropriate set of time lags is not an obvious task; we use the set advised in [START_REF] Tang | Validation of SOBI components from high-density EEG[END_REF]. Unlike JDIAG, the joint diagonalization criterion is ad-hoc, and does not correspond to a principled statistical criterion.

For instance, since JDIAG follows the maximum-likelihood principle, it is asymptotically Fisher-efficient and reaches the Cramer-Rao lower bound, unlike SOBI. This is why several articles argue for the use of JDIAG rather than SOBI [START_REF] Doron | Asymptotically optimal blind separation of parametric Gaussian sources[END_REF][START_REF] Congedo | On the blind source separation of human electroencephalogram by approximate joint diagonalization of second order statistics[END_REF].

Estimating fewer sources than sensors using PCA. Contrarily to SMICA, algorithms described above can only estimate as many sources as sensors. Therefore, in order to estimate fewer sources, a dimension reduction step should be performed beforehand. Principal Component Analysis is the algorithm of choice for this task. The components are chosen to explain as much variance in the data as possible. Since this method is blind to higher order interactions, it might discard some sources which are important but of low power. As a consequence, Artoni et al. [START_REF] Artoni | Applying dimension reduction to EEG data by principal component analysis reduces the quality of its subsequent independent component decomposition[END_REF] argue that applying PCA before ICA leads to degraded decomposition.

Numerical setup. In our experiments, we take Infomax as the reference non-Gaussian ICA algorithm with the fast and robust optimization algorithm Picard [2, 1] and using tanh(•) as the non-linear activation function. The jointdiagonalization algorithm for SOBI is a combination of [START_REF] Ziehe | A fast algorithm for joint diagonalization with non-orthogonal transformations and its application to blind source separation[END_REF] with a backtracking line-search. The joint-diagonalization algorithm for JDIAG uses the fast implementation described in [START_REF] Ablin | Beyond Pham's algorithm for joint diagonalization[END_REF]. In all experiments we set the frequency bins F b of SMICA and JDIAG as uniform in the range 1 -70 Hz, with 40 bins. The M/EEG analysis is conducted using the Python package MNE [START_REF] Gramfort | MNE software for processing MEG and EEG data[END_REF][START_REF] Gramfort | MEG and EEG data analysis with MNE-Python[END_REF]. Figures are made using Matplotlib [START_REF] Hunter | Matplotlib: A 2d graphics environment[END_REF].

The python code for SMICA is available online at https://github.com/ pierreablin/smica.

Qualitative comparisons

Comparisons on a MEG dataset

We start by showing the decomposition found by SMICA, JDIAG, SOBI and Infomax on a MEG dataset, where the subject was presented checkerboard patterns into the left and right visual fields, as well as monaural auditory tones to the left or right ears. Stimuli occurred every 750 ms (See Gramfort et al. [START_REF] Gramfort | MNE software for processing MEG and EEG data[END_REF] for a description of the dataset). MEG is acquired with 102 magnetometers and 204 gradiometers.

For this experiment, we only consider the 102 magnetometer channels. Each ICA algorithm returns 40 sources (after PCA for JDIAG, SOBI and Infomax).

We hand pick 10 sources to display in Figure 2, which shows their time-course, power spectrum and topography. SMICA isolates heartbeats (source 1), eye blinks (source 7) and brain (sources 8-10) from line noise: only source 9 shows a very small peak at 60 Hz. Source 7 and 9 in JDIAG's decomposition, source 1, 7, 9, 10 for SOBI and source 1, 8, 9 for Infomax do show a higher peak at 60 Hz suggesting that line noise leaks 1 Source Spectrum Topo into these others sources. SOBI fails to separate properly the eye-blinks, and gives decompositions that differ substantially from other algorithms. This first experiment demonstrates that SMICA does reveal expected artifactual sources, both physiological and environmental, and that they potentially leak less into the valuable neural ones. More quantitative evidence is provided below.

Comparisons on an EEG dataset

We run SMICA and JDIAG on a 69-channel EEG data coming from the dataset described in [START_REF] Delorme | Independent EEG sources are dipolar[END_REF]. Both algorithms return 20 sources, which are displayed in Figure 3. Differences between SMICA and JDIAG are now more striking, probably due to the greater noise level compared to the MEG recording.

Although decompositions differ in many aspects, we want to focus on source 2 recovered by SMICA, which is not found by JDIAG. There is a sharp peak at 60 Hz, indicating that it likely corresponds to line noise. The second peak at 50

Hz may seem spurious but is most probably due to spectral aliasing. Indeed the fifth harmonic of a line at 60 Hz when sampled at 250 Hz appears at frequency 5 * 60 -250 = 50 Hz. To test whether it is a plausible line source, we resort to a separate experiment: we determine the spatial filter w ∈ R 59 such that the time series wX is of power 1 with a maximal power at 60 Hz. This method is called Spatio-Spectral Decomposition (SSD) [START_REF] Nikulin | A novel method for reliable and fast extraction of neuronal EEG/MEG oscillations on the basis of spatio-spectral decomposition[END_REF]. To do so, we filter X in a narrow band around 60 Hz, yielding signals X f . Then, we find w by maximizing the power of wX f under the constraint that the power wX is 1 which is done by maximizing the Rayleigh quotient between the covariance of X f and of X.

The power spectrum of the corresponding source wX is displayed in Figure 3, along with the power spectrum of the source number 2 found by SMICA. The 50 Hz aliased harmonic is also recovered by SSD, suggesting that the source recovered by SMICA isolates the line signal correctly.

Quantitative results on large datasets

Experiment on MEG Phantom data

In this section, we experiment with MEG Phantom data. The recording comes from a fake plastic head with electric dipoles. Dipoles emit sinusoidal pulses at a fixed frequency 20 Hz for 0.5 second (10 periods) and are then turned off for 1 second. This is repeated for 150 seconds.

We have 24 datasets. Each dataset corresponds to one dipole location among 8 and one amplitude (either 1000, 200 or 20 nAm). The amplitude corresponds to the peak-to-peak difference. We cut each dataset in half to obtain twice as many datasets.

As the true locations of the sources in the phantom are known one can map the origin of an ICA source by fitting an equivalent current dipole (ECD) to the source topographies and evaluate the localization errors. On each dataset, we accurately recover the eye blinks (source 1). For SMICA, sources 4-20 correspond to brain activity. SMICA finds two dipolar beta-rhythm sources (sources 7 and 10). Each algorithm recovers a source corresponding to line noise, with a large peak at 60 Hz (source 2). For SMICA, there is an additional peak at 50 Hz, which is not an artifact. Spatio-Spectral Decomposition (SSD) tuned to the source of maximal power at 60 Hz yields a similar peak at 50 Hz (c).

apply SMICA, Infomax and JDIAG to obtain 40 sources. For each source, we fit an ECD, and only keep the source corresponding to the closest location to the true dipole.

Besides ICA, Maxwell filtering can also be used for dipole localization. After Maxwell filtering, we compute the evoked potential and fit a dipole at the peak.

Finally, we employ Spatio-Spectral Decomposition (SSD), as described in section 3, to find a linear combination of sensors with maximal amplitude at 20 Hz. This method incorporates more knowledge of the problem than others, because we provide it with the dipole frequency.

In Figure 4, we display the average distance to the true dipole and the residual variance in the dipole fit, for each dipole amplitude (1000, 200 and 20 nAm). The localization errors increase as the sources amplitude diminishes for each algorithm. Yet, for the smallest amplitude (20 nAm), which is the most challenging, SMICA outperforms all other methods in terms of localization (note the logarithmic scale).

SMICA finds highly dipolar and independent source subspaces

In this section, we illustrate the ability of SMICA to capture meaningful brain sources. We use the same datasets as in [START_REF] Delorme | Independent EEG sources are dipolar[END_REF]. It contains the EEG recording of 15 subjects, with 69 EEG channels. For a target number of independent sources, different ICA procedures described in the article are applied to the datasets. The Wiener and Pseudo-Inverse methods correspond to the combination of SMICA with Infomax, as described in the paragraph 'Combining SMICA with another ICA algorithm' of section 2.3. Wiener corresponds to computing SMICA's sources with Wiener filtering, Pseudo-Inverse corresponds to computing SMICA's sources with pseudo-inversion of the mixing matrix A, as described in the paragraph 'Source estimation by Wiener filtering' of Sec. 2.3. For each decomposition, we compute the dipolarity of each source, as well as the pairwise mutual-information between each pair of sources.

Results for the dipolarity are displayed in Figure 5, results for the pairwise mutual information are displayed in Figure 6. The figure should be understood in the following way. Looking at the first plot, corresponding to 10 sources, we see that about 80% of components found by the method 'Wiener' have an equivalent dipole map residual variance lower than 10%. About 65% of components found by 'Infomax' have an equivalent dipole map residual variance lower than 10%. Overall, the method 'Wiener' finds more dipolar components, for every number of sources considered. The results of raw SMICA are not displayed here, as it recovered less dipolar sources than Infomax. That is not unexpected if one follows the discussion of [START_REF] Delorme | Independent EEG sources are dipolar[END_REF], arguing that maximum-likelihood non-Gaussian methods, like Infomax, are the ICA methods that recover the most dipolar sources.

For this experiment, using SMICA as a subspace identifier to perform dimension reduction yields the best results, regardless of the number of sources that are recovered. The benefits of the Wiener filtering over pseudo-inversion are also illustrated: even without dimension reduction, it denoises the signals, which leads to improved ICA decompositions.

Discussion

We now discuss the advantages and shortcomings of SMICA for M/EEG processing.

Noise model. Results in Sec. 3.1.1 (Figure 2) have shown that SMICA can better isolate stationary artifacts (environmental line noise, muscule activity, eye blinks, heartbeats,...) than non-Gaussian ICA methods. This can be explained by the fact that artifacts and brain sources usually have very different spectra. Artifacts tends to have peaks or bumps in the spectra in higher frequencies while neural sources have spectra with power laws or exponential decays [START_REF] Buzsáki | Neuronal oscillations in cortical networks[END_REF][START_REF] Dehghani | Comparative power spectral analysis of simultaneous elecroencephalographic and magnetoencephalographic recordings in humans suggests non-resistive extracellular media[END_REF][START_REF] Rocca | Self-similarity and multifractality in human brain activity: A wavelet-based analysis of scale-free brain dynamics[END_REF]. The rate of decay of the different spectrum of neural sources also depends on the underlying source. This spectral diversity can be exploited by SMICA to improve source separation. Overall, this method is well suited to separate artifacts from neural activity. We would also like to stress that these artifacts are captured as sources by SMICA: they correspond to "S" in the equation X = AS + N , and not to the modelling noise "N ". For instance, if an external signal, independent from other sources, has a correlated spectral power over the sensors, then SMICA will capture it as a source. By assumption, the modelled noise N only identifies sensor noise: noise that is decorrelated across sensors. Because we estimate the spectral power of noise, the model is flexible enough to capture any such signal.

SMICA as a dimension reduction tool. Thanks to the noise model, we have a principled way to perform dimension reduction and to recover the source time courses. In [START_REF] Artoni | Applying dimension reduction to EEG data by principal component analysis reduces the quality of its subsequent independent component decomposition[END_REF], it is argued that PCA is suboptimal for EEG data, and that even channel subsampling is to be preferred. By contrast, the subspace identified by SMICA contains sources that are estimated by incorporating spectral information of the original signals, and these sources have spectral diversity. This benefit is observed in practice. In Sec. 3.3 (Figure 5 and 6) we show that SMICA identifies a source subspace that contains more dipolar and independent components than PCA. As such, SMICA can be a useful tool for dimension reduction.

Comparison with noiseless model. Importantly, as illustrated by Figure 2 Other spectral approaches. This work is not the first one using spectral ICA and highlighting its benefits for M/EEG recordings. Some works have focused on convolutional mixtures (for which using the Fourier domain turns convolutions into products). In [START_REF] Dyrholm | Model selection for convolutive ICA with an application to spatiotemporal analysis of EEG[END_REF], authors use convolutive ICA in time-domain, while [START_REF] Anemller | Complex independent component analysis of frequency-domain electroencephalographic data[END_REF] use a complex Infomax to find travelling waves. In the end, these methods estimate one mixing matrix for each frequency bin. However, independent models in each band might fail to recover brain rhythms with several frequency peaks. For instance, mu-rhythm is characterized by concurrent activity near 10 Hz and 20 Hz [START_REF] Hari | Human cortical oscillations: a neuromagnetic view through the skull[END_REF][START_REF] Niedermeyer | Electroencephalography : basic principles, clinical applications, and related fields[END_REF].

Automatic estimation of the number of sources?. Some improvements to the current SMICA algorithm can be investigated. First, since it comes from a principled statistical framework, it would be interesting to implement a data-driven way of computing the most likely number of sources in the data: an algorithm to automatically select the correct number of sources. However, preliminary experiments show that usual statistical criteria like Akaike Information Criterion or Bayesian Information Criterion are not satisfactory in this setting, likely because the model is not complex enough to explain fully M/EEG signals. The EM algorithm for fitting SMICA is also quite slow, some improvements could be possible by further studying the geometry of the cost function and proposing quasi-Newton algorithms, as done recently for Infomax [START_REF] Ablin | Faster independent component analysis by preconditioning with Hessian approximations[END_REF].

Possible extensions of SMICA. SMICA could be extended in several interesting ways. In MEG acquisition, the empty room is sometimes recorded before the experiment. In this case, we could assume that the noise term in Eq. In order to obtain a pure subspace identification method, one could also drop the hypothesis of independence between sources and minimize the spectral matching criterion with respect to any positive spectral covariance matrices P b , instead of constraining them to be diagonal. 2 .

Finally, the proposed method could also be extended to non-stationary signals,

where the spectral covariance matrices are replaced by time-lagged covariance matrices. The resulting algorithm would resemble the procedure proposed in [START_REF] Parra | Convolutive blind separation of non-stationary sources[END_REF],

but with a proper likelihood-based estimation rather than the ad-hoc criterion proposed by the authors. More generally, the spectral model with noise and EM algorithm can be employed to recover the parameters of a noisy ICA model X = AS + N using any kind of second-order statistics.

Conclusion

In

  k T ∈ I b } denotes the number of Fourier coefficients available in band b. The set θ = (A, P 1 , . . . , P B , Σ 1 , . . . , Σ B ) of all unknown parameters can be estimated by adjusting the model C b = AP b A + Σ b to the data as summarized

Figure 1 :

 1 Figure 1: The SMICA method: empirical spectral covariance matrices Ĉ1 , . . . , ĈB computed from the M/EEG data are matched by the model Ĉb AP b A + Σ b , where the P b and Σ b are diagonal positive matrices. Matching is performed in a statistically sound way, by using a matching criterion derived from a simple likelihood.

  Wiener filtering ([START_REF] Brown | Introduction to random signals and applied Kalman filtering[END_REF], Chapter 4). In noiseless ICA models of the form X = AS, the natural source estimates simply is S = A -1 X. SMICA can employ the same technique for recovering the source, albeit replacing A -1 by A † , the Moore-Penrose pseudo-inverse of A. However, the availability of a noise term in SMICA offers a more attractive possibility: one can compute the expected value of the sources given the parameters: S = E[S|X, θ] which yields the lowest noise contamination among all unbiased estimators. This is the well known Wiener filtering which has a simple explicit expression in the SMICA model: in each frequency band, it reads

  of SMICA. Estimation of the parameters of SMICA is a two step process: first, the empirical spectral covariances Ĉ1 , . . . , ĈB are computed, and then the EM algorithm is used to infer the parameters, based solely on those covariances. Only the first step scales with the length of the signal T : the second only depends on the dimensions of the problem p, d and the number of bins B. In practice, we find that in the setting of our experiments, computing the covariances takes a negligible time compared to the time it takes to infer the parameters with the EM algorithm. The complexity of the EM algorithm does not depend on the number of samples T ; only the computation of the covariances does. This differs from non-Gaussian ICA algorithms, for which the estimation time is roughly proportional to the length of the recordings.

Figure 2 :

 2 Figure 2: Different ICA decompositions on MEG data. Source 1 corresponds to heartbeats, sources 2-6 to environmental noise, with strong peaks around 60 Hz. Source 7 corresponds to eye blinks. Sources 8, 9 correspond to occipital alpha rhythm (excect for SOBI which did not find such sources). Source 10 corresponds to a dipolar auditory source. Sources 7, 9 found by JDIAG present artifactual 60 Hz components, like source 1, 9 for Infomax. SOBI does not properly identify the eye-blink source.

  Spectrum of source 2 of SMICA and source obtained by Spatio-Spectral Decomposition (SSD) at 60Hz.

Figure 3 :

 3 Figure 3: Comparison of SMICA and JDIAG on an EEG dataset. Both algorithms return 20 sources. The ordering of the sources of SMICA is made by hand; the ordering of the sources of JDIAG is made by maximizing the correlation with SMICA's sources. Both algorithms

Figure 4 :Figure 5 :

 45 Figure 4: Dipole localization on MEG Phantom dataset. Top: Distance between the true dipole and the estimated dipole. Bottom: Residual variance in the dipole fit. Each column corresponds to a different source amplitude (from left to right, 1000, 200 and 20 nAm). Dipole fitting is applied on 24 datasets. Each black dot corresponds to a dataset. When the source dipole amplitude is low (20 nAm), SMICA has better localization performance than the other methods.

Figure 6 :

 6 Figure 6: Pairwise mutual information (PMI) for different ICA decompositions returned by each algorithm on the 15 datasets of 69 sensors, where each algorithm returns 10, 20, 30, 40, 50 or 69 sources. For Infomax, SOBI and JDIAG, PCA is first applied on the data matrix to obtain the desired number of channels. Wiener and Pseudo-Inverse correspond to applying Infomax on the sources recovered by SMICA, either by the Wiener or Pseudo-Inverse method. PMI is displayed in a Q-Q plot showing the quantiles of the distribution of PMI found by each algorithm against the quantiles of distribution of PMI found by Infomax. Infomax therefore corresponds to the line x = y. Algorithms above the line x = y have less PMI remaining than Infomax, and hence obtain more independent sources than Infomax as quantified by PMI. Here again the the Wiener approach is particularly competitive, especially when few sources are estimated.

  (1) shares the same spectral signature as the empty room recording: the matrices Σ b are no longer estimated by the model, but are now taken as the spectral covariance matrices of the empty room. Only the mixing matrix A and the source powers P b are left to estimate. The sources estimated by the algorithm should then only correspond to biological sources; in particular they should automatically be cleaned from line noise.

= 2

 2 this work, we have introduced a novel ICA algorithm, SMICA, which adds a noise model to the standard ICA model. By assuming a model of Gaussian stationary sources, we obtained a tractable closed-form likelihood, which can then be maximized with the expectation-maximisation algorithm. Blind identifiability stems from spectral diversity: this likelihood permits the separation of sources with different power spectra. The model can estimate fewer sources than sensors, and the sources can then be recovered by Wiener filtering, which takes noise into account. We then demonstrated the promises of our method for M/EEG signal processing: compared to noise-free ICA algorithms, SMICA extracts sources that are more dipolar and have less pairwise mutual information. Besides when applied on controlled recordings using an MEG phantom, SMICA is able to locate dipoles more precisely, especially for low amplitude sources. These results indicate that SMICA is a promising alternative to other more standard ICA algorithms for M/EEG signal denoising and exploration. invariance on the structure of C and give the joint pdf of a and b when they are jointly Gaussian. We start with:C = E[xx H ] = E[(a + ib)(a + ib) H ] = E[(aa + bb ) -i E[ab -ba ],andE[xx ] = E[(a + ib)(a + ib) ] = E[(aa -bb ) + i E[ab + ba ]. If x is changed into e iφ x, then matrix E[xx H ] is unchanged but E[xx ] is changed into e 2iφ E[xx ]. However, by the phase invariance, there should be no change. That is only possible if E[xx ] = 0. Therefore E[aa -bb ] = 0 and E[ab + ba ] = 0. Combining that with the assumption of a real covariance matrix C = E[xx H ] which implies E[ab -ba ] = 0 yields E[aa ] = E[bb ] = C/2 E[ab ] = E[ba ] = 0 Therefore the joint (2p) × (2p) covariance matrix for the pair (a, b) If a and b are jointly Gaussian, their probability density p(a, b) is given by log p(a, b)The block structure of Σ yields log det(2πΣ) = 2 log det(πC) and also Tr C -1 (aa + bb ) .Since C is a real (by assumption) symmetric matrix, we haveTr C -1 (aa + bb ) = Tr C -1 xx H = Tr C -1 Re(xx H )Combining the previous two results yields the joint density:log p(a, b) = -Tr C -1 Re(xx H ) -log det(πC) (B.1)algorithm, we first fit the model with a fixed noise level for each bin: we estimate Σ subject to Σ b = Σ for all i. In this setting, the M-step is much simpler and computationally quicker. Then, the core SMICA algorithm with free noise starts with Σ b all equal to the estimated noise level, and A and P b start from the same initial value.

In this setting, it is possible, but not mandatory, to constrain matrix A to have orthonormal columns: (A A = Iq).

Acknowledgements

This project has received funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (Grant agreement SLAB ERC-StG-676943). P.A. acknowledges support from the French government under management of Agence Nationale de la Recherche as part of the Investissements d'avenir program, reference ANR19-P3IA-0001 (PRAIRIE 3IA Institute).

Competing interests statement

Jean-François Cardoso: Conceptualization, methodology, software, investigation, writing -review and editing, supervision.

Alexandre Gramfort: Conceptualization, methodology, software, writingreview and editing, supervision, funding acquisition.

Appendix A. Spectral mismatch and likelihood

We show the statistical origin of the measure (8) of spectral mismatch. The starting point is that, for a zero-mean p-variate stationary times series with spectral covariance matrix C(f ), the Fourier coefficients xk defined at Eq.( 6) have zero mean and covariance E[x k

xH k ] = C(k/T ). Next, asymptotically (for large T ), these coefficients also are normally distributed and pairwise uncorrelated [START_REF] Brillinger | Time series: data analysis and theory[END_REF]. The Whittle approximation to the likelihood consists in assuming that these properties hold even in finite sample size. In that case, the probability density for the Fourier coefficients (x 1 , . . . , xT/2 ) is the product over frequencies of Gaussian densities: log p(x 1 , . . . , xT/2 ) = T /2 k=1 log CN (x k ; C(k/T )) where CN (x; C) denotes the complex circular Gaussian density for a zero-mean vector x of covariance matrix

In Appendix B, we give a simple expression (Eq. (B.1)) for CN (x; C) in the special case where C(f ) is real-valued, yielding:

It is worth noting that, in general, the covariance matrix of a circular complex Gaussian vector is complex-valued. However, the noisy stationary ICA model considered in this paper does enjoy real-valued spectral covariance matrices, as seen from Eq. ( 3).

The final step is to approximate the spectra as constant over thin spectral bands, that is,

Then, the sum (A.1) over frequencies can be expressed as a sum over spectral bands:

Hence, with the definition (9) of the KL divergence, we obtain

which shows that, up to a constant term, the spectral matching criterion is minus the logarithm of the Whittle likelihood (A.2).

The presence of the factor 2 in (A.2) can be traced back to the fact that a Fourier component at frequency k has 2 degrees of freedom: its real and its imaginary parts, or equivalently, that in real space any Fourier component has a sine and a cosine terms.

Appendix B. Likelihood and complex vectors

We give joint pdf for the real and imaginary parts of a zero-mean complex random p × 1 vector x = a + ib with a real covariance matrix:

and whose distribution is invariant under any phase change, that is, x has the same distribution as e iφ x for any angle φ. We look at the consequence of this

M-step. At the M-step, the parameters of the model θ should be modified in order to decrease the loss function, using the sufficient statistics obtained in the E-step.

To do so, we compute the EM functional, Φ(A, (P b ), (Σ b ), S) = -log(p(X, S|θ)).

To begin, we assume that there is only on frequency band (B = 1), which gives a Gaussian "white" model: X = AS + N , with S ∼ N (0, P ) and N ∼ N (0, Σ).

We find on the one hand

and on the other hand

Then, in the Whittle approximation, the EM functional for a spectral model is simply the weighted sum of the previous EM functionals in each band, which gives the EM functional for SMICA:

which should be minimized with respect to the parameters θ.

• Optimizing P b : the source powers are decoupled from the other parameters in (C.7). Minimization of Φ w.r.t. P b is easily obtained by canceling the gradient, yielding:

• Optimizing Σ b : the mixing matrix A and the noise covariance are entangled in eq. (C.7), rendering the analytic minimization of Φ impossible. Therefore, we first minimize Φ w.r.t Σ b , keeping A constant. This yields:

• Optimizing A: keeping the noise levels fixed, minimizing Φ w.r.t. A yields, by canceling the gradient:

This can be seen as a system of equations for the rows of A which, thanks to the diagonality of Σ b , is easily seen to decouple across the rows. For each row, simple algebra yields the close form solution: for r = 1 . . . p, A r: = Q r: M -1

Therefore, the EM update of A only requires solving p linear systems of size q × q.

Implementation details. The EM algorithm iterates the E and M step until a certain convergence criterion is reached. In practice, iterations are stopped when the difference between two consecutive values of the log-likelihood is below a threshold: L t+1 > L t -ε. In order to have a good initialization for the


domain: sdv info

text:
Baptiste Morel
email: bamorel@univ-tours.fr

Pierre Bertault

Géraldine Favrais

Elsa Tavernier

Barthelemy Tosello

Nathalie Bednarek

Laurent Barantin

Alexandra Chadie

Maia Proisy

Yongchao Xu

Isabelle Bloch

Dominique Sirinelli

Catherine Adamsbaum

Clovis Tauber

Elie Saliba

Automated brain MRI metrics in the EPIRMEX cohort of preterm newborns: Correlation with the neurodevelopmental outcome at 2 years

Keywords:

Premature;

Magnetic resonance imaging; Brain; cohort studies

List of abbreviations

ASQ: age and stage questionnaire; BGT: basal ganglia and thalami; DEHSI: diffuse excessive high signal intensity; GA: gestational age; MRI: magnetic resonance imaging; SD: Standard deviation; TEA: term equivalent age; WA: weeks of amenorrhea; WM: white matter

Introduction

Owing to progress in neonatal care, the preterm newborn survival rate has increased during recent years. However, prematurity remains a leading cause of morbidity and conveys high risk of developmental delay [START_REF] Saigal | An overview of mortality and sequelae of preterm birth from infancy to adulthood[END_REF]. Clinical perinatal factors related to adverse outcomes have been identified, such as lower gestational age (GA) [START_REF] Wolke | Universal gestational age effects on cognitive and basic mathematic processing: 2 cohorts in 2 countries[END_REF], intrauterine growth restriction [START_REF] Murray | Differential effect of intrauterine growth restriction on childhood neurodevelopment: a systematic review[END_REF], bronchopulmonary dysplasia [START_REF] Doyle | Long-term outcomes of bronchopulmonary dysplasia[END_REF], infection [START_REF] Van Vliet | Perinatal infections and neurodevelopmental outcome in very preterm and very low-birth-weight infants: a metaanalysis[END_REF], necrotizing enterocolitis, severe brain lesions on cranial ultrasound or brain magnetic resonance imaging (MRI) [START_REF] Hintz | Neuroimaging and neurodevelopmental outcome in extremely preterm infants[END_REF], but are not yet sufficient to predict accurately all disability developments. Early identification of newborns at risk of later developmental disabilities may increase access to specific evaluation and subsequent adapted care, potentially influencing the course of otherwise persistent difficulties.

Brain MRI is a very powerful complementary tool to diagnose and follow injury in preterm infants [START_REF] Maalouf | Comparison of findings on cranial ultrasound and magnetic resonance imaging in preterm infants[END_REF]. Periventricular leukomalacia associated with cystic lesions, white matter reduced volume, stigmata of hemorrhage and delayed myelination are possible brain sequelae of prematurity.

Simple brain metric such as biparietal width or interhemispheric distance did not reliably predict neurodevelopmental outcome in a cohort with low prevalence of high-grade brain injury [START_REF] Dewan | Are simple magnetic resonance imaging biomarkers predictive of neurodevelopmental outcome at two years in very preterm infants?[END_REF]. Associations of conventional brain MRI at term equivalent age (TEA) with long-term neurodevelopmental impairments in preterm children are of particular interests with high specificity [START_REF] Anderson | Associations of newborn brain magnetic resonance imaging with long-term neurodevelopmental impairments in very preterm children[END_REF][START_REF] Martinez-Biarge | Neurodevelopmental outcomes in preterm infants with white matter injury using a new MRI classification[END_REF]. Structural and volumetric MRI at TEA and neurological examinations could predict long-term neuromotor outcome in premature newborns [START_REF] Hintz | Neuroimaging and neurodevelopmental outcome in extremely preterm infants[END_REF][START_REF] Setänen | Prediction of neuromotor outcome in infants born preterm at 11 years of age using volumetric neonatal magnetic resonance imaging and neurological examinations[END_REF][START_REF] Jeon | Neurodevelopmental outcomes in preterm infants: comparison of infants with and without diffuse excessive high signal intensity on MR images at near-term-equivalent age[END_REF][START_REF] Mürner-Lavanchy | Thirteen-year outcomes in very preterm children associated with diffuse excessive high signal intensity on neonatal magnetic resonance imaging[END_REF]. But neonatal brain MRI segmentation is difficult to perform due to the imaging characteristics of the developing brain (smaller head size and shorter scanning duration).

Moreover, quantitative analysis of white matter in premature newborns remains a challenging issue [START_REF] Morel | A challenging issue: detection of white matter hyperintensities in neonatal brain MRI[END_REF]. Maalouf et al. have described diffuse excessive high signal intensity (DEHSI) of the white matter in the centrum semiovale on conventional T2-weighted images [START_REF] Maalouf | Magnetic resonance imaging of the brain in a cohort of extremely preterm infants[END_REF]. Increasing severity of white-matter abnormalities on MRI at term equivalent is associated with poorer performance on the cognitive and psychomotor scales of the BSID-II [START_REF] Woodward | Neonatal MRI to predict neurodevelopmental outcomes in preterm infants[END_REF]. When focusing on DEHSI, these results were supported by Dyet et al. [START_REF] Dyet | Natural history of brain lesions in extremely preterm infants studied with serial magnetic resonance imaging from birth and neurodevelopmental assessment[END_REF] and Parikh et al. [START_REF] Parikh | Automatically quantified diffuse excessive high signal intensity on MRI predicts cognitive development in preterm infants[END_REF][START_REF] Parikh | Objectively diagnosed diffuse white matter abnormality at term is an independent predictor of cognitive and language outcomes in infants born very preterm[END_REF]. The prognosis was secondarily controverted. Normal neurologic outcomes were observed at a corrected age of 18 months to 36 months in premature with DEHSI [START_REF] Jeon | Neurodevelopmental outcomes in preterm infants: comparison of infants with and without diffuse excessive high signal intensity on MR images at near-term-equivalent age[END_REF][START_REF] Calloni | Neurodevelopmental outcome at 36 months in very low birth weight premature infants with MR diffuse excessive high signal intensity (DEHSI) of cerebral white matter[END_REF][START_REF] De Bruïne | Clinical implications of MR imaging findings in the white matter in very preterm infants: a 2-year follow-up study[END_REF].

Today, the prognostic value of DEHSI is still controversial, ranging from a normal transient feature to a predictive factor for mental retardation [START_REF] Jeon | Neurodevelopmental outcomes in preterm infants: comparison of infants with and without diffuse excessive high signal intensity on MR images at near-term-equivalent age[END_REF][START_REF] Mürner-Lavanchy | Thirteen-year outcomes in very preterm children associated with diffuse excessive high signal intensity on neonatal magnetic resonance imaging[END_REF][START_REF] Dyet | Natural history of brain lesions in extremely preterm infants studied with serial magnetic resonance imaging from birth and neurodevelopmental assessment[END_REF][START_REF] Parikh | Automatically quantified diffuse excessive high signal intensity on MRI predicts cognitive development in preterm infants[END_REF][START_REF] Parikh | Objectively diagnosed diffuse white matter abnormality at term is an independent predictor of cognitive and language outcomes in infants born very preterm[END_REF][START_REF] Calloni | Neurodevelopmental outcome at 36 months in very low birth weight premature infants with MR diffuse excessive high signal intensity (DEHSI) of cerebral white matter[END_REF][START_REF] De Bruïne | Clinical implications of MR imaging findings in the white matter in very preterm infants: a 2-year follow-up study[END_REF]. However, apart from Parikh et al. study that has used a T2 relaxometry algorithm to detect white matter abnormality in centrum semiovale, the DEHSI detection was performed with a qualitative visual assessment by one or two observers, in a limited size of cohort.

To address these issues, we have developed a specific algorithm to automatically obtain brain MRI metrics in preterm neonates, based on a mathematical morphology approach with max-tree representations and a supplementary quantitative evaluation of diffuse high signal intensity regions of the white matter on T2 weighted images [START_REF] Xu | The challenge of cerebral magnetic resonance imaging in neonates: A new method using mathematical morphology for the segmentation of structures including diffuse excessive high signal intensities[END_REF]. It has been already validated on 1.5T and 3T MRI devices with a high reproducibility and robustness [START_REF] Xu | The challenge of cerebral magnetic resonance imaging in neonates: A new method using mathematical morphology for the segmentation of structures including diffuse excessive high signal intensities[END_REF].

Using this approach, the aim of our study was to identify in the large EPIRMEX cohort the correlations between MRI brain metrics, including DEHSI and the neurodevelopmental outcomes at 2 years.

Materials and methods

Population

Between June 28th, 2011 and October 31st, 2012, 519 preterm infants from 26 to 31 completed weeks of gestation from 16 university hospital centers in 12 regions were recruited at birth into a national prospective longitudinal cohort study (EPIRMEX cohort study associated with EPIPAGE 2). The inclusion was the same as the EPIPAGE2; eligible participants include all infants liveborn or stillborn between 22 and 31 completed weeks of gestation in all the maternity units in French during the inclusion period. The population eligible for follow-up includes all children alive at discharge from neonatal intensive care or special units or maternity wards whose parents have not declined to participate [START_REF] Ancel | EPIPAGE 2: a preterm birth cohort in France in 2011[END_REF]. Noninclusion criteria were children with severe karyotypic abnormalities, antenatal severe central nervous system malformations.

All survivors were enrolled for longitudinal follow-up and included in the study at 2 years corrected age if parents consented. Perinatal and maternal details were collected by chart review at the time of recruitment. Ninety-nine patients with incomplete clinical data were excluded from the statistical analysis. Thirty brain MRI at TEA had severe artifacts and were non-interpretable. Three hundred ninety brain MRI with complete clinical data and ages and stages questionnaire (ASQ) score were analyzed. Comparisons between excluded and included patients were performed. Two sub-groups have been defined depending on the magnetic field used (1.5 or 3T). The study was approved by the Human Research Ethic Committee (2011-10_9). Written informed parental consent was obtained for all infants.

MRI protocol

MRI brain scans were performed in natural sleep at TEA (i.e., GA of 39-41 weeks), using a 1.5T or 3T MRI system with a dedicated 8-channel head coil.MR devices with a magnetic field of 1.5T were Philips Achieva ® , Philips Intera ® , Toshiba MRT 200 ® , GE SignaHdxt ® (General Electric Healthcare), Siemens Avanto ® , Siemens Symphony ® , and Siemens SymphonyTim ® (Siemens Heathineers). The MRI device with a magnetic field of 3T was a Philips Achieva ® (Philips Healthcare). T2 datasets were obtained using an axial T2 morphological sequence (fast spin echo/turbo spin echo with a 90 flip-back pulse); slice thickness, 3 mm; pixel size, 0.39 × 0.39 mm 2 ; field of view, 192 mm; repetition time, 6680 ms; echo time, 142 ms; flip angle, 120°.

The axial MRI reference plane was the bi-commissural plane. MRI protocol also included also an axial T2*-weighted sequence, an axial diffusionweighted sequence (b = 1000 s/mm²), and a three-dimensional T1 MPRage sequence (spatial resolution, 0.8 × 0.8 × 1.2 mm 3 ; inversion time, 110 ms; repetition time, 2200 ms; echo time,

ms; flip angle, 8°).

A medical engineer visited all participating centers to check the sequence parameters.

The infants were fed, swaddled and had earplugs. No child has received medicated sedation.

Throughout the scan, infants were monitored using an apnea monitor and an oxygen saturation probe, and if required, oral sucrose was administered with parental consent.

Automated brain metrics

Anatomical structures

The previously described validated semi-automated mathematical morphology segmentation method was used by one observer on a single axial T2 section on which the caudate heads, lentiform nuclei, and thalami were maximally visible, as defined by Kidokoro et al. [START_REF] Kidokoro | New MR imaging assessment tool to define brain abnormalities in very preterm infants at term[END_REF]. A quality control was performed by a secondary blind observer to determine the most appropriate slice. A qualitative assessment of the absence of morphologic sequelae of hemorrhage and cystic lesions of the segmented MRI axial slice was performed by a pediatric radiologist prior to segmentation.

The segmentation method sequentially extracted each brain tissue from the axial T2weighted neonatal brain images. The algorithm relied on some morphological methods, notably the max-tree representation proposed by Salembier et al. [START_REF] Salembier | Antiextensive connected operators for image and sequence processing[END_REF], which is a hierarchical representation of the image based on threshold decomposition. First the intracranial cavity (ICC) was extracted as the region of interest and then the cerebrospinal fluid and ventricles using the max-tree and the context-based energy based on the selection strategy relying on markers. The basal ganglia and thalami (BGT) were extracted using a modified max-tree built. Then the gray matter and white matter (WM) were separated from the remaining tissues based on the optimal histogram-based thresholding. White matter hyperintensities (i.e., DEHSI) were segmented using the max-tree and the context-based energy again [START_REF] Xu | The challenge of cerebral magnetic resonance imaging in neonates: A new method using mathematical morphology for the segmentation of structures including diffuse excessive high signal intensities[END_REF].

Comparisons on the segmentation results on premature scanned at 1.5T and 3T were performed.

Detection and quantification of white matter signal abnormalities (DEHSI)

Within the segmented white matter, DEHSI were automatically detected with three different thresholds of intensity of 80%, 85% and 90% of the maximal intensity signal of the MRI (where 100% corresponds to the cerebrospinal fluid intensity signal). The corresponding surfaces were measured in mm². A qualitative visual inspection of the segmentation results has been performed by a second observer. No delineation error was reported. The considered method was observed to perform well in the presence of common brain abnormalities in preterm infants, including enlarged ventricles and WM signal intensity abnormalities. Results of the segmentation were automatically displayed and tabulated. Further analyses were performed on the results obtained on two-subgroups: premature scanned at 1.5T and 3T.

Two-year developmental assessment

All parents were invited to attend a developmental assessment at 2 years' corrected age. Children's development was assessed with the second version of the 24-month ASQ [START_REF] Squires | Revision of a parent-completed development screening tool: ages and stages questionnaires[END_REF][START_REF] Flamant | Parent-completed developmental screening in premature children: a valid tool for follow-up programs[END_REF]. Data were analyzed when completed between 22-and 26-months corrected age in children without deafness, blindness, or severe congenital brain malformations. Each questionnaire includes 30 items covering five developmental domains: communication abilities, gross motor skills, fine motor skills, problem solving abilities, and personal-social skills. Each developmental domain comprises 6 items scored on a three-point scale depending on whether the child performs the task: yes (10 points), sometimes (5 points), or not yet (0 points). Responses are summed to give a score of 0 to 60 for each domain and an overall maximum ASQ score of 300 points. Analyses were based on domain specific scores [START_REF]Council on Children with Disabilities, Section on developmental behavioral pediatrics, bright futures steering committee, medical home initiatives for children with special needs project advisory committee. Identifying infants and young children with developmental disorders in the medical home: an algorithm for developmental surveillance and screening[END_REF].

Statistical analysis

Qualitative variables were expressed as raw number and percentages. Percentages were reported with their exact 95% binomial confidence intervals. Quantitative variables were expressed as means ± standard deviations (SD) and ranges. Qualitative variables were compared using Chi squared test, and quantitative variables were compared using Student ttest after the evaluation of their distribution by a histogram.

The associations between brain MRI measures with global ASQ score and between BGT and ASQ sub-domains were explored using adjusted linear regression analyses (age of gestation, gender, corrected age during MRI, severe brain lesion such as grade III hemorrhage and cystic periventricular leukomalacia observed at cranial ultrasound). All tests were two sided. A P value < 0.05 was considered to indicate statistical significance. Statistical analyses were performed using R software [START_REF] Core | R: a language and environment for statistical computing[END_REF].

Results

Population and MRI

Of the 519 newborn participants included in this study, 489 premature newborns had brain MRI axial slice with normal findings. Nineteen (19/519; 3.66%) had sequelae of hemorrhage, 10 (10/519; 1.93%) had ventricular dilatation and 1 (1/519; 0.19%) had sequelae of cystic periventricular leukomalacia. In all newborns, segmentation was obtained.

Of the 519 newborn participants included in this study, 390 (390/519; 75.14%) patients had complete clinical data and brain MRI suitable for analysis and DEHSI detection. Participants were similar in all perinatal variables compared with excluded patients (Table 1).

The mean birth weight of our cohort was 1174 ± 339 (SD) g (range: 443-2380 g), the mean adjusted gestational age was 28.57 ± 1.82 (SD) (range: 24-33) weeks of amenorrhea (WA) and the mean corrected age for the MRI at TEA was performed at 39.8 ± 1.40 (SD) WA (range: 38-42 WA).

Distribution and characteristics of newborn scanned at 1.5T (n = 338) and 3T (n = 52) are detailed in the Table 1. Figure 1 shows the flow chart of inclusion and exclusion of patients

Two-year developmental assessment

Associations between MRI measures and 2-year development

Brain metrics

Metrics were obtained in all preterm infants (Figure 2). The comparison of the metrics obtained in newborns scanned at 1.5T and 3T are reported in the Table 3. Distribution of grey and white matter was different in the subgroups (Table 3). Correlations between the segmented brain MRI structure surfaces with the ASQ score evaluated by adjusted linear regressions are reported in Table 4.

Significant positive correlations with global ASQ score were found for intracranial cavity areas (β = 0.010; P = 0.008), basal ganglia and thalami (β = 0.091; P = 0.0001), gray (β = 0.014; P = 0.043) and white matter surfaces (β = 0.012; P = 0.015), where β represents the correlation between the area and the ASQ score, and P the degree of significance of the correlation. The more parenchyma was present, the better are the ASQ scores. Cerebrospinal fluid and lateral ventricles significantly negatively correlated with global ASQ score (β = -0.021, P = 0.047 and β = -0.088, P = 0.005 respectively). The more cerebrospinal fluid was present; the worse were the ASQ scores. For all ASQ sub-domains (communication abilities, gross motor skills, fine motor skills, problem solving abilities, and personal-social skills), BGT surfaces significantly correlated with the scores (Table 3).

Diffuse excessive high signal intensity

DEHSI was automatically detected in 289 of 390 (74%) included premature newborns with a threshold of 80% of highest intensity of the white matter, 219 (56%) with a threshold of 85% and 144 (37%) with a threshold of 90%. Metrics of DEHSI were higher at 3T, whatever the thresholds. An illustration of the segmentation results with a progressive increase of thresholds (80, 85 and 90%) is provided in Figure 3. With the highest threshold of 90%, the DEHSI corresponded to the periventricular crossroad areas.

Regardless of the intensity threshold chosen (80%, 85% or 90%) or the magnetic field used, no significant correlation between DEHSI detected by the software and ASQ at 2 years was observed in our cohort (β = -0.003, P = 0.664; β = -0.005, P = 0.607; β = -0.003, P = 0.834). Whisker plots showed mean ASQ score depending on the DEHSI threshold detection (Figure 4).

Discussion

Automated brain metrics provide reliable prognostic information with significant correlation with the neurodevelopmental evolution evaluated with the ASQ score at 2 years of age.

Correlations were particularly observed with the surface of BGT; the higher the surface, the better were the neurodevelopmental outcomes. By contrast, the focus on DEHSI did not show any correlation with the ASQ score, regardless of the intensity threshold chosen or the magnetic field, based on our quantitative detection.

Deep gray matter evaluated by the BGT surfaces significantly correlated with the ASQ score and all its sub-domains. Deep grey matter growth rates in preterm newborns are promising biomarkers of long-term outcomes [START_REF] Kidokoro | New MR imaging assessment tool to define brain abnormalities in very preterm infants at term[END_REF][START_REF] Young | Deep grey matter growth predicts neurodevelopmental outcomes in very preterm children[END_REF]. Reduced BGT volumes and slower growth rates in very preterm newborns correlated with distinct neurological dysfunctions at 7 years [START_REF] Loh | Longitudinal growth of the basal ganglia and thalamus in very preterm children[END_REF].

Even in our preterm newborns with a "normal" size of lateral ventricles on imaging, the more cerebrospinal fluid was observed, the lower the ASQ score was. Fox et al. also found a relationship between ventricular size at 1 month and outcome in infants less than 30 weeks' gestation [START_REF] Fox | The relationship between ventricular size at 1 month and outcome at 2 years in infants less than 30 weeks' gestation[END_REF]. In this regard, larger lateral ventricles in the parietal region were related to poorer motor development at 2 years, and to slower early language development [START_REF] Fox | The relationship between ventricular size at 1 month and outcome at 2 years in infants less than 30 weeks' gestation[END_REF].

Cortical gray matter and white matter surfaces also correlated with the ASQ score. Pittet et. al. have suggested that the brain maturity at term age explored by the radiological aspect of the transient fetal compartments (periventricular crossroad areas, von Monakow segments, subplate compartments) can be assessed on regular T2-weighted MRI [START_REF] Pittet | Newborns and preterm infants at term equivalent age: a semi-quantitative assessment of cerebral maturity[END_REF]. They found a maturity difference in newborn compared to premature newborn at TEA [START_REF] Pittet | Newborns and preterm infants at term equivalent age: a semi-quantitative assessment of cerebral maturity[END_REF], but with the limitation of a subjective visual assessment [START_REF] Morel | Neonatal brain MRI: how reliable is the radiologist's eye?[END_REF]. Our results suggest that an automated software would be a useful tool to analyze brain MRI more objectively and that it would provide interesting data to predict the neurodevelopmental outcome at two years. While the different thresholds of detection used have influenced the measured DEHSI rates, they had not influenced the lack of correlation observed with the ASQ score. The predictive value of DEHSI did not vary with the intensity threshold. Based on the presented results, without topographical analysis, we hypothesize that the DEHSI might no longer be described as a pathological process with poor prognosis. As all MRI sequence parameters used were homogenized by a medical physicist, no influence of the type of MRI scanner should be expected. However, we did observe an impact of the magnetic field on some particular metrics. Differences in distribution between white and gray matter could be explained by a higher resolution at 3T that would limit the partial volume effect observed at the proximity of the cortical gray matter. This phenomenon was already described in segmentation review.

Another effect of the magnetic field was greater metrics of DEHSI segmented at 3T. It could probably be explained by an increased signal-to-noise ratio and a greater contrast resolution that might have led to a higher sensitivity for revealing DESHI. However, the increased metrics DEHSI did not correlate with a lower ASQ in our cohort. A meta-analysis has recently suggested that DEHSI on TEA-MRI does not predict future development of neurodevelopmental disabilities [START_REF] Rath | Diffuse excessive high signal intensity on term equivalent MRI does not predict disability: a systematic review and meta-analysis[END_REF]. This meta-analysis was mainly based on studies that diagnosed DEHSI on visual impression on conventional T2 imaging [START_REF] Rath | Diffuse excessive high signal intensity on term equivalent MRI does not predict disability: a systematic review and meta-analysis[END_REF].

Agreed by the American Academy of Paediatrics, the ASQ is a worldwide valid parent completed developmental screening test [START_REF]Council on Children with Disabilities, Section on developmental behavioral pediatrics, bright futures steering committee, medical home initiatives for children with special needs project advisory committee. Identifying infants and young children with developmental disorders in the medical home: an algorithm for developmental surveillance and screening[END_REF]. ASQ, which is available in France [START_REF] Flamant | Parent-completed developmental screening in premature children: a valid tool for follow-up programs[END_REF], has been presented with arguments to support the cross-cultural validity for other European countries [START_REF] Janson | Parent-completed developmental screening in a Norwegian population sample: a comparison with US normative data[END_REF]. It allows the parents to be genuinely and centrally involved in the follow-up of their child. Nevertheless, the parents involved in our study may have under-or over-estimated the capabilities of their child. ASQ is a less resource-consuming tool to assess neurodevelopmental outcome but is less sensitive and specific than Bailey III. Potential social and cultural differences between parents are not factors deemed to alter the validity of their neurologic assessment of their child [START_REF] Flamant | Parent-completed developmental screening in premature children: a valid tool for follow-up programs[END_REF]. However, the ASQ questionnaire enabled us to detect more subtle neurological abnormalities by exploring each area of neurologic development.

In the present study, we have segmented one slice per newborn, which was representative of the whole brain, as all main supratentorial brain structures were analyzed, as proposed by Kidokoro et al. [START_REF] Kidokoro | New MR imaging assessment tool to define brain abnormalities in very preterm infants at term[END_REF]. Our metrics results have revealed a correlation between the measured surfaces and the neurodevelopmental evolution, as it was observed in three-dimensional segmentation studies. Our algorithm provides robust metrics results with a processing time less than a minute, without a specific MRI sequence or an experienced user, that could be performed in daily practice. It encourages continuing to use validated automated brain metrics during childhood [START_REF] Serru | Improving diagnosis accuracy of brain volume abnormalities during childhood with an automated MP2RAGE-based MRI brain segmentation[END_REF].

A limitation of our study would be a partial exploration of the supratentorial brain stage. Our metrics were a complementary tool to analyze the brain MRI, but a conventional qualitative interpretation of the brain MRI is still required. The evaluation of the volume of the corpus callosum would be of interest but not available in our segmentation algorithm. To analyze neurodevelopmental outcomes further, a detailed segmentation of each specific basal ganglia and the cerebellum evaluation would be necessary. Reproducibility was not evaluated in this study. However, we have obtained good to excellent intra-and inter-observer Dice similarity coefficients while using the algorithm on 1.5T and 3T brain MRI in a previous study [START_REF] Xu | The challenge of cerebral magnetic resonance imaging in neonates: A new method using mathematical morphology for the segmentation of structures including diffuse excessive high signal intensities[END_REF].

The image quality might influence the Dice similarity coefficients. Indeed, lower quality exams induced a slightly lower coefficient and a higher variance [START_REF] Serru | Improving diagnosis accuracy of brain volume abnormalities during childhood with an automated MP2RAGE-based MRI brain segmentation[END_REF]. As observed in a recent data challenge, although 3D segmentation is useful clinically, the choice of a bi-dimensional segmentation had simplified data collection, annotation, and storage without leading to significant difference of the Dice similarity coefficient among the participants [START_REF] Blanc-Durand | Abdominal musculature segmentation and surface prediction from CT using deep learning for sarcopenia assessment[END_REF]. Another limitation is that the key period for neuropsychometric assessment is around 6 years. As the initial research protocol defined a two-year examination, we have used these data. However, some of the children included in the cohort still have a clinical follow-up. Complementary research studies will be performed on the available long-term follow-up data, despite the possibility of a smaller cohort due to patients lost to follow-up.

In conclusion, in our large prospective, multicentric cohort of preterm newborns, the use of our software revealed a significant correlation between reproducible metrics of the main brain structures with the neurodevelopment at two years of age. Furthermore, the quantitative detection and segmentation of isolated DEHSI of the white matter have allowed excluding any correlation with the neurodevelopment at two years of age. Automated brain MRI metrics are now suitable in daily practice for neonatologists and radiologists.

At 2

 2 years of age, 390 (390/519; 75.14%) children completed the ASQ score. The mean overall ASQ score was 221.5 ± 53.4 (SD) (range: 30-300) out of 300 points. The detailed results of the ASQ sub-domains are presented in Table 2. The gender was an important adjusted coefficient: boys had a significantly lower ASQ than girls (linear adjusted coefficient of -16 for the BGT surface associated with the global ASQ score; P=0.002). There was no difference in the ASQ score of newborns scanned at 1.5T (222 ± 53 [SD]; range: 30-300) or 3T (215 ± 53 [SD]; range: 20-295) (P = 0.404).

Acknowledgments:

The authors thank John Sheath for English language assistance and his friendly support and the EPIRMEX study group: Catherine ADAMSBAUM, Pierre-Yves ANCEL, Catherine ARNAUD, Laurent BARANTIN, Olivier BAUD, Nathalie BEDNAREK, Valérie BIRAN, Olivier BRISSAUD, Aude CHAROLLAIS, Katia CHAUMOITRE, Thierry DEBILLON, Denis DEVICTOR, Patrice GILLET, Catherine GIRE, Bernard GUILLOIS, Sophie GUYETANT, Jean-Michel HASCOET, Petra HUPPI, Stéphane MARRET, Sylvie NGUYEN, Véronique PIERRAT, Patrick PLADYS, Jean Christophe ROZE, Umberto SIMEONI, Dominique SIRINELLI, Véronique ZUPAN

Funding source: This study did not receive any specific funding Author contributions. All authors attest that they meet the current International Committee of Medical Journal Editors (ICMJE) criteria for Authorship.

Guarantors of integrity of entire study, BM, CA, CT, ES; Study concepts/study design or data acquisition or data analysis/interpretation, all authors; Manuscript drafting or manuscript revision for important intellectual content, all authors; Approval of final version of submitted manuscript, all authors; Agrees to ensure any questions related to the work are appropriately resolved, all authors; Literature research, BM, PB, CT, ES; Clinical studies, CG, SM; NK; Statistical analysis, BM, ET; Manuscript editing, all authors.

Human rights. All procedures performed were in accordance with the ethical standards of the institutional and/or national research committee and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. The authors declare that the work described has not involved experimentation on humans. The local institutional review board approved this single multicentric prospective cohort study.

Informed consent and patient details. The authors declare that this report does not contain any personal information that could lead to the identification of the patient. Informed consent was obtained from all individual participants included in the study. The local institutional review board approved this single multicentric prospective cohort study.

Conflict of interest:

The authors have no conflicts of interest or financial relationships relevant to this article to disclose.

Legends