S Chevallier
email: sylvain.chevallier@uvsq.fr

E K Kalunga

Q Barthélemy

E Monacelli

Review of Riemannian distances and divergences, applied to SSVEP-based BCI

Keywords: Riemannian geometry, covariance matrices, distances, divergences, EEG, SSVEP

Background The first generation of brain-computer interfaces (BCI) classifies multi-channel electroencephalographic (EEG) signals, enhanced by optimized spatial filters. The second generation directly classifies covariance matrices estimated on EEG signals, based on straightforward algorithms such as the minimum-distance-to-Riemannian-mean (MDRM). Classification results vary greatly depending on the chosen Riemannian distance or divergence, whose definitions and reference implementations are spread across a wide mathematical literature.

Methods This paper reviews all the Riemannian distances and divergences to process covariance matrices, with an implementation compatible with BCI constraints. The impact of using different metrics is assessed on a steady-state visually evoked potentials (SSVEP) dataset, evaluating centers of classes and classification accuracy.

Results and Conclusions Riemannian approaches embed crucial properties to process EEG data. The Riemannian centers of classes outperform Euclidean ones both in offline and online setups. Some Riemannian

Introduction

Brain-computer interfaces (BCI) provide humans and animals with the ability to interact with machines without the use of their motor functions. They offer alternatives for the compensation of physical limitations (Mc-Farland et al. 2017), with applications in rehabilitation and assistive technologies [START_REF] Hosni | EEG-EOG based virtual keyboard: Toward hybrid brain computer interface[END_REF]. Brain signals could be recorded with various systems, with the electroencephalography (EEG) being widely adopted [START_REF] Niedermeyer | Electroencephalography: Basic Principles, Clinical Applications, and Related Fields[END_REF] as it has a low cost and requires low computational power. BCI designs commonly rely on one of the three following paradigms [START_REF] Nam | Braincomputer interfaces handbook: technological and theoretical advances[END_REF][START_REF] Quitadamo | Describing different brain computer interface systems through a unique model: A UML implementation[END_REF]: the event-related (de)-synchronization (ERD/S), the event-related potentials (ERP) and the steady-state visually evoked potentials (SSVEP). ERD/S are detected in the premotor cortex and are correlated with motor imagery. ERP are the observed brain responses to an external stimulation of the sensory nervous system, often auditory or visual in BCI. SSVEP are elicited by the attentional focus on a steady pattern of a visual stimulus [START_REF] Wang | Lead selection for SSVEP-based brain-computer interface[END_REF]). Auditory or somatosensory patterns could be considered, but visual stimulation is the most common in BCI [START_REF] Zhang | Ssvep recognition using common feature analysis in braincomputer interface[END_REF]. In SSVEP paradigm, a stimulus flickering at constant frequency is generated by LEDs or by alternating images on a computer screen. When a user's gaze is directed toward a flickering stimulus, it elicits cortical oscillations with matching fre-quency. Those cortical oscillations are modulated by the selective attention [START_REF] Vialatte | Steady-state visually evoked potentials: focus on essential paradigms and future perspectives[END_REF]. The experiments conducted in this paper focus on SSVEPbased BCI due to their general robustness and efficiency [START_REF] Volosyak | BCI demographics II: How many (and what kinds of) people can use a high-frequency SSVEP BCI[END_REF].

BCI face important challenges that are mainly caused by the poor spatial resolution of EEG. This limitation is due to the volume conduction effect [START_REF] Niedermeyer | Electroencephalography: Basic Principles, Clinical Applications, and Related Fields[END_REF]; the bone structure filters out the higher frequency part, mixing the brain sources signals and reducing the signal-to-noise ratio (SNR). Consequently, spatial filtering methods are developed to enhance SNR, such as common spatial patterns (CSP) [START_REF] Johannes | Designing optimal spatial filters for single-trial EEG classification in a movement task[END_REF], independent component analysis (ICA) [START_REF] Wang | Enhancing Evoked Responses for BCI Through Advanced ICA Techniques[END_REF], xDAWN [START_REF] Rivet | xDAWN Algorithm to Enhance Evoked Potentials: Application to Brain-Computer Interface[END_REF] or canonical correlation analysis (CCA) [START_REF] Yang | Removal of electrooculogram artifacts from electroencephalogram using canonical correlation analysis with ensemble empirical mode decomposition[END_REF]. Filtered signals define a feature space where machine learning methods are used to classify trials. While spatial filters are very efficient on clean datasets obtained in constrained environments, they are easily distorted by artifacts and outliers [START_REF] Tomioka | Logistic regression for single trial EEG classification[END_REF], and are specific to subjects and sessions.

All the above-mentioned spatial filters rely on complex computations based on covariance matrices estimated from the EEG signal. Covariance matrices have a specific structure and are defined on a subspace of the Euclidean space. Classical matrix computationssuch as Euclidean distance, 2 norm or arithmetic mean -results in inaccurate results when applied directly on covariance matrices. It also induces numerical artifacts, like the swelling effect [START_REF] Arsigny | Geometric means in a novel vector space structure on symmetric positive-definite matrices[END_REF]). The Riemannian framework accounts for the specific geometry of the space of covariance matrices. It embeds specific properties, like affine-invariance, resulting in equivalent distances between sensor-based covariance matrices or between source-based covariance matrices (Congedo et al. 2017a;[START_REF] Yger | Riemannian approaches in brain-computer interfaces: a review[END_REF]. This invariance to full-rank spatial filtering provides good generalization capabilities to classification algorithms. Recent works [START_REF] Barachant | Multiclass brain-computer interface classification by Riemannian geometry[END_REF][START_REF] Kalunga | Online SSVEP-based BCI using Riemannian geometry[END_REF][START_REF] Meinel | Characterizing regularization techniques for spatial filter optimization in oscillatory EEG regression problems[END_REF]) demonstrate the benefit of Riemannian framework. Riemannian approaches are also applied in the winning strategy of five BCI competitions (Congedo et al. 2017a). For EEG classification, a very straightforward approach introduced in [START_REF] Barachant | Multiclass brain-computer interface classification by Riemannian geometry[END_REF] is the known as minimum distance to Riemannian mean (MDRM). The class centers are defined as the average of covariance matrices for each class and an unseen trial is assigned to the class with the nearest center, in the sense of a Riemannian distance. This straightforward approach is robust and yields better accuracy than more sophisticated classifiers of the litera-ture [START_REF] Barachant | Multiclass brain-computer interface classification by Riemannian geometry[END_REF], especially on systems with few electrodes [START_REF] Yang | Subject-specific channel selection using time information for motor imagery brain-computer interfaces[END_REF].

The space of covariance matrices is curved and proper distances are geodesics, i.e. curves that flow along the surface of the space. Different distances are defined in the mathematical literature, each with specific properties that could be applied to the space of covariance matrices [START_REF] Arsigny | Geometric means in a novel vector space structure on symmetric positive-definite matrices[END_REF]. Still, there is no prior work that indicates the most appropriate distances to process EEG signals, or even to propose a thorough and reproducible analysis.

While most applications of Riemannian geometry on EEG are dedicated to motor imagery and ERP paradigms (Congedo et al. 2017a;[START_REF] Yger | Riemannian approaches in brain-computer interfaces: a review[END_REF], this work focuses on SSVEP-based BCI. The experimental part is applied on a dataset acquired in the context of assistive robotics, aiming toward applications for cerebral palsy and post-stroke rehabilitation. The task follows a shared control paradigm, where an arm exoskeleton is operated through a controller shared between an SSVEP-based BCI and a 3D touchless interface based on infrared sensors [START_REF] Martin | Fast calibration of hand movement-based interface for arm exoskeleton control[END_REF][START_REF] Kalunga | Hybrid interface: Integrating BCI in multimodal humanmachine interfaces[END_REF].

The relevant properties to cope with brain signal variability have not been identified, and thus EEG signal processing is lacking the proper metrics to conduct deeper investigation. To address several open questions in EEG applications, such as handling specific noise, change in electrode positioning, inter and intraindividual variability, advances in Riemannian BCI could open new opportunities to develop applicative solutions. Such opportunities require a thorough analysis of existing distances and their properties, in a reproducible context.

Advances in Riemannian BCI have opened new opportunities to address several open questions in EEG applications, such as handling specific noise, change in electrodes positioning, coping with inter and intraindividual variability. However, the relevant metrics and their properties that could encompass the brain signal variability have not been formally identified. It is therefore important to identify the relevant metrics through an extended analysis in a reproducible context.

While we start to investigate these issues in a previous work (Kalunga et al. 2015), this paper extends this former approach and the contributions are the following:

a review of classification methods for SSVEP-based BCI, showing remaining challenges such as restingstate class, a complete survey of the existing distances and divergences for covariance matrices for the first time in the literature, an evaluation of the computational cost, numerical stability and a qualitative analysis of results of those distances and divergences, an open framework which includes the algorithms in Matlab (Kalunga 2018(Kalunga , 2015) ) and Python (Chevallier 2020; Bertrand-Lalo 2020), and the datasets (Chevallier 2017) to ensure the reproducibility of experiments.

This paper is organized as follows. Section 2 details the neurophysiological ground of the steady-state visually evoked potentials and reviews the existing methods for their classification. Section 3 describes the framework for the classification of covariance matrices, including the considered distances and divergences, and presents the minimum distance to mean (MDM), which is the classification algorithm used in this work. The real EEG dataset used for the experiment is presented in Section 4. In Section 5, the results of both offline and online analyses are shown and discussed. Section 6 concludes this paper.

SSVEP-based BCI

Stimulus design

Steady-state visual evoked potential is a brain response to a repeated visual stimulus. It is measured in the primary visual cortex as oscillations at a frequency which corresponds to the visual stimulus frequency and its harmonics [START_REF] Vialatte | Steady-state visually evoked potentials: focus on essential paradigms and future perspectives[END_REF]. In BCI, while various possible stimuli design are possible [START_REF] Gergondet | SSVEP stimuli design for object-centric BCI[END_REF], SSVEP is triggered by a flashing visual target using either LEDs or a computer screen. The LEDs induce larger oscillation amplitude than a computer screen [START_REF] Zhu | A survey of stimulation methods used in SSVEP-based BCIs[END_REF]. A classical setup with LEDs is visible on Fig. 1. Computer screens also have a limitation in terms of the flashing frequencies that could be produced, as it should be a multiple of R, the screen refresh rate. This limitation is overcome by the sampled sinusoidal stimulation method that allows to simulate intermediate frequencies. These intermediate frequencies are obtained by modulating the luminance of the screen using the following equation [START_REF] Chen | A high-ITR SSVEPbased BCI speller[END_REF]:

s(t, f, θ) = 1 2 1 + sin 2πf (t/R) + θ ,
where s(t, f, θ) is the luminance of the screen at time t, varying with frequency f and phase θ. In this stimulation mode, a phase shift θ is added to differentiate between target stimulus with identical oscillation frequencies, allowing for more possible targets. Nonetheless, paradigms that require the use of phase information add a constraint on the BCI system: the SSVEP stimulation device should be synchronized with the EEG acquisition system, usually with wire cables. Any frequency, between 2 Hz and 50 Hz, induces visible oscillations in the visual cortex [START_REF] Herrmann | Human EEG responses to 1-100 Hz flicker: resonance phenomena in visual cortex and their potential correlation to cognitive phenomena[END_REF], the amplitude of the oscillations decreases as the frequency increase. Frequencies as high as 90 Hz can be used for SSVEP stimulation and could improve the visual comfort of the subject (Chen et al. 2015a), but are difficult to detect. Visual stimuli between 12 Hz and 25 Hz are the most common in SSVEP-BCI as this range provides a good compromise between detectability and reactivity [START_REF] Zhu | A survey of stimulation methods used in SSVEP-based BCIs[END_REF]).

Classification with CCA

Classification of SSVEP brain signals relies on the identification of stimulus-specific frequency in the visual cortex. Before advances introduced by Riemannian frameworks, methods relying on CCA [START_REF] Lin | Frequency recognition based on canonical correlation analysis for SSVEP-based BCIs[END_REF]) achieved the highest classification performances. CCA finds a projection space that maximizes the cross-covariance of two sets of input signals while jointly minimizing their covariance [START_REF] Kalunga | SSVEP enhancement based on Canonical Correlation Analysis to improve BCI performances[END_REF]. The common methodology is to find the canonical space between a multichannel EEG trial X, on the one hand, and the reference signal Y f , on the other hand. Y f is usually defined as:

Y f =        sin(2πf t) cos(2πf t) . . . sin(2πN h f t) cos(2πN h f t)       
, where N h is the number of harmonics, and frequency f belongs to the set of stimulation frequencies. The goal of the CCA is to maximize the correlation ρ f between w T X X and w T Y Y f as follows:

ρ f = max w X ,w Y w T X Σ XY f w Y w T X Σ X w X w T Y Σ Y f w Y , (1)
where Σ X and Σ Y f are the covariance matrices of X and Y f , and Σ XY f is the cross-correlation matrix of X and Y f . The estimated frequency f is computed from the values of ρ f , as in [START_REF] Lin | Frequency recognition based on canonical correlation analysis for SSVEP-based BCIs[END_REF]:

f = arg max f {ρ f } . (2)
It is possible to optimize the results by carefully selecting the reference signals and its harmonics, but there is no specific recommendation. With Filter Bank Canonical Correlation Analysis (FBCCA) (Chen et al. 2015a), an input signal is decomposed in n = 1 . . . N b sub-bands with overlapping frequency bands. The correlation coefficients ρ n f are computed for each sub-band n and are combined as:

ρf = N b n=1 w n (ρ n f ) 2 with w n = n -a + b , (3)
where a and b are constants defining the weight w n of each sub-band. To maximize the classification performance, an offline grid search is needed to select the correct hyper-parameters a, b, and N b . The classification decision is reached with Eq. ( 2) as soon as the all the ρf are computed. There is a significant improvement in terms of accuracy from the classical CCA with this filter bank scheme.

An important milestone has been reached with studies conducted by [START_REF] Nakanishi | High-speed brain speller using steady-state visual evoked potentials[END_REF] and Chen et al. (2015b) for SSVEP-based spellers. These SSVEP-based BCI achieve the highest information transfer rate (Mc-Farland et al. 2003) reported in the EEG literature: [START_REF] Nakanishi | High-speed brain speller using steady-state visual evoked potentials[END_REF] brought the information transfer rate to an average of 166.91 bits/min and later, Chen et al. (2015b) raised it to an average of 270 bits/min.

These very high performances rely mostly on the improvement of two key aspects: the stimulation protocol, using phase modulation and combination of frequencies, and stimulus presentation, using a scheme similar to SMS writing. These approaches comes with two strong limitations: the first drawback is that there is no resting-state (also known as reject class, or no-SSVEP class). Even without any intent from the subject, the BCI system is selecting letters without any pause. The ability for a user to act at their own pace is crucial in Human-Machine Interface and in BCI, a system without any resting-state is hardly compatible for a real implementation of a BCI system. It should also be emphasized that handling the diversity of brain signals corresponding to resting-state is a challenging issue as stimulus frequency could still be observed in the EEG. The second limitation is the complexity of the protocol, which requires a hardware synchronization between the stimulation device and the EEG acquisition system for precise phase measurements. This limits potential applications.

A complete review of classification methods for SSVEPbased BCI can be found in [START_REF] Lotte | A review of classification algorithms for EEG-based brain-computer interfaces: a 10 year update[END_REF].

3 Riemannian BCI Spatial filters and classifiers in BCI rely on covariance matrices, as it could be seen in Eq. ( 1) for CCA. It is possible to operate directly on the covariance matrices of EEG signals, with formal tools of differential geometry and algorithms borrowed from the mathematical community of information geometry (Nielsen and Bhatia 2012). This section introduces the covariance matrices, their geometry and the two most important notions for EEG signal processing methods, which are the computation of distance between two covariance matrices and the estimation of the mean of a set of covariance matrices. Endowed with a distance and an average, it is possible to define all the required methods for processing and classifying EEG signals in the space of covariance matrices.

From EEG signals to covariance matrices

In the context of EEG signal, we will denote by X an EEG signal in R C×N where C is the number of channels, i.e. sensors that are placed on the subject's scalp, and N is the number of samples, which are the observations recorded by the system, with N > C. The covariance matrix is estimated with

Σ = 1 N XX T . (4)
For noisy or truncated signal is possible to use more advanced estimator, an experimental analysis of covariance estimators applied to EEG signal is detailed in [START_REF] Kalunga | Online SSVEP-based BCI using Riemannian geometry[END_REF]. Covariance matrices are symmetric and positive-definite (SPD) and it is possible to define a manifold M of SPD C × C matrices [START_REF] Bhatia | Positive definite matrices[END_REF] as:

M = Σ ∈ R C×C | Σ = Σ T and z T Σz > 0, ∀z ∈ R C \0 .
The eigenvalue decomposition of Σ is where λ 1 , . . . , λ C are the eigenvalues and U the matrix of eigenvectors of Σ. The unique symmetric square root Σ 1 2 and symmetric inverse square root Σ -1 2 are defined as:

Σ = U diag(λ 1 , . . . , λ C ) U T ,
Σ 1 2 = U diag(λ 1 2 1 , . . . , λ 1 2 C ) U T , Σ -1 2 = U diag(λ -1 2 1 , . . . , λ -1 2 C ) U T ,
and the exponential Exp and logarithm Log matrix operators are defined as:

Exp(Σ) = U diag(exp(λ 1 ), . . . , exp(λ C )) U T , Log(Σ) = U diag(log(λ 1 ), . . . , log(λ C )) U T .
Following Congedo et al. (2017a), specific covariance matrices for SSVEP are designed to embed frequency information. Let X be an EEG trial measured during an SSVEP experiment with F stimuli blinking at different frequencies. The covariance matrices are estimated from an extended version of the input signal X:

X ∈ R C×N →    X freq 1 . . . X freq F    ∈ R F C×N , (5)
where a band-pass filter around the frequency freq f , with f = 1, . . . , F , is applied on the input signal X to obtain X freq f . Fig. 2 show the extended signal of one electrode, filtered at F = 3 different frequencies.

Hereafter, EEG signals always refer to the augmented signal of Eq. (5).

Distances and divergences

In Riemannian framework, several distances have been proposed for positive-definite matrices with different properties, yielding various impacts on signal processing and classification of brain signals. There is a large literature in differential geometry, probability and matrix calculus describing a great variety of distances for positive-definite matrices, whose definitions are outside the scope of this paper. In this large literature, there are many reformulations of existing distances and generalization. Some distances cannot be computed while others have very efficient implementations available; some are very sensitive to the perturbations while others have useful invariances for robustness to noise. This paper summarizes all this literature and provides the precise bibliographical references for the distance that could be of use for Riemannian BCI. Divergences and distances are measures of dissimilarity between two points in a space. A distance function d has the following properties for all Σ 1 , Σ 2 , Σ 3 ∈ M:

1. Non-negativity: d(Σ 1 , Σ 2 ) ≥ 0 , 2. Identity: d(Σ 1 , Σ 2 ) = 0 iff Σ 1 = Σ 2 , 3. Symmetry: d(Σ 1 , Σ 2 ) = d(Σ 2 , Σ 1 ) , 4. Triangular inequality: d(Σ 1 , Σ 3 ) ≤ d(Σ 1 , Σ 2 ) + d(Σ 2 , Σ 3 ) .
Divergences, which are generalization of squared distance, do not necessarily satisfy properties (3) and ( 4).

In this work, to compare covariance matrices, we consider several existing distances, namely the Euclidean distance, the harmonic distance [START_REF] Lim | Matrix power means and the Karcher mean[END_REF], the affine-invariant Riemannian distance [START_REF] Pennec | A Riemannian Framework for Tensor Computing[END_REF], the Log-Euclidean distance [START_REF] Arsigny | Geometric means in a novel vector space structure on symmetric positive-definite matrices[END_REF]), the Wasserstein distance [START_REF] Villani | Optimal transport: old and new[END_REF], and divergences, such as the Kullback-Leibler divergence [START_REF] Nielsen | Sided and symmetrized Bregman centroids[END_REF]) (equivalent to the log-det divergence [START_REF] Dhillon | Matrix nearness problems with Bregman divergences[END_REF] for centered multivariate Gaussian distributions), the Jeffreys divergence [START_REF] Jeffreys | An invariant form for the prior probability in estimation problems[END_REF], the S-divergence [START_REF] Sra | Positive definite matrices and the S-divergence[END_REF]) (equivalent to the Bhattacharyya divergence [START_REF] Chebbi | Means of Hermitian positivedefinite matrices based on the log-determinant αdivergence function[END_REF])), and the α-divergence family [START_REF] Cichocki | Families of alpha-beta-and gamma-divergences: Flexible and robust measures of similarities[END_REF][START_REF] Nielsen | On clustering histograms with k-means by using mixed α-divergences[END_REF]. All these distances and divergences are detailed in Supplementary Material.

Table 1 provides consistent definitions of all these distances along with pointers to most efficient implementations. It is important to select the required properties carefully for a specific problem, as it could provide a measure that is robust to noise and undesired changes.

Riemannian mean

Considering the set of covariance matrices {Σ i } i=1...I , the mean (also called barycenter or center of mass) Σ

dE = Σ1 -Σ2 F ΣE = 1 I I i=1 Σi Harmonic dH = Σ -1 1 -Σ -1 2 F ΣH = 1 I I i=1 Σ -1 i -1 (
Jeffreys DJ = 1 2 tr(Σ -1 2 Σ1) + tr(Σ -1 1 Σ2) -C ΣJ = Σ 1 2 E Σ-1 2 E ΣH Σ-1 2 E 1 2 Σ 1 2 E ( Moakher
Wasserstein dW = tr Σ1 + Σ2 -2 Σ 1/2 1 Σ2Σ 1/2 1 1/2 1/2
Eq. ( 19) in minimizes the dispersion of matrices Σ i and could be expressed as:

Σ = µ({Σ 1 , . . . , Σ I }) = arg min Σ∈M I i=1 d p (Σ i , Σ) , (6)
for p = 2 d p (•, •) is a distance and for p = 1 it is a divergence. Σ is unique and is denoted with various names in the literature, such as geometric mean, or Frechet mean. It is sometimes called Cartan mean or Karcher mean [START_REF] Ando | Geometric means[END_REF][START_REF] Lim | Matrix power means and the Karcher mean[END_REF]. Cartan had shown that a unique solution to (6) exists if all Σ i lie in a convex ball [START_REF] Cartan | Groupes simples clos et ouverts et géométrie riemannienne[END_REF].

The choice of the divergence or the distance defined in Eq. ( 6) has an important impact on the resulting mean, that will be explored hereafter. In particular, for the affine-invariant Riemannian (AIR) distance, the geometric mean is the unique matrix ΣAIR of the manifold M satisfying:

I i=1 Σ1/2 AIR Log( Σ-1/2 AIR Σ i Σ-1/2 AIR ) Σ1/2 AIR = I i=1 Log ΣAIR (Σ i ) = 0 ,
where Log ΣAIR (Σ i ) is the logarithmic map projecting matrix Σ i from manifold to tangent space at point ΣAIR . Existence and unicity of this Riemannian mean can be proved [START_REF] Moakher | A differential geometric approach to the geometric mean of symmetric positive-definite matrices[END_REF]Afsari 2011;[START_REF] Massart | Inductive means and sequences applied to online classification of eeg[END_REF]. A summary of distances and divergences considered in this study, as well as their associated means, is proposed in Table 1, and their properties are listed in Table 2.

Multi-class classification of covariance matrices for SSVEP

Minimum Distance to Mean (MDM) [START_REF] Barachant | Multiclass brain-computer interface classification by Riemannian geometry[END_REF] allows to directly classify covariance matrices in

Σk = µ({Σ i } i∈I(k) ) , (7)
where I(k) is the set of indices of training data belonging to class k, and µ(.) is the mean estimation function.

For generalization, a new matrix Σ is assigned to the class of the closest mean, as illustrated on Fig. 3:

k = arg min k=1...K d(Σ, Σk ) , (8)
where d(.) is one of the distances from Table 1. The classification of SSVEP signals discriminates between K = F + 1 classes, to take into account the resting-state class. From I labeled training trials {X i } I i=1 recorded per subject, K centers of classes Σk are estimated. To obtain a correct estimation of the centers of classes, the artifacts are removed with the Riemannian potato [START_REF] Barachant | The Riemannian potato: an automatic and adaptive artifact detection method for online experiments using Riemannian geometry[END_REF]) used with a z-score threshold z th = 2.0 and a number of iterations L = 3.

The MDM algorithm described above is suitable for offline and block-online BCI implementations [START_REF] Chevallier | Brain Computer Interfaces Handbook: Technological and Theoretical Advances[END_REF]. Offline implementation requires a reference time, that is aligned with the cue onset. This reference time is used to estimate the covariance matrices, thus focusing on the SSVEP occurrences. In online and asynchronous setups, it is not possible to rely on such cue onset: the EEG epochs should be classified iteratively. An online computation scheme is applied on overlapping epochs, as described in [START_REF] Kalunga | Online SSVEP-based BCI using Riemannian geometry[END_REF], and the length of the trial is determined with a dynamic stopping approach [START_REF] Verschore | Dynamic stopping improves the speed and accuracy of a P300 speller[END_REF]). This avoids false detection induced by inattention of the subject or by external noise. The online computation scheme makes use of a curve direction criterion that shortens the synchronization latency by reducing the delays of the epoch-based decision process.

To conclude this section, MDM and MDM-online algorithms will be used in subsequent experiments to evaluate the classification accuracy of different distances and divergences.

Experimental SSVEP dataset

This section describes how the experiments have been conducted to obtain the data for the experimental analysis presented in the next section. The subject sits on an electric-powered wheelchair and an exoskeleton supports their right upper limb [START_REF] Martin | Fast calibration of hand movement-based interface for arm exoskeleton control[END_REF][START_REF] Kalunga | Hybrid interface: Integrating BCI in multimodal humanmachine interfaces[END_REF]. During the experimental recording, the exoskeleton is powered on but not used. On the left side of the subject, a panel shows 3 groups of 4 LEDs, blinking each at specific frequency (see Fig. 1). Although the panel is on the left side, the user could see it without moving their head. The subject is asked to follow audio instruction, while moving and blinking freely. In the protocol, a trial starts by an audio cue prompting the subject to gaze either at a group of LEDs or at a fixation point positioned outside the LED panel. This fixation point is used to characterize the reject class, which is a reference state without any associated SSVEP command.

A g.Mobilab+ is used for recording EEG, with C = 8 channels. The electrodes are positioned in the cap at PO7, PO8, PO3, PO4, POz, O1, O2 and Oz, with an earlobe reference and a ground electrode in Fz. The impedance is kept below 10 kΩ throughout the experiments. The EEG sampling rate is 256 Hz, with a bandpass filtering between 0.5 Hz and 100 Hz, and no spatial referencing is performed to avoid obtaining low-rank matrix. An example of EEG recording, extended as per Eq. ( 5), is shown on Fig. 2.

Three distinct visual stimuli (F = 3) characterized by their respective blinking frequencies of 13, 21, and 17 Hz, are presented to the subject. Thus, the classification task is set up for K = 4 classes, that is the F = 3 classes of the stimuli and the resting-state (no-SSVEP). Each session is made of 32 trials, divided in 8 visual stimuli for each class including the resting-state class. A trial is 4 second long. These experiments gather 12 subjects and, for each subject, there is between 2 and 5 sessions. The data are split into I = 32 trials for training set and the 32 to 128 trials left are for the test set. This dataset is accessible at [START_REF] Chevallier | SSVEP data[END_REF], with all the information to extract and process the EEG markers.

Experimental evaluation of distances and divergences

This section presents the analysis and experimental results of applying Euclidean and Riemannian distances and divergences in SSVEP classification task.

Assessment of centers of class

Covariance matrices of extended EEG signal, defined in Eq. ( 5), capture discriminative features of EEG generated by same brain sources but with different oscillating frequencies. The center of mass Σk for the K classes is shown on Fig. 4, with left (resp. right) part showing the subject with the highest classification results (resp. lowest). The diagonal blocks are the 8×8 covariance matrices computed for the F = 3 target frequencies, the offdiagonal blocks are almost null as there is a null interfrequency covariance. The block with the largest values is almost always indicating the target frequency. In the case of the resting state, all F blocks have equivalent values, with small absolute value. Indeed, this is most visible for the subject with the highest classification rate and more difficult to see for the one with the lowest classification rate. The mean covariance for 13 Hz is very similar to the resting-state, and the first block of the mean covariance for 21 Hz displays non-negligible values. For this subject, the mean covariance matrices are very informative but are not discriminative. The subject with the highest classification accuracy has relatively lower within frequency correlations: each class has a higher within frequency correlation in its corresponding target frequency, and the resting-state class has all within frequency correlations equally reduced. Consequently, the 4 classes are thus easily separable.

Classifiers commonly used in BCI, such as linear discriminant analysis (LDA) or support-vector machine (SVM) [START_REF] Hastie | The Elements of Statistical Learning[END_REF], can appear as black-boxes with decisions that are difficult to interpret. The presented classifier relies on features associated with a simple representation. This allows an intuitive understanding for the domain experts, such as medical doctors, neuroscientists, etc. One of the most interesting aspect is that there is a physiological meaning to the observed covariance matrices. Most of the complexity is embedded in a dedicated distance, allowing a simple interpretation, whereas machine learning algorithms produce complex decisions that are difficult to explain. Understanding why a subject shows some poor classification results is easily inferred from the values of the mean covariance matrices.

Fig. 5 shows the classification accuracy and CPU time, for all distances and divergences listed in Table 1, except for the Wasserstein which has a very high computational cost and a low accuracy. Firstly, Fig. 5a shows the computation time of averaging, to estimate the K = 4 centers of class Σk during MDM offline training. Secondly, Fig. 5b shows the computation time of distances and divergences elapsed on the trial classification of test data. The computation times are measured on every distance and divergence, then averaged across all subjects and trials. In an online setting where the target center of mass is updated at every trial, a classification using a specific metric would take the time to compute one distance/divergence plus the time to update the center of mass. These figures show that S-divergence and Bhattacharyya divergence, which are formally equivalent [START_REF] Bhattacharyya | On a measure of divergence between two statistical populations defined by their probability distribution[END_REF]; Sra 2016) but with different implementations, yield similar accuracy but different computation times. They also constitute a reference for compatibility in applications with realtime constrains.

Offline classification

On Fig. 6, the performances achieved with offline MDM using the above-mentioned distances and divergences are compared to the three state-of-the-art methods: [START_REF] Lin | Frequency Recognition Based on Canonical Correlation Analysis for SSVEP-Based BCIs[END_REF], [START_REF] Nakanishi | High-speed brain speller using steady-state visual evoked potentials[END_REF], and FBCCA of Chen et al. (2015a). All these methods are described in Section 2. Only classes corresponding to the stimulation frequencies (F = 3) are considered for the offline classification, K = 3, as CCA-based methods are not able to identify resting-state class. The harmonic number is set to N h = 6 for the CCA based methods, this value is selected after a grid search on the training data. The signal is decomposed in N b sub-bands for the FBCCA, where the sub-band n ranges from n × 8 Hz to 93 Hz. The maximum value is ten times the frequency span of the stimulation: in the described experiment, 8 Hz corresponding to the span from 13 Hz to 21 Hz. The parameters N b , a and b are selected by a grid search in the following range: N b = i N , with i N = 1, . . . , 7, a = 0.25 × i a , with i a = 0, 1, . . . , 40, b = 0.25 × i b , and with i b = 0, 1, . . . , 4. The selected values are N b = 3, a = 2 and b = 0.

Centers of class and classification performances

When evaluating the results of offline MDM, the arithmetic distance yields the poorest performance, far lower than the results obtained with the state-of-the art methods, as it could be seen on Fig 6. This poor performance is attributed to the lack of invariance under inversion and to the swelling effect (see Supplementary Material). The determinant value is tightly linked to the dispersion of multivariate variables, here the frequency bands and the EEG electrodes, thus leading to low classification accuracy. The swelling effect of the arithmetic mean is shown on the TraDe plot (Congedo et al. 2017b) of Fig. 7, which displays the log-trace as a function of the log-determinant for the different covariance matrices taken from the trials of the 13 Hz class of (b) Subject with lowest BCI performance. the subject with the highest BCI performance. Markers represent the different means and the color map represents the distribution of covariance matrices used in the computation of the means. The determinant of the arithmetic mean is strictly larger than any other mean. The Log-Euclidean, the AIR, the S-divergence, the Jeffreys, and the Bhattacharyya means yield similar determinants, close to the mean value of determinants of individual matrices. The determinant of the harmonic mean shows a reverse phenomenon: there is a shrinking effect on the determinant of the individual matrices. In fact, the lack of self-duality of arithmetic mean, i.e. ΣE (Σ i ) = Σ-1 E (Σ -1 i ), can be reformulated as ΣE = ΣH , where ΣH is the dual of ΣE (Congedo et al. 2017b). With the generalization capability of the power means Σt with t ∈ [-1, +1] (see Eq. ( 4) in Supplementary Material), the geometric mean (t → 0) is right between arithmetic (t = 1) and harmonic (t = -1) means. This justifies that the shrinking factor of the harmonic mean is of the same order as the swelling factor of the arithmetic mean. The same phenomenon can be observed for the trace, with a shrinking effect stronger for the Harmonic mean than for the AIR mean. The Log-Euclidean mean and the AIR one conserve their determinant, because their determinant is the geometric mean of the determinants of individual matrices, see Supplementary Material, whereas the arithmetic mean conserves the trace, because its trace is the arithmetic mean of the traces of individual matrices. The equivalence between the S-divergence and the Bhattacharyya is also verified on this plot. Finally, the Log-Euclidean mean appears as a good trade-off for the estimation of such properties, close to the center of the distribution.

Riemannian distances and divergences evaluation

Using Riemannian distances significantly improves the classification performances, in comparison with the stateof-the-art CCA-based methods which have an average performance across subjects of 87.5%, 81.2% and 87.4% for [START_REF] Lin | Frequency Recognition Based on Canonical Correlation Analysis for SSVEP-Based BCIs[END_REF], [START_REF] Lin | Frequency Recognition Based on Canonical Correlation Analysis for SSVEP-Based BCIs[END_REF] and Chen et al. (2015a) respectively. These methods rely on CCA coefficients for the classification: this is a strong limitation, as they could not account for the reject class. As the reject class is not associated with any reference signal, no correlation coefficient could be computed for this class. Due to this limitation, CCA-based methods avoid to confront with the most challenging case, that is to identify when the user does not look at any stimulus. These approaches are thus not compatible with realistic BCI use case.

Various Riemannian distances and divergences achieve different classification performances. It is possible to determine the Riemannian distance or divergence that is best suited for the classification of EEG covariance matrices, and for SSVEP in particular.

The α-divergence, being a generalization of the Kullback-Leibler divergence (see Supplementary Material), offers Fig. 6 Offline classification accuracy (in %) for K = 3 classes (no resting-state), by subject. Classification is performed with MDM using either Euclidean or Riemannian means (see Table 1). The results are benchmarked against the state-of-the-art CCA-based methods.

a particular advantage as it allows the assessment of a family of divergences by changing the value of α. It is possible to find the optimal value of α with respect to the classification accuracy for the investigated dataset. Through a cross-validation process with -1 ≤ α ≤ 1 (see Fig. 9), α = 0.6 is the optimal α for the current study. This procedure lasted 225.42 seconds and makes α-divergence the most costly method. This assessment goes through some of the divergences that are known to be particular cases of the α-divergence: the left-and right-sided Kullback-Leibler at α = -1 and α = 1 respectively, and the Bhattacharyya and S-divergence at α = 0. The results of α-divergence reported on Fig 6 is obtained with the optimized value of α = 0.6, and represents the best accuracy achievable with the family of α-divergence. The α-divergence yields the highest accuracy across distances and divergence considered in the current study. Unlike other distances and divergences, the α-divergence is not parameter-free and require an optimization step to determine the value of α. This is a limitation of the α-divergences for real BCI application, as it needs an extra processing step which is computationally demanding.

Despite the apparent good properties of the Wasserstein mean (see Fig. 7), the Wasserstein distance has the lowest accuracy among all Riemannian methods. While the Bhattacharyya distance and the S-divergence have different expressions in Table 1, they obtain similar results in terms of accuracy, and differ only in their computational time. For this reason, the Fig 6 only reports one of them.

The best trade-off regarding computation time and accuracy is provided by the Jeffreys divergence and mean, which achieves results close to the accuracy of the α-divergence or affine-invariance distance with a com-putation time that is an order of magnitude lower. The reported accuracy show that the Jeffreys divergence is a better estimator for classifiers than affine-invariant distance for all subjects but two.

To analyze the statistical significance of the improvement brought by Riemannian distances and divergences, we perform a paired t-test between pairs of results and we evaluate the p-values. The analysis reveals that the difference in performance of the two Euclidean distances, arithmetic and harmonic, are not statistically significant, with p-value p > 0.05. Also, despite the improvement in accuracy from state-of-the-art method, the test reveals that the improvement is not statistically significant (p > 0.05). However, it reveals a significant improvement between the Euclidean distances and the Riemannian ones, with p-values p < 0.05 (in the order of 10 -5 ). With the exception Wasserstein distance, the analysis shows no significant difference between the different Riemannian approaches. The most important point is that all the Riemannian approaches, with the exception of the Wasserstein distance, significantly improve the state-of-the-art (p < 0.05), that is MDM with any Riemannian distance or divergence outperforms the CCA-based methods for classifying SSVEP.

Online classification

Fig. 8 indicates the performances of online MDM with the different distances and divergences, using the K = 4 classes, including the resting-state class. Firstly, we observe that the online results are lower than the offline ones displayed on Fig. 6. There are two reasons: (i) Fig. 6 displays the results for only K = 3 classes, without the resting-state class; (ii) offline implementations capture more information than online ones (Chevallier Fig. 7 The TraDe plot displays the log-trace as a function of the log-determinant of the different covariance matrices: the markers indicate the type of the mean and the color map shows the distribution of covariance matrices used in the computation of the means. Most of the means align diagonally, away from the individual matrices. A swelling effect is visible on the determinant for the arithmetic mean and a shrinking effect for the harmonic mean. A shrinking effect on the trace could also be seen for the harmonic mean and the AIR mean. et al. 2018), leading to better accuracy. Secondly, online and asynchronous MDM with Riemannian metrics outperform those with Euclidean metrics. Wasserstein distance has the lowest accuracy among all Riemannian methods. These two observations are qualitatively similar to the observation made for the offline implementation.

The paired t-test shows that there is no statistically significant difference between MDM approaches with Riemannian metrics, apart from the Wasserstein that obtains significantly lower performances. It is to be noted that CCA-based methods, used for comparison in offline classification, are not adequate for online and asynchronous applications as they cannot separate SSVEP classes from the resting class. This limitation disqualifies them for any real implementation of BCI.

With these experimental results on real EEG, the benefit of dedicated Riemannian tools for processing covariance matrices is clearly demonstrated. As seen on offline experiments the Jeffreys divergence is a good choice, as it is a good trade-off between the accuracy and the computation time.

Conclusion

The objective of this study is to investigate how fast and accurate systems that could be applied in a loosely controlled environment, instead of carefully tuned systems that are designed in strongly constrained environ- ments. Riemannian approaches have been successfully applied on EEG signals for brain-computer interfaces, demonstrating their robustness to noise and to the intrinsic variation of EEG signals. A direct approach, like Minimum Distance to Riemannian Mean, yields comparable results with state-of-the-art without requiring a complex choice of parameters or a lengthy optimization phase. Covariance matrices belong to a curved space and, expressed in a Riemannian framework, several distances are defined in the literature. Most of these distances embed desirable invariances, avoiding costly calibration of spatial filters for each session and yielding estimation that are robust to artifacts and outliers. Different distances and divergences, along with their implementation are reviewed and investigated, including the Log-Euclidean distance, the affine-invariant Riemannian distance, the Kullback-Leibler divergence, the Jeffreys divergence, the S-divergence, the α-divergence and the Wasserstein distance.

This study shows that, for SSVEP classification, Riemannian approaches outperform previous state-ofthe-art methods based on CCA. Riemannian approaches also address the issue of resting class, where the user is not interacting with the system, which is a challenging task in machine learning. A thorough analysis of computational cost and of numerical stability is proposed in the experimental section. The best accuracy is obtained with the α-divergence, after finding the optimal value of α. The Jeffreys divergence offers the best accuracy for the lowest computational cost. At last, to ensure the reproducibility of the results and to propose a framework for evaluation of new approaches to the communities of brain-computer interfaces and information geometry, all the data, including the algorithms to process them are available under an open source license.

The presented Riemannian framework can be extended to MI or ERP paradigms (Congedo et al. 2017a). Future works on this topic could consist of investigating whether distances and divergences behave similarly with different paradigms, in terms of classification accuracy.

On a methodological level, these first results call for an extension of covariance-based features to other VEP paradigms: phase-based SSVEP and code-modulated VEP (c-VEP). These paradigms achieve very high information transfer rates and could greatly benefit from the accuracy improvement provided by Riemannian geometry.

Finally, another opportunity ahead is to evaluate the possibility of making transfer learning between frequencies for a given subject. This approach would allow to calibrate the system on single frequency and generalize the classifier's training to a large number of frequen-cies. This is particularly useful as it could help reduce the calibration time and allow subjects to keep their concentration for the real task.

Fig. 1

 1 Fig.1SSVEP-based BCI setup for controlling an arm exoskeleton on an electric-powered wheelchair. The SSVEP stimulation board is made of 3 groups of LEDs that blink at different frequencies.

Fig. 2

 2 Fig. 2 Samples of EEG trials from Eq. (5). For each sub-band X freq f only the channel Oz is shown. Each subplot shows the first second of a trial from classes: resting-state (a), 13 Hz (b), 21 Hz (c), 17 Hz (d).

  17-20) in (Cherian et al. 2011) (Sra 2016; Cherian et al. 2011) α-divergence D α LD from Eq. (3) in (Chebbi and Moakher 2012) Algorithm 1 in (Chebbi and Moakher 2012) (Chebbi and Moakher 2012)

Fig. 3

 3 Fig. 3 Illustration of the MDM principle for a biclass problem in a curved space.

Fig. 4

 4 Fig. 4 Representation of covariance matrices: each image is the covariance matrix mean Σk of the class k, for one session of the recording. The diagonal blocks show the covariance in different frequency bands, i.e. 13 Hz in the upper-left block, 21 Hz in the middle, and 17 Hz in the bottom-right.

Fig. 5

 5 Fig. 5 Processing time and accuracy (mean and standard deviation across subjects and trials) for all distances and divergences. Computation of averages during offline training (a) and of distances/divergences during classification (b). For computation times, methods were run on Matlab 7.11.0.584 (R2010b), on a Linux Intel i7 machine.

Fig. 8

 8 Fig.8Online classification accuracy (in %) for K = 4 classes, by subject. Classification is performed with MDM using either Euclidean or Riemannian means (see Table1).

Fig. 9

 9 Fig.8Online classification accuracy (in %) for K = 4 classes, by subject. Classification is performed with MDM using either Euclidean or Riemannian means (see Table1).

Table 1

 1 Formula of distances, divergences and means considered in this article.

	Distance/Divergence	Mean	References
	Arithmetic

  [START_REF] Álvarez-Esteban | A fixed-point approach to barycenters in Wasserstein space[END_REF])[START_REF] Agueh | Barycenters in the Wasserstein space[END_REF][START_REF] Álvarez-Esteban | A fixed-point approach to barycenters in Wasserstein space[END_REF]

Table 2

 2 Properties of distances and divergences of SPD matrices.

	Distance/Divergence Restriction on eigenvalues	Invariances
	Arithmetic	None	Rotation
	Harmonic	None	Rotation
	Log-Euclidean	Strictly positive	Rotation, inversion, similarity
	Affine-invariant	Strictly positive	Rotation, inversion, similarity, affine transform
	Kullback-Leibler	Strictly positive	Rotation, affine transform
	Jeffreys	Strictly positive	Rotation, inversion, affine transform
	S-divergence	Strictly positive	Rotation, inversion, affine transform
	α-divergence	Strictly positive	Affine transform
	Wasserstein	Positive	Rotation

Acknowledgements This work has been done without any support of the ANR or the ERC.

Compliance with Ethical Standards

At the time this study was carried out, Emmanuel Kalunga was an employee of VASTech company, Stellenbosch, South Africa; and Quentin Barthélemy was an employee of Mensia Technologies S.A., Paris, France. Sylvain Chevallier and Eric Monacelli declare that they have no conflict of interest. This study was conducted by the authors without any specific funding or grant.

All procedures performed in studies involving human participants 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.

Information Sharing Statement

All the data used in this paper could be freely downloaded from [START_REF] Chevallier | SSVEP data[END_REF] or using the MOABB framework [START_REF] Barachant | MOABB framework[END_REF]. The code for the experiments of this article is written in Matlab: the code is available on Kalunga (2018) for offline analysis, and on Kalunga (2015) for online analysis. A code written in Python is available on Chevallier (2020) for offline analysis, and on Bertrand-Lalo (2020) for online analysis.


Arnaud F Klein

Ludovic Arandel

Joëlle Marie

Denis Furling

Joelle Marie

Joelle Marie

FISH protocol for Myotonic Dystrophy type 1 cells

Keywords: myotonic dystrophy, nuclear RNA foci, CUG expansion, DMPK transcripts, MBNL1, fluorescent in situ hybridization, immunofluorescence

come

Introduction

Myotonic dystrophy type 1 (DM1 or Steinert's disease) is a multisystemic disorder characterized by myotonia, progressive muscle weakness and atrophy, cardiac conduction defects and cognitive impairment together with other endocrine dysfunctions (1, 2). The disease is caused by an abnormal expansion of a CTG repeated sequence within the 3'untranslated region (UTR) of the Dystrophia Myotonica Protein Kinase (DMPK) gene (3-5). The size of the CTG tract is lower than 37 repeats in healthy population, but can reach several thousands in DM1 patients. Moreover, CTG expansions in patients are instable at both the somatic and intergenerational levels (6). Since the identification of this expanded microsatellite repeats, DM1 has been described as the archetype of RNA gain-of-function diseases. The mutant DMPK allele is transcribed, but messenger RNAs containing expanded CUG repeats (CUGexp-DMPK) are retained into the nucleus and form discrete aggregates or foci (7, 8). These ribonuclear foci sequester factors of the Muscleblind-like (MBNL) RNA binding protein family (9-11), resulting in their functional loss and thereby the disruption of RNA processing events including the misregulation of the alternative splicing of many effector genes (12-15). More than forty misregulated splicing events have been described in the skeletal muscle of DM1 patients and some of them were associated with symptoms (12, 16-20).

The finding of ribonuclear aggregates or RNA foci in DM1 cells was a major breakthrough in the understanding of the molecular disease process (7) and their presence in affected tissues of DM1 patients is currently a hallmark of the disease (11,[START_REF] Jiang | Myotonic dystrophy type 1 is associated with nuclear foci of mutant RNA, sequestration of muscleblind proteins and deregulated alternative splicing in neurons[END_REF][START_REF] Michel | Sense and Antisense DMPK RNA Foci Accumulate in DM1 Tissues during Development[END_REF]. Nuclear CUGexp-RNA foci were also detected in both in vivo and in vitro models of the disease, and their presence, number or intensity is often used as a readout to assess efficacy of therapeutic approaches (23). RNA foci are composed of mutant CUGexp-DMPK transcripts that are normally spliced, capped and polyadenylated (8). However the abnormal binding of MBNL proteins to expanded CUG repeats contributes to the nuclear retention and/or impaired cytoplasmic export of these pathologic RNA (24-26).

In this chapter, we describe in details how to perform a fluorescent in situ hybridization (FISH) staining on human DM1 or WT cells to visualize CUGexp-RNA and DMPK transcripts. We have adapted protocol from previous works using fluorescent probes for either FISH or single molecule inexpensive FISH (smiFISH). We also describe the preparation and the fixation of the cells prior to hybridization, and an immunofluorescence protocol to detect MBNL1 in colocalization analysis purpose. The culture conditions for myoblasts and fibroblasts are adapted from Arandel et supplemented with 20% FBS and 50μg/ml gentamycin. For immortalized human myoblasts: mix of 199:DMEM (1:4 ratio) supplemented with 20% FBS, 50μg/ml gentamycin, 25 μg/ml fetuin, 0.5 ng/ml bFGF, 5 ng/ml EGF and 0.2 μg/ml dexamethasone. For fibroblast and immortalized fibroblasts: DMEM supplemented with 15% FBS and 50μg/ml gentamycin. 7. Differentiation medium (DM): For primary and immortalized myoblasts: DMEM supplemented with 5 μg/ml insulin. For trans-differentiated MYOD1 fibroblasts: DMEM supplemented with 5 μg/ml insulin and 4 μg/ml of doxycycline. -100 mM NaCl -300 mM Sucrose -10mM PIPES, pH 6.8 -3 mM MgCl 2 10. Blocking buffer: PBS, tween 0.1%, 3% BSA 11. Hoechst solution: 2mg/mL of Hoechst in PBS-T

Methods

Preparation of cells for hybridization and immunofluorescence

DM1 cells (myoblasts or fibroblasts) are grown on glass coverslips coated with gelatin:

Preparation of coverslips and coating

1. Place a glass coverslip in the bottom of a multi-well plate or a petri dish using forceps. 2. Add a gelatin solution to cover the coverslip (ex. 1ml/35mm petri) and incubated for 1h at RT. 3. Rinse twice with PBS. 4. Plate the cells at the required density (see Note 2).

Cell growth and fixation

1. Cells are grown at 37°C and 5%CO 2 . 2. Let the cells grow to the desired confluence (80% for proliferation or 100% for differentiation) (see Note 3). 3. Remove medium (see Note 4). 4. Carefully rinse the cells once with PBS. 5. Fix the cells with a solution of 4% PFA under chemical hood for 10 min (see Note 1). 6. Carefully rinse the cells twice with PBS. 7. Rinse the cells with EtOH 70%. 8. Incubate the cells overnight in EtOH 70% at 4°C (For convenience, cells/coversplips can be transferred to another multi-well plate using forceps).

Fluorescent in situ hybridization

Before starting:

1. See Note 5 for the standard incubation temperature. 2. Do not forget to work with RNAse-free solutions.

3. Turn on a hybridization oven at 37°C and another at 45°C. 4. Warm 50 mL of PBS-T at 45°C. 5. Prepare the wet chamber for incubation/hybridization.

Hybridization:

1. Carefully wash twice the cells with PBS-MgCl2 for 10 min. 2. Under the chemical hood: Prepare the hybridization buffer (ex. consider 20µl per 12mm coverslip, 50µL for 30mm)(see Note 1). 3. Place a piece of parafilm on the bottom of the chamber.

4. Dispose a drop of hybridization buffer on the parafilm . 5. Place the coverslip over the drop with the cells down. 6. Seal the incubation chamber with parafilm to avoid formamide fumes. 7. Incubate the cells in the chamber at 37°C for 90 min in a hybridization oven, in the dark (see note 6). 8. Under the chemical hood, replace the cells/coverslip in a multi-well plate using forceps. 9. Wash the cells once with PBS-T for 5 min at RT. 10. Wash the cells one with warmed PBS-T in a hybridization oven at 45°C for 30 min. 11. Incubate the cells once with Hoechst solution to counterstain the nucleus for 10 min. 12. Wash the cells twice with PBS-T for 10 min. 13. Mount the cells with an anti-fading mounting media on slide and stock the cells at least one night before visualizing the cells at the microscope (Figure 1).

smiFISH

Single Molecule inexpensive FISH or smiFISH is a technique allowing the visualization and quantification of single mRNA molecules using gene-specific unlabelled probe set hybridized (in vitro) to an fluorochrome-labelled secondary detector oligonucleotide (FLAP)(Figure 2A).

Design and preparation of DMPK probe set

1. The probes used to target DMPK transcripts using smiFISH protocol were designed using Oligostan script included in Fish_Quant package available at https://bitbucket.org/muellerflorian/fish quant (29). The probes consist of 26 to 30 nucleotide gene-specific sequences coupled to a TTACACTCGGACCTCGTCGACATGCATT sequence, complementary to the FLAP secondary oligonucleotide (Table 1). 2. The 24 probes target the coding region of DMPK transcripts (Figure 2B). 3. The 24 oligonucleotide probes are resuspended TE buffer at a stock concentration of 100µM. 4. Prepare an equimolar mixture of the 24 oligonucleotide probes (concentration of each individual probe 4.16µM). 5. Dilute the probe mixture 5 times in TE buffer to get a 20µM probe solution (final concentration of each individual probe 0.83µM). 6. The probe-set is kept at -20°C for long-term storage.

Preparation of smiFISH probe

Figure 1 :

 1 Figure 1: RNA aggregates (or foci) of CUGexp-DMPK transcripts were detected by FISH using a Cy3-(CAG)7 (red) in the nucleus (Hoechst, blue) of DM1 muscle cells .

Figure 2 :

 2 Figure 2: single molecule inexpensive FISH. A) schematic representation of the smiFISH principle. An ensemble of 24 gene specific probe are hybridized in vitro to a fluorescent secondary probe via the FLAP sequence. The resulting smiFISH probes are subsequently hybridized to the targeted transcripts in cells (adapted form Tsanov et al. 2016). B) Schematic representation of the 24 oligonucleotides probes on the DMPK transcripts. C) Representative image of DMPK transcripts (red) observed using smiFISH technique in differentiated human WT muscle cells (nucleus in blue).

Figure 3 :

 3 Figure 3: Localisation of MBNL1 in WT and DM1 myoblasts. A) Visualization of MBNL1 (in green) using immunofluoresence technique in differentiated human WT muscle cells (nucleus in blue). B) Visualization of both mutated DMPK transcripts using FISH technique (CAG7 proble, red) and MBNL1 using immunofluorescence technique (MB1a antibody, green) in the nucleus of DM1 myoblast (Hoechst , blue)

al. 2017 (27), the

  FISH protocol is adapted and optimized from[START_REF] Taneja | Foci of trinucleotide repeat transcripts in nuclei of myotonic dystrophy cells and tissues[END_REF][START_REF] Taneja | Localization of trinucleotide repeat sequences in myotonic dystrophy cells using a single fluorochrome-labeled PNA probe[END_REF]

7, 28), the smiFISH protocol from Tsanov et al 2016 (29); and

  MBNL1 immunofluorescence protocol from Arandel et al., and Byron et al.

(27, 30). 2. Materials 2.1. Preparation of cells for hybridization and immunofluorescence

	1. Petri dish or multi-well plates for cell culture
	2. Sterile glass coverslips (12, 14, 18 or 30 mm diameter)
	3. Microscope forceps
	4. Autoclaved solution of gelatin 0.5% (in water)
	5. Sterile PBS for cell culture
	6.

Growth medium (GM): For primary human myoblasts: mix of 199:DMEM (1:4 ratio)

Table 1 : sequence and modification of the probes used in FISH and smiFISH protocols

 1

	Protocol Name	Sequence 5' -> 3'
	FISH	CAG sequence	Cy3-CAGCAGCAGCAGCAGCAGCAG
	smiFISH DMPK-Probe 1	TCCGGGGGACACAGCAACCGCTGAATTTACACTCGGACCTCGTCGACATGCATT
		DMPK-Probe 2	ACAGGTAGTTCTCATCCTGGAAGGCGAAGTTTACACTCGGACCTCGTCGACATGCATT
		DMPK-Probe 3	CAGTCGGACCTCCTTAAGCCTCACCATTACACTCGGACCTCGTCGACATGCATT
		DMPK-Probe 4	CAACTGCCGGACGTGTGCCTCTAGGTTTACACTCGGACCTCGTCGACATGCATT
		DMPK-Probe 5	TTCTGGTTGTCCGTGCGGATGGCCTCTTACACTCGGACCTCGTCGACATGCATT
		DMPK-Probe 6	CGCTGGAACTGCCACTTCAGCTGTTTTTACACTCGGACCTCGTCGACATGCATT
		DMPK-Probe 7	CTGTCCCTGAGGGCCATGCAGGAGTATTACACTCGGACCTCGTCGACATGCATT
		DMPK-Probe 8	GTGAGCCCGTCCTCCACCAAGTCGAATTACACTCGGACCTCGTCGACATGCATT
		DMPK-Probe 9	TGCATGTGTCGGTGGCACCTTCGAAATCTTACACTCGGACCTCGTCGACATGCATT
		DMPK-Probe 10	GCTCCTTGTAGTGGACGATCTTGCCATTTACACTCGGACCTCGTCGACATGCATT

Acknowledgments :

This research was supported by Sorbonne Université, Institut National de la Santé et de la Recherche Médicale (INSERM), Centre National de la Recherche Scientifique (CNRS) and the Association Institut de Myologie (AIM).