P L C Rodrigues

C Jutten

M Congedo

Abstract-Objective: This paper presents a Transfer Learning approach for dealing with the statistical variability of EEG signals recorded on different sessions and/or from different subjects. This is a common problem faced by Brain-Computer Interfaces (BCI) and poses a challenge for systems that try to reuse data from previous recordings to avoid a calibration phase for new users or new sessions for the same user. Method: We propose a method based on Procrustes analysis for matching the statistical distributions of two datasets using simple geometrical transformations (translation, scaling and rotation) over the data points. We use symmetric positive definite matrices (SPD) as statistical features for describing the EEG signals, so the geometrical operations on the data points respect the intrinsic geometry of the SPD manifold. Because of its geometry-aware nature, we call our method the Riemannian Procrustes Analysis (RPA). We assess the improvement in Transfer Learning via RPA by performing classification tasks on simulated data and on eight publicly available BCI datasets covering three experimental paradigms (243 subjects in total). Results: Our results show that the classification accuracy with RPA is superior in comparison to other geometry-aware methods proposed in the literature. We also observe improvements in ensemble classification strategies when the statistics of the datasets are matched via RPA. Conclusion and significance: We present a simple yet powerful method for matching the statistical distributions of two datasets, thus paving the way to BCI systems capable of reusing data from previous sessions and avoid the need of a calibration procedure.

Index Terms-Brain-Computer Interface, Riemannian geometry, Transfer Learning, Covariance Matrices, EEG.

I. INTRODUCTION

A Brain-Computer Interface (BCI) is a system that allows a person to interact with a machine without any physical interaction. It works by extracting features from neuro-physiological signals (e.g., the power spectral densities on certain frequency bands) and assigning them to different classes. These classes may be associated to cognitive states, sensory responses, etc., and the features are chosen so that they are discriminative for each class. Among the many types of neural signals that have been used in BCI systems, electroencephalographic (EEG) recordings have received particular attention and are the focus of this paper. The success of EEG-based systems comes from several reasons, such as its low price as compared to other neuro-physiological modalities, its non-invasiveness, and its high temporal resolution. However, they also have to cope with low signal-to-noise ratio, low spatial resolution [START_REF] Congedo | EEG Source Analysis[END_REF], and high cross-session and cross-subject variability [START_REF] Lotte | A review of classification algorithms for EEGbased brain-computer interfaces: a 10 year update[END_REF].

While many signal processing techniques have been developed to ameliorate the quality of EEG signals and improve its spatial resolution [START_REF] Congedo | A new generation of braincomputer interface based on riemannian geometry[END_REF], the cross-session and cross-subject variability problem have received much less attention. The standard approach in BCI consists of re-training a statistical classifier at the beginning of every experimental session with the help of a sequence of calibration trials, a procedure that can be time consuming and is clearly suboptimal, since it does not leverage any information from past experiments. Instead, "second-generation BCIs" initialize classifiers using Transfer Learning and update their parameters along sessions, avoiding the calibration step altogether [START_REF] Congedo | EEG Source Analysis[END_REF].

We present in this paper a Transfer Learning approach based on geometrical transformations for matching the statistical distributions coming from two different experimental sessions, a source and a target session. For simplicity of exposition, hereafter we use the term "session" generically, however source and target may refer to different sessions as well as different subjects. We perform simple linear transformations, such as translation, rotation, and scaling, on the data points of both sessions with the goal of making the shape of their statistical distributions as similar as possible. Once the distributions are matched, one can expect that a classifier optimized for the data of the source session will work well enough with the data from the target session.

The branch of Machine Learning that studies the effects of mismatches between statistical distributions is called Transfer Learning [START_REF] Pan | A survey on transfer learning[END_REF]. It has been of great interest in several domains besides BCI, such as in computer vision, where the statistics of the data may vary due to changes in lighting conditions and acquisition devices, or in speech processing systems, where the changes in background noise and the differences in speaker genders and voice tonalities may affect the statistics of the signals. In [START_REF] Shimodaira | Improving predictive inference under covariate shift by weighting the log-likelihood function[END_REF], the phenomenon responsible for the drift in statistical distributions of two datasets was termed covariate shift and modelled by assuming that the distributions of the data points can be different for the source and target datasets, but the conditional distributions of the labels are the same. Ref. [START_REF] Sugiyama | Direct importance estimation for covariate shift adaptation[END_REF] presented real data examples on BCI experiments and showed that the covariate shift describes well the changes in statistics for this kind of application. A recent attempt in the literature on Transfer Learning has been to apply the theory of Optimal Transport to determine which geometrical transformations one should perform to match two statistical distributions, as proposed in [START_REF] Courty | Optimal transport for domain adaptation[END_REF] and applied to the BCI P300 paradigm in [START_REF] Gayraud | Optimal Transport Applied to Transfer Learning For P300 Detection[END_REF]. This approach begins by first defining a function that models the cost of transporting a point at location x to a location y. Then, it solves an optimization problem that minimizes the total cost of moving each point of the source dataset and make its statistical distribution as similar as possible to the distribution of the target dataset. Note that this is a completely unsupervised procedure and does not rely on any assumption regarding the statistical distributions of the datasets.

Another geometrical approach [START_REF] Zanini | Transfer learning: a Riemannian geometry framework with applications to braincomputer interfaces[END_REF] matches the statistical distributions from the source and target sessions by means of a re-centering of their data points to the origin of the SPD space, the identity matrix. During the process of review of this paper, we came across the work in [START_REF] Yair | Parallel transport on the cone manifold of SPD matrices for domain adaptation[END_REF], which also proposes a Transfer Learning procedure for SPD matrices in the spirit of [START_REF] Zanini | Transfer learning: a Riemannian geometry framework with applications to braincomputer interfaces[END_REF], but with the main difference that the data points are re-centered to the midpoint between the geometric means of the source and target datasets Besides of distribution matching based on geometrical transformations, there has been mainly two other kinds of proposals for doing Transfer Learning in the BCI literature [START_REF] Lotte | A review of classification algorithms for EEGbased brain-computer interfaces: a 10 year update[END_REF]. One is based on the concept of ensemble classifiers [START_REF] Congedo | A new generation of braincomputer interface based on riemannian geometry[END_REF], [START_REF] Blankertz | Invariant common spatial patterns: Alleviating nonstationarities in brain-computer interfacing[END_REF]- [START_REF] Waytowich | Spectral transfer learning using information geometry for a user-independent brain-computer interface[END_REF], where the information from multiple source datasets are combined into a "global" classifier, which is then used to label the trials from any other target dataset. There has also been works using Bayesian models to describe the variability of the statistics on the source datasets and gather information from multiple datasets [START_REF] Jayaram | Transfer Learning in Brain-Computer Interfaces[END_REF]. A recent approach that builds upon such Bayesian methods are the works in [START_REF] Kindermans | True zero-training brain-computer interfacing ? an online study[END_REF] and [START_REF] Hübner | Unsupervised learning for brain-computer interfaces based on event-related potentials: Review and online comparison [research frontier[END_REF], which use a special form of the P300 experimental paradigm to do classification with no calibration. There are two main differences between these approaches and our proposal in this work. First, our approach is based on matching as much as possible the information from each source-target pair of datasets. Thus, it can be used in addition to an ensembling approach as in [12], which will combine the information from multiple matched-source subjects. Second, our approach is paradigmagnostic and does not rely on any special modification of the experimental setup where the EEG signals are collected (as opposed to [START_REF] Hübner | Unsupervised learning for brain-computer interfaces based on event-related potentials: Review and online comparison [research frontier[END_REF]), a feature that is appealing to a great number of practionners.

Our approach for distribution matching is based on the concept of Procrustes Analysis (PA) [START_REF] Gower | Procrustes Problems[END_REF], a tool often employed in statistical shape analysis [START_REF] Kendall | A survey of the statistical theory of shape[END_REF] with applications in text analysis [START_REF] Wang | Manifold alignment using procrustes analysis[END_REF], protein alignment [START_REF] Wang | Manifold alignment using procrustes analysis[END_REF], and many other fields. PA works by first selecting a set of pairs of landmark points from two different shapes and then performing geometrical transformations to get these landmarks as close as possible to each other. In our context, the shapes to be matched are actually point clouds consisting of high-dimensional statistical features describing EEG signals, and the landmarks are points that can be used to describe these statistical distributions (e.g. the mean, the farthest point from the mean, etc.). Because of its linear nature in the Euclidean case, and the fact that the operations are always global (the same rotation/translation/scaling is applied to all points each time), the space of transformations that one can cover using Procrustes Analysis does not include all possible transformations between the statistics of datasets. Nevertheless, the results that we have obtained on BCI applications indicate that the set of transformations that we apply are rich enough to model the covariate shift between experimental sessions.

A particularity of our work is that we use the spatial covariance matrices (SCM) of the EEG signals as statistical features to discriminate between classes. SCMs are symmetric positive definite matrices (SPD) and as such they are defined in a Riemannian manifold [START_REF] Bhatia | Positive definite matrices[END_REF] whose geometry is taken into consideration during the classification procedure. Such an approach has proved sucessful in BCI applications in recent years [START_REF] Barachant | Multiclass brain-computer interface classification by Riemannian geometry[END_REF]- [START_REF] Yger | Riemannian approaches in braincomputer interfaces: A review[END_REF] and has become part of the state of the art of the field [START_REF] Lotte | A review of classification algorithms for EEGbased brain-computer interfaces: a 10 year update[END_REF]. The geometric transformations that we use for matching two statistical distributions are done taking into account the intrinsic geometry of the SPD manifold as well. Because of its geometry-aware nature, we call our method Riemannian Procrustes Analysis (RPA).

RPA can be seen as an evolution of the aforementioned procedures [START_REF] Zanini | Transfer learning: a Riemannian geometry framework with applications to braincomputer interfaces[END_REF] and [START_REF] Yair | Parallel transport on the cone manifold of SPD matrices for domain adaptation[END_REF], with the re-centering step corresponding to the first of a series of geometrical transformations. Furthermore, the procedures in [START_REF] Zanini | Transfer learning: a Riemannian geometry framework with applications to braincomputer interfaces[END_REF] and [START_REF] Yair | Parallel transport on the cone manifold of SPD matrices for domain adaptation[END_REF] are completely unsupervised, since they do not use any information from the labels of the data points, whereas RPA benefits from the labels in the source session (which are all known in advance) as well as from (at least part of) the labels that become sequentially available in the target session trial after trial.

We compared the performance of the RPA procedure for Transfer Learning to that of other distribution-matching methods proposed in [START_REF] Zanini | Transfer learning: a Riemannian geometry framework with applications to braincomputer interfaces[END_REF], [START_REF] Yair | Parallel transport on the cone manifold of SPD matrices for domain adaptation[END_REF] and [START_REF] Courty | Optimal transport for domain adaptation[END_REF]. The results demonstrate that the RPA yields a superior classification accuracy in both simulated and real datasets. We also observe that, in general, RPA needs a very small amount of labeled trials from the target dataset to work well. Moreover, it is always superior (or at least equal) to the performance obtained with calibration, i.e., simply using the available trials on the target dataset and not do any transfer from the source.

The remainder of the paper goes as follows: Section II gives a mathematical formulation for the paradigm of Transfer Learning that we consider in this work. Section III introduces the tools for manipulating and classifying SPD matrices in their Riemannian manifold, and Section IV presents the method of Riemannian Procrustes Analysis. Sections V and VI discuss the results on simulated and real data, respectively. Lastly, Section VII concludes the paper.

II. PROBLEM FORMULATION

We consider two datasets, the source (S) and the target (T ) datasets. They are comprised of couples

S = (C i , y i ) for i = 1, . . . , K S , T = ( C i , y i ) for i = 1, . . . , K T , (1)
with C i , C i ∈ R n×n being data points, and y i , ỹi ∈ {1, . . . , L} their corresponding class labels; K S and K T are the number of trials in the source and target sessions respectively. In this paper, the data points in S and T are always symmetric positive definite (SPD) matrices and are used to parametrize EEG multivariate time series [START_REF] Congedo | Riemannian geometry for EEG-based brain-computer interfaces; a primer and a review[END_REF]. As mentioned before, SPD matrices are defined in a Riemannian manifold, so it is important that operations involving them respect the intrinsic geometry of this space (see Section III for more information).

Our distribution matching method considers the semisupervised Transfer Learning paradigm [START_REF] Pan | A survey on transfer learning[END_REF], where one has knowledge of all the labels from the source dataset and access to a small subset of labels from the target dataset. Put in mathematical terms, we assume knowledge of all the labels from the elements in S and of a small subset T ⊂ T with

T = T ∪ T u and T ∩ T u = ∅, (2)
where stands for labeled and u for unlabeled. We further assume that T has at least one example from each class. This setup describes well applications where a few calibration points from the target dataset can be used to guide the Transfer Learning procedure. Another relevant case is online algorithms, where labels are available sequentially and augment the T dataset after each time step. Our goal is to train a classifier that leverages the available information from S and T and has good accuracy on the classification of data points from T u .

III. BACKGROUND AND NOTATION

This section begins with a brief introduction to concepts of Riemannian geometry on SPD matrices. Then, we define the notion of statistical distributions of SPD matrices and review a simple way to do classification on datasets containing this type of matrix.

A. The Symmetric Positive Definite manifold

Let P(n) be the set of n × n symmetric positive definite (SPD) matrices, which is defined as

P(n) = C ∈ R n×n C T = C, x T Cx > 0, ∀x ∈ R n .
(3) Matrices in P(n) lie in a manifold [START_REF] Bhatia | Positive definite matrices[END_REF], a set of points with the property that the neighborhood of each C ∈ P(n) can be bijectively mapped to an Euclidean space, also known as its tangent space T C P(n). In particular, because P(n) is an open subspace of the set Sym(n) of symmetric matrices in R n×n , we can identify its tangent space as simply being Sym(n) [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]. If we endow every tangent space of a manifold with a metric that changes smoothly along its elements, we say that we have a Riemannian manifold [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]. In this case, fundamental geometric notions are naturally defined, such as geodesic (shortest curve joining two points), distance between two points (length of the geodesic connecting them), the center of mass of a set of points, etc.

There are several possible choices of metric for P(n) and each one induces a different geometry that can be more or less adequate according to the applications we are interested in. A metric that is particularly relevant is the Affine-Invariant Riemannian metric (AIRM) [START_REF] Bhatia | Positive definite matrices[END_REF], defined for η, ξ ∈ T C P(n) as the inner product

η, ξ C = tr C -1 ηC -1 ξ , (4)
where C ∈ P(n) and tr(•) denotes the trace operator. It is well known [START_REF] Bhatia | Positive definite matrices[END_REF] that the distance between matrices C i , C j ∈ P(n) induced by ( 4) is

δ 2 R (C i , C j ) = n k=1 log 2 (λ k ) , (5)
where the λ k 's are the eigenvalues of

C -1 i C j or of similar matrix C -1/2 i C j C -1/2 i
. Because of the many interesting properties of ( 5), such as invariance to affine transformations by any invertible matrix A ∈ R n×n ,

δ 2 R (C i , C j ) = δ 2 R (AC i A T , AC j A T ) , (6)
invariance under inversion, etc. [START_REF] Bhatia | Positive definite matrices[END_REF], the AIRM has found great popularity in geometry-aware algorithms for processing SPD matrices [21] [25]. From now on, whenever we refer to P(n) we will be implicitly considering it has been equipped with the AIRM.

The center of mass according to distance (5) of a set of SPD matrices {C 1 , . . . , C K } is defined as [START_REF] Congedo | EEG Source Analysis[END_REF]

G {C i } K i=1 = argmin X∈P(n) K i=1 δ 2 R (X, C i ) . (7)
In words, G {C i } K i=1 is the point in the manifold minimizing the dispersion (variance) of the set of matrices. Note that when the C k 's are actually strictly positive scalars, G {C i } K i=1 is their geometric mean. This explains why many researchers adopt the term "geometric mean" to refer to the center of mass of a set of SPD matrices. For three or more matrices, there is no closed form solution for [START_REF] Courty | Optimal transport for domain adaptation[END_REF], so one has to resort to iterative algorithms as the one given in [START_REF] Congedo | Fixed point algorithms for estimating power means of positive definite matrices[END_REF]. The above definitions suffice for the intents of this paper, but the interested reader will find a thorough treatment of the subject in [START_REF] Bhatia | Positive definite matrices[END_REF].

B. Statistical distributions in the SPD manifold

In this work, we assume that the data points come from statistical distributions that can be parametrized just by their geometric mean and the dispersion of points around it. More precisely, we assume that the statistical distribution generating the samples of the dataset is a mixture of Riemannian Gaussian distributions on the SPD manifold (one for each class) [START_REF] Said | Riemannian Gaussian distributions on the space of symmetric positive definite matrices[END_REF].

Under this assumption, we parametrize the statistics of each dataset defined in (1) using a set consisting of L+2 elements : the geometric mean M of the dataset, the geometric means M k of each of the L classes, and the dispersion d around M . We have then

Θ S = M, M 1 , . . . , M L , d , Θ T = M , M 1 , . . . , M L , d , where M = G({C i | C i ∈ S}) , M = G({ C i | C i ∈ T })
are the center of mass or geometric means of the datasets (all classes combined),

d = Ci∈S δ 2 R (M, C i ) , d = Ci∈T δ 2 R ( M , C i )
are the dispersions around the geometric mean, and

M k = G({C i | C i ∈ S and y i = k}) , M k = G({ C i | C i ∈ T and y i = k})
are the geometric mean of the trials belonging to each class. Note that the parameters for the target dataset are estimated using only the data points from T . The law of great numbers from statistics also applies to datasets containing SPD matrices [START_REF] Bhatia | Positive definite matrices[END_REF]. In particular, if the elements of a dataset come from a statistical distribution with geometric mean M , the center of mass of a set of K matrices will converge to M as K grows. This implies that in experimental paradigms where the trials come in sequentially, it is reasonable to expect that with more trials one will obtain better and better estimates of the geometric mean.

C. Classification in the SPD manifold

There are many known classifiers for datasets whose elements are SPD matrices [START_REF] Congedo | Riemannian geometry for EEG-based brain-computer interfaces; a primer and a review[END_REF]. In this paper, we use the Minimum-Distance to Mean (MDM) classifier [START_REF] Barachant | Multiclass brain-computer interface classification by Riemannian geometry[END_REF], which works as follows: in the training phase, calculate the geometric means for each class of a training dataset (D train ). In the testing phase, each matrix C i from the testing dataset (D test ) is associated to the label y i of the class mean that is the closest to C i . Note that the statistical model described in Section III-B fits well the MDM classifier, since it assumes that the statistics of the dataset can be described by the geometric means of its classes.

IV. TRANSFER LEARNING VIA PROCRUSTES ANALYSIS

This section presents the Riemannian version of the Procrustes Analysis or RPA. We first introduce the concept of Procrustes Analysis on an Euclidean setting and then describe how to perform equivalent transformations on the SPD manifold. Then, we justify such operations with the help of a model relating the statistical distributions of the source and target datasets. Such model is the main theoretical contribution of this paper, since it gives a concrete justification for the geometric operations done in the RPA procedure and allows for a better comprehension of the assumptions that one has to make regarding the statistical distributions of the datasets.

A. Matching statistical distributions

A common approach to match two shapes in Euclidean space is the Procrustes Analysis [START_REF] Gower | Procrustes Problems[END_REF], which works as follows : suppose we have two sets of landmark points

X = x i ∈ R n m i=1 and X = x i ∈ R n m i=1 , (8)
and assume there is a linear relationship relating the m pairs of landmark points as in

x i -m = s U x i -m , (9)
where s ∈ R, m, m ∈ R n , and U ∈ R n×n is an orthogonal matrix. The goal of the procedure is to determine the values of {s, m, m, U } so to obtain a new set X (PA) containing points x(PA) i that matches exactly with x i , where

x(PA) i -m = 1 s U T x i -m . ( 10
)
Note that the operations transforming x i can be interpreted as a re-centering to zero (subtracting m) followed by a stretching or compression (division by s), and a rotation (multiplication by U T ); the final re-centering to m is optional, since it is often more interesting to re-center X to the origin and consider only the zero-mean matched shapes.

B. Riemannian Procrustes analysis (RPA)

In order to do Procrustes Analysis on SPD matrices, we have to adapt the steps of re-centering, stretching and rotation according to the intrinsic geometry of P(n). We call such procedure Riemannian Procrustes Analysis (RPA) and describe its steps here below :

1) Re-center matrices to identity: In P(n) the Identity matrix plays the role of the origin of the space. Therefore, the first step of RPA is to transform the matrices in S and T so they are both centered around I n (see Fig. 1B). This amounts to the transformation proposed in [START_REF] Zanini | Transfer learning: a Riemannian geometry framework with applications to braincomputer interfaces[END_REF] if the covariance matrices used to describe the resting activity of each session were chosen to be the geometric mean of the trials of each dataset.

Due to the affine-invariance of ( 5) and ( 7), the geometric mean of the set of re-centered matrices

C (rct) i = M -1/2 C i M -1/2 (11)
is I n . Moreover, since M is estimated from a subset of points in T ⊂ T , the geometric mean of the set of matrices

C (rct) i = M -1/2 C i M -1/2 (12)
is approximately the identity matrix (it tends to the identity as the number of elements in T grows).

We have then two new datasets consisting of re-centered matrices

S (rct) = (C (rct) i , y i ) for i = 1, . . . , K S , T (rct) = ( C (rct) i , y i ) for i = 1, . . . , K T , (13)
with the indices of the partition T (rct) = T (rct) ∪ T (rct) u being the same as in [START_REF] Congedo | EEG Source Analysis[END_REF].

2) Equalize the dispersions on each dataset: The next step of RPA consists of rescaling the distributions on both datasets so that their dispersions around the mean are the same (see Fig. 1C). To do so, we can see from (5) that

δ 2 R C (rct) i s , I n = s 2 δ 2 R C (rct) i , I n , (14)
which implies that one can modulate the dispersion of T (rct) by simply moving each of its matrices along the geodesic that links it to the identity matrix. Note that the parameter s plays the same role of the scaling factor in [START_REF] Zanini | Transfer learning: a Riemannian geometry framework with applications to braincomputer interfaces[END_REF]. We match the dispersions from source and target by building a new dataset T (str) containing the stretched matrices

C (str) i = C (rct) i s , (15)
where we require s ∈ R to verify

s 2 = d/ d . (16)
Note that the re-centering of matrices in Step 1 does not alter the dispersion of the matrices around their geometric mean, which means that the stretching step could have been done before re-centering the matrices in T . However, in this case, the geodesic move in [START_REF] Kindermans | True zero-training brain-computer interfacing ? an online study[END_REF] would have to be done with respect to M , that is, we would have to use a more involved relation

C (str) i = M 1/2 M -1/2 C i M -1/2 s M 1/2 . ( 17
)
Note that up to this point no information from the trials' classes has been used. We say then that the recentering and stretching operations form the unsupervised part of the RPA method.

3) Rotate around the geometric mean: The last step consists of rotating the matrices from T (str) around the origin and matching the orientation of its point cloud with that of S (rct) (see Fig. 1D). To do so, we note that if U is an orthogonal matrix, then

δ 2 R (U T C (str) i U, U T U ) = δ 2 R ( C (str) i , I n ) , (18)
which indicates that the effect of an orthogonal matrix over a set of matrices centered at the identity is that of a rotation around their mean. We form a new dataset T (rot) containing rotated matrices with

C (rot) i = U T C (str) i U , (19)
where U is an orthogonal matrix to be determined from the data. As we will see in the next sub-section, matrix U is determined using the labels from the trials, so it corresponds to the supervised part of the procedure.

C. The orthogonal matrix U

The procedure to determine matrix U comes up naturally once the assumptions of the RPA method are written in mathematical form. For simplicity, we will first assume that T = T . We use the geometric means of the source and target datasets as landmarks to be matched, so one can write

M = A M A T , (20)
and

M k = A M k A T k ∈ {1, . . . , L} , (21)
where M, M , M k , M k are all defined in Section III-B, and A ∈ R n×n is an unknown invertible matrix that models the discrepancies between the statistics of the source and target datasets. We can rewrite the relation in [START_REF] Bhatia | Positive definite matrices[END_REF] as

M 1/2 M 1/2 = A M 1/2 M 1/2 A T , (22)
I n = M -1/2 A M 1/2 M 1/2 A T M -1/2 , (23)
I n = M -1/2 AM 1/2 M -1/2 AM 1/2 T , (24)
U U T = M -1/2 AM 1/2 M -1/2 AM 1/2 T , ( 25
)
where U is the n × n orthogonal matrix that we want to determine. Matrix U can then be simply written as

U = M -1/2 AM 1/2 , ( 26
)
where M and M are directly estimated from the data points as explained in Section III-A, and A remains unknown. To determine an expression for U only in terms of variables that can be estimated from the data, we use ( 26) in ( 21) to get

M k = M 1/2 U M -1/2 M k M 1/2 U M -1/2 T , (27)
M -1/2 M k M -1/2 = U M -1/2 M k M -1/2 U T . (28)
Defining the matrices

G k = M -1/2 M k M -1/2 (29)
and

G k = M -1/2 M k M -1/2 , ( 30
)
and using the expression in [START_REF] Bhatia | Positive definite matrices[END_REF] for M , one can rewrite [START_REF] Bhatia | Procrustes problems in riemannian manifolds of positive definite matrices[END_REF] as

G k = M 1/2 A -T M -1/2 G k M 1/2 A T M -1/2 (31)
and conclude that G k and G k are related via a similarity transform. Similar matrices having the same eigenvalues, one can write the eigendecompositions

G k = Q k ΛQ T k and G k = Q k Λ Q T k so that (28) becomes Q k Λ Q T k = U Q k ΛQ T k U T . (32)
Solving [START_REF] Lafon | Diffusion maps and coarse-graining: a unified framework for dimensionality reduction, graph partitioning, and data set parameterization[END_REF] for U we obtain

U = Q k Q T k , (33)
which is ultimately an expression in terms of variables estimated from the dataset. Note that [START_REF] Rodrigues | Multivariate timeseries analysis via manifold learning[END_REF] is also the solution to the following optimization problem minimize

U T U =In δ 2 R G k , U G k U T , (34)
and so it can be interpreted as the orthogonal matrix that acts to minimize the distance between a modified version of the class means of the source and target datasets.

D. Determining U from data

In the previous section, we assumed T = T . In practice, however, we have T ⊂ T , so the estimation of the class means of the target dataset are only approximations of the real class means of the statistical distribution. Because of this, instead of giving preference to a particular noisy estimate of a class mean to determine U via [START_REF] Rodrigues | Multivariate timeseries analysis via manifold learning[END_REF], we obtain it as a solution to the following optimization problem on the manifold of orthogonal matrices : minimize

U T U =In L k=1 w k δ 2 R G k , U G k U T , (35)
where the w k ∈ [0, 1] are coefficients allowing to balance the optimization according to the quality of the estimators for the mean of each class. Note that problem ( 35) is a generalization of the Procrustes problem in the SPD manifold proposed in [START_REF] Bhatia | Procrustes problems in riemannian manifolds of positive definite matrices[END_REF]. We solve (35) using a special form of the steepestdescent algorithm adapted for optimization procedures on manifolds, as described in [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]. To do so, we first rewrite each term of the cost function in [START_REF] Jayaram | Moabb: Trustworthy algorithm benchmarking for bcis[END_REF] as

f k (U ) = δ 2 R G k , U G k U T (36)
and express its Jacobian as

D U L(U ) = L k=1 w k D U f k (U ) , (37)
with

D U f k (U ) = 4 log G k U G k U T U , (38)
where the derivative of the AIRM distance was obtained from [START_REF] Moakher | A differential geometric approach to the geometric mean of symmetric positive-definite matrices[END_REF]. On each iteration of the gradient descent procedure, the vector D U f k (U ) is projected onto the tangent space of the manifold of orthogonal matrices (see [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF] for details). We used the pymanopt package [START_REF] Townsend | Pymanopt: A python toolbox for optimization on manifolds using automatic differentiation[END_REF] for carrying out the optimization procedure.

E. Classification on the transformed datasets

Once the datasets have been transformed with RPA, we use the MDM classifier (see Section III-C) to infer the unknown labels from the elements in T u . In the training phase we have

D train = S (rct) ∪ ( C (rot) i , y i ) ( C i , y i ) ∈ T , (39)
and in the testing phase we infer the labels from

D test = ( C (rot) i , ỹi ) ( C i , ỹi ) ∈ T u . (40)
F. Summary of the RPA method Algorithm 1 recapitulates the steps of a classification task using RPA for matching the statistical distributions of the source and target datasets. 2) Output: accuracy of classification using MDM on Tu 1 Estimate M and M from the data in S and T 2 Re-center the matrices in S and T using ( 11) and (12), and form new datasets S (rct) and T (rct) = T (rct) ∪ T (rct) u 3 Calculate the ratio of dispersions in S (rct) and T (rct) as in [START_REF] Hübner | Unsupervised learning for brain-computer interfaces based on event-related potentials: Review and online comparison [research frontier[END_REF] and use it to form the new dataset T (str) = T (str) ∪ T (str) u with matrices as described in [START_REF] Kindermans | True zero-training brain-computer interfacing ? an online study[END_REF] 4 Estimate matrices M k and M k for k ∈ {1, . . . , L} and obtain the orthogonal matrix U as a solution from (35) 5 Rotate the matrices from T (str) as in ( 19) and obtain G. An interpretation of the steps in RPA

T (rot) = T (rot) ∪ T (rot)
We give now an interpretation of the steps of RPA in terms of the statistical distributions of the datasets. Without loss of generality, we will consider that d = d, since the dispersions can always be made equal prior to the transformations. We will also assume T = T for simplicity.

The relations in ( 20) and ( 21) define which landmark data points we should match in the RPA procedure, an approach that is justified from the fact that we parametrize the statistics of S and T using their geometric means, as described in Section III-B. From this observation, one can also conclude that a simple approach for matching the statistical distributions of S and T would be to estimate matrix A from the available data points and apply A -1 to all the elements of T , as in

C i → A -1 C i A -T . (41)
From ( 26) we can write

A = M 1/2 U M -1/2 , ( 42
)
where M and M are estimated directly from the dataset, and the orthogonal matrix U is determined as discussed in Section IV-C. Applying A -1 to the matrices in T , we get

A -1 C i A -T = M 1/2 U T M -1/2 C i M -1/2 U M 1/2 , (43)
which describes the same steps of RPA: re-center to identity, stretch (with s = 1) and rotate, followed by a translation of the mean back to M . From the expressions above, we see that the sequence of operations in RPA are nicely justified by the assumptions of our statistical model for the data points.

V. NUMERICAL ILLUSTRATIONS : SIMULATED DATA

A. The dataset

Our example on simulated data consists of a source and a target datasets containing 2×2 SPD matrices and belonging to two classes. Data points from the source dataset are generated as follows:

1) Generate a random SPD matrix M 1 ∈ P(2) and define it to be the geometric mean of class 1 2) Generate N t = 100 random SPD matrices around M 1 by mapping small random tangent vectors (norm fixed to ε = 10) from T M1 P(2) back to the SPD manifold. We associate to each of these matrices the label y i = 1 3) Generate a random SPD matrix M 2 whose distance to M 1 is s = 5. This is the geometric mean for class 2 4) Generate N t = 100 random SPD matrices around M 2 by mapping small random tangent vectors (norm fixed to ε = 10) from T M2 P(2) back to the SPD manifold. These matrices have a label y i = 2 associated to them We generate the data points from the target dataset (T ) exactly as for the source dataset, but add an extra translation step that ensures that the geometric mean of all the matrices from S will be at a distance d = 8 from the geometric mean of T .

B. Visualization of the steps of RPA

We use the algorithm of Diffusion maps [START_REF] Coifman | Diffusion maps[END_REF] to obtain new representations of our data points using only two axis. This nonlinear dimensionality reduction algorithm works by first building a matrix with all pairwise distances between data points. Then, it calculates a new matrix called Laplacian which contains important information on the geometry of the low-dimensional manifold where the data points are assumed to live. The axis for the new data representations are then obtained from the spectral decomposition of the Laplacian matrix (see [START_REF] Lafon | Diffusion maps and coarse-graining: a unified framework for dimensionality reduction, graph partitioning, and data set parameterization[END_REF] for more details on the algorithm and [START_REF] Rodrigues | Multivariate timeseries analysis via manifold learning[END_REF] for an application on multivariate time series).

Figure 1 illustrates the distribution of data points after each step of RPA applied to the source and target datasets. In this example, we consider that we know the labels of all matrices from the target dataset, i.e., T = T . Figure 1A shows the point clouds of each dataset, which are clearly unmatched. After recentering (Fig. 1B), stretching (Fig. 1C) and rotating (Fig. 1D), the statistical distributions get matched and the same classifier can be used on both datasets.

C. Classification accuracy after RPA

We compare the classification accuracy results on the simulated dataset for six different methods of Transfer Learning. In each of them, the training and testing datasets are different but the classifier is always the MDM :

• direct (DCT): direct use of the points from the source dataset to do classification on the target dataset (no transformation whatsoever). • recentering (RCT): transfer learning considering only the data points of each dataset recentered to I n . This corresponds to step (1) in the RPA procedure and is similar to what has been done in [START_REF] Zanini | Transfer learning: a Riemannian geometry framework with applications to braincomputer interfaces[END_REF]. • parallel transport (PRL): transfer learning using the method proposed in [START_REF] Yair | Parallel transport on the cone manifold of SPD matrices for domain adaptation[END_REF]. The procedure is analogous to RCT, but with the points being re-centered to the halfway point along the geodesic path linking the geometric means of each dataset instead of the Identity matrix. • optimal transport (OPT): transfer learning using the optimal transport approach proposed in [START_REF] Courty | Optimal transport for domain adaptation[END_REF] and adapted to take into account the fact that we have data points defined in the SPD manifold instead of Euclidean vectors. • RPA: transfer learning with matrices transformed using RPA, as described in Section IV-B. • calibration: classification using only the labeled trials available in the target dataset, with no help from the data in the source dataset. We assess the performance of each method via a randomized cross-validation procedure consisting of:

1) Select 2n random elements from T (n from each class). These data points define T 2) Define the test dataset T u containing the other 200 -2n elements of T 3) Obtain the accuracy of the classification via MDM for this particular partition of T 4) Repeat the three preceeding steps 10 times and get the mean accuracy for each method. The results in Fig. 2 show that the DCT pipeline gives classification results at chance level (0.5) independently of the number of matrices available in T . We also observe that simply using RCT already greatly improves classification accuracy, as reported in [START_REF] Zanini | Transfer learning: a Riemannian geometry framework with applications to braincomputer interfaces[END_REF]. Our RPA method further improves the results. We also observe that PRL has virtually the same performance as RCT, which is not surprising, since they are both unsupervised methods based on the idea of re-centering the datasets to a common point in the SPD manifold. The results with OPT are equivalent to RCT and PRL as well. The accuracy with calibration improves when the number of available labels in the target dataset increases, eventually converging to the same performance as RPA. This result is not surprising, because with a sufficient amount of data in T it is already possible to train a good classifier without the need of doing transfer learning.

Our observations in this session are in accordance with the theoretical results of [START_REF] Ben-David | A theory of learning from different domains[END_REF], which says that "if there is enough target data, then no source data are needed (...). This is because the possible reduction in error due to additional source data is always less than the increase in error caused by the source data being too far from the target data". Such a result points to the existence of a certain "saturation" effect in the quality of transfer learning when too many trials are available in the target session, a behavior that could be exploited to decide when to stop transferring information from previous experimental sessions.

VI. NUMERICAL ILLUSTRATIONS : REAL DATA

The analysis on real data were all done using EEG recordings from experiments with Brain-Computer Interfaces (BCI)

(A) (C) (B) (D)
original recentering stretching rotation Fig. 1. Representation of the sequence of operations of RPA applied to a dataset simulated as described in Section V-A (better visualized with colors). Each point on the scatter plot represents a SPD matrix and the axes for the figures were obtained using Diffusion Maps [START_REF] Lafon | Diffusion maps and coarse-graining: a unified framework for dimensionality reduction, graph partitioning, and data set parameterization[END_REF]. The filled dots (degree of transparency set to α = 0.30) represent the target dataset whereas the circles are the source dataset. Each color represents a class and the black star is the Identity matrix. (A) Distribution of the SPD matrices in the source and target datasets as they are originally available and (B) after re-centering their geometric means to the Identity. In (C) the distribution after the stretching operation and (D) after the rotation. experiments. As usual, the EEG signal is a n-dimensional multivariate time-series and is denoted by x(t) where each dimension represents an electrode. Each experimental trial i lasts a few seconds and is associated to a matrix X i ∈ R n×T , where T is the number of time samples defining the trial.

To every trial we associate a SPD matrix C i describing its multivariate statistics and a label y i indicating what was the task performed during the trial. The dataset for each subject is composed of a set of couples (C i , y i ). Our investigation focuses on the classification accuracy of a MDM classifier that is trained with the data from a source subject and is used to classify the signals from a target subject. We compare the performance of such classifier using the different Transfer Learning strategies described in Section V-C.

A. The datasets

The investigations were carried out on eight datasets covering three different BCI paradigms. All Motor Imagery datasets are publicly available and were downloaded and pre-processed using the MOABB framework [START_REF] Jayaram | Moabb: Trustworthy algorithm benchmarking for bcis[END_REF]. The P300 dataset comes from experiments performed in our laboratory on the P300based game Brain Invaders [START_REF] Congedo | Brain invaders: a prototype of an open-source p300-based video game working with the openvibe platform[END_REF]. The SSVEP dataset was the same as the one presented in [START_REF] Kalunga | Online SSVEP-based BCI using Riemannian geometry[END_REF]. See Table I for a brief overview of each dataset's features.

We estimated the SPD matrices for each BCI paradigm differently. For the MI datasets, the SPD matrices were the spatial covariance matrices of the multivariate EEG recordings. The signals of each trial in the SSVEP paradigm were filtered using bandpass filters around certain frequencies of interest and its SPD matrices were diagonal blocks concatenating the spatial covariance matrices of the filtered signals [START_REF] Congedo | EEG Source Analysis[END_REF]. For the P300, the SPD matrix of each trial was obtained using the approach from [START_REF] Barachant | A plug & play P300 BCI using information geometry[END_REF], where one estimates a special form of covariance matrix that captures the influence of event-related potentials in each trial.

B. Comparing cross-subject classification accuracies

In this section, we compare the accuracy of the classification of trials for each pair of target and source subject. The classification is done using the MDM classifier and the values of the accuracies are assessed using the same cross-validation scheme explained in Section V-C.

We begin with a qualitative analysis of the cross-subject classification accuracy using a tool from combinatorial data analysis called seriation [START_REF] Liiv | Seriation and matrix reordering methods: An historical overview[END_REF]. This procedure sorts the lines and columns of a data matrix in order to make relevant patterns appear. In our case, the matrix S to be re-ordered contains at its (i, j) coordinate the accuracy of the classification using subject i as target and subject j as source. The rows of S are then sorted in decreasing order of columns sum, and the columns of this row-sorted-matrix S are rearranged in decreasing order of their rows sum. The output of this procedure is a new representation where the pairs of source-target subjects with the best accuracy are located at the top-left region of the matrix, while the worst pairs are at the bottom-right region. Figure 3 shows the results of this seriation procedure on the PhysionetMI dataset for two sizes of T (the number of labeled target trials) and three different pipelines: DCT, RCT, and RPA. We observe that with RCT and RPA there are more pairs of subjects with high values of cross-subject classification than with DCT. In particular, for RCT and RPA we note that there are a few target subjects that have very good accuracy on classification for almost all possible source subjects, a feature that is possibly related to the performance of each target subject to classify its own trials (intra-subject accuracy). This can been interpreted as : subjects that are "good" for classifying their own data should be "good" for receiving information from other source subjects. We also observe a clear improvement in the average value of the cross-subject classification accuracies when more points are available in T .

Our next analysis consisted in calculating the mean over all cross-subject AUC (Area Under the ROC Curve) for each Transfer Learning pipeline on each dataset. We used these values as quantitative measures for assessing whether one pipeline is better than the other on cross-subject classification. The scores are shown in Table II. We should mention that only the subjects whose intra-score AUC (i.e., classification of its own data) was above chance level were used in these calculations.

Figure 4 shows the results of statistical tests performed on each pair of methods, allowing for a more substantiated assessment of the performance of the methods. The statistical tests comparing method A versus method B were carried out in the following way: (1) For each source subject j, we perform a signed paired t-test comparing the scores of method A to method B along all target subjects. Each of these tests yield a statistic T j and a p-value p j is obtained via permutations tests [START_REF] Edgington | Randomization Tests[END_REF]. [START_REF] Lotte | A review of classification algorithms for EEGbased brain-computer interfaces: a 10 year update[END_REF] We combine the p-values of all the source subjects using Stouffer's Z-score method [START_REF] Zaykin | Optimally weighted z-test is a powerful method for combining probabilities in meta-analysis[END_REF]. This yields a single p-value for the comparison between methods as well as the direction to which the null hypothesis has been rejected (i.e., whether method A is better than B or vice-versa).

(3) We adjust p-values of each pairwise comparison using Holm's step-down procedure [START_REF] Holm | A simple sequentially rejective multiple test procedure[END_REF] to account for the multiple comparison problem.

The results in Figure 4 indicate that when there are enough points in T ("enough" depending on each dataset), transforming the data points with RCT, PRL or RPA is always better than not doing any distribution matching (DCT). We also observe that most of the time there is no statistical significance between the results with PRL and RCT, as expected, since they both amount to re-centering the datasets to a new point in the SPD manifold. For increasing values of N (the number of labeled trials in the target dataset), RPA gets better in comparison to almost all other methods, as expected and observed in Figure 2 for simulated data. Interestingly, OPT has very poor results in comparison to all other methods, probably because it does not use any prior hypothesis on the statistical distributions of the datasets and has to solve a difficult optimization problem to determine its transportation plan. Lastly, during our statistical analysis of the results, we have observed that better results on Transfer Learning via RPA are often associated to good intra-subject accuracy, since in this case the estimation of the class means is more stable and thus the rotation matrix U is better estimated. This explains why for some databases (e.g. PhysionetMI) the RPA is not necessarily the best method for Transfer Learning and an unsupervised approach like RCT has better results.

C. The role of the size of T

As pointed out in the simulation results from Section V-C, when the size of T increases, using Transfer Learning is no longer relevant, since one may already have enough data to build a good classifier for the target subject. To investigate this behavior on our real datasets, we compared the cross-subject classification accuracy of RPA to that of the calibration method (see Section V-C for details).

Figure 5 shows a scatter plot comparing the classification accuracies on the MunichMI dataset. We see that when T grows, there are more pairs of subjects for which the classification using the calibration method on the target subject is better than doing transfer learning via RPA (28% to 34% of all the pairs of subjects). However, the location of the cloud of points in the figure indicates that the Transfer Learning with RPA is still superior to the Calibration method for most pairs of subjects. We used a one-sided paired t-test with random data permutation [START_REF] Edgington | Randomization Tests[END_REF] to compare the accuracies of RPA and calibration on each dataset for different sizes of the T . The null hypothesis of equivalency between the two methods was rejected (p < 0.01) on almost all tests, the only exception being for those on the BNCI2014001 dataset. Moreover, for All conclusions are with p < 0.05 corrected via Holm's adjustement [START_REF] Holm | A simple sequentially rejective multiple test procedure[END_REF]. For instance, we see that for dataset Cho2017, the method RPA is inferior to RCT when N = 1, but RPA becomes superior when N = 25. Furthermore, for this same dataset, there's no statistical difference in the comparison between PRL and RCT for any N . the tests where H 0 was rejected, we observed a superiority of RPA in comparison to calibration.

N 3 N 4 N 2 L > R L < R L = R (L) (R)

D. Combining information from multiple subjects

We investigated how the matching of statistical distributions via RPA affects the performance of two baseline methods for gathering information from the data of multiple subjects: pooling and ensembling. The classification for each target subject is done using information coming from all other source subjects available in the database. Following the same approach as in previous sections, we only considered source subjects featuring an intra-subject accuracy above chance level, i.e., subjects in which it is meaningful to use transfer learning. The experiments were done on the PhysionetMI dataset with |T | = 15 labeled trials available for each target subject and the Cho2017 dataset with |T | = 25.

The pooling strategy consists of gathering for each target subject the data from all other source subjects into one big dataset. Then, a classifier is trained on the pooled dataset and used to infer the trials from the target subject. We compared the performance of a MDM classifier when the source subjects were pooled with no transformation (DCT) to when the statistical distributions of each source subject were matched to that from the target subject using RCT or RPA (PRL is not fit for pooling, since the matrices are not all recentered to the same place in the SPD manifold). The boxplots in Figure 6 show the distributions of the classification scores of each of the target subjects. Using pairwise one-sided paired t-test with random permutations, the null hypothesis of equivalency between the scores of DCT, RCT, and RPA were all rejected Fig. 6. Box plots with the distribution of the classification scores for the ensembling and pooling strategies for different methods of statistical matching between datasets. For the Cho2017 dataset we had 25 labeled trials in the target dataset and for the PhysionetMI there were 15 labeled trials in it. with p < 10 -6 (adjusted for multiple comparisons). The results show a clear improvement in the average score for the pooling strategy when using a method for matching the statistics of the source and target datasets, with differences of at least 15% between RPA and DCT for both datasets.

Our second analysis considered an ensembling strategy, where the trials of each target subject were classified using a majority voting scheme. These votes came from MDM classifiers trained on all other source subjects and were weighted equally. The results in Figure 6 show the scores with each method (including the PRL approach this time). To compare the scores of each method, we used pairwise one-sided paired t-tests with random permutations (corrected for multiple comparisons). The results of the statistical tests indicate that the ensembling strategy with RPA is superior as compared to DCT in the PhysionetMI dataset (p < 0.05) but they are equivalent for the Cho2017 dataset (p = 0.23). The RCT method is superior to DCT for both datasets (p < 0.01) whereas the scores with PRL are equivalent to DCT for both datasets. We see then that the ensembling strategy can also be improved when adding an extra step for matching the statistics of the datasets of each pair of source-target subjects.

VII. CONCLUSION AND FUTURE WORKS

In this paper, we have presented a new method for overcoming the negative effects of statistical distribution mismatch between subjects in BCI classification. Our proposal consists of a sequence of geometrical transformations on the data points from two sessions with the intention of making the shapes of their point clouds in a high-dimensional space similar as possible. The inspiration for this method comes from Procrustes Analysis, however, here the method has been adapted to the case where the data points are SPD matrices and live in a Riemannian manifold, which is here proposed for the first time. Another relevant theoretical contribution is the mathematical framework proposed in Section IV-B, which includes the methods in [START_REF] Zanini | Transfer learning: a Riemannian geometry framework with applications to braincomputer interfaces[END_REF] and [START_REF] Yair | Parallel transport on the cone manifold of SPD matrices for domain adaptation[END_REF] and extends them, leading to our RPA method. Such formalism allows for a better understanding of the intrinsic assumptions regarding the statistics of the data points during the distribution matching procedure.

An important aspect of our proposal is that we exploit the availability of supervised information in the source session as well as the sequential nature of the trials in the target session. It should be noted, however, that when no labels are available for the target session, a re-centering of data points based solely on the geometric means of each dataset (which does not rely on any supervised information) already greatly improves the cross-session and cross-subject classification, as first noted in [START_REF] Zanini | Transfer learning: a Riemannian geometry framework with applications to braincomputer interfaces[END_REF]. This would be the case, for instance, in BCI applications for people with extreme motor disability, where the labeling of classes is very challenging. In this kind of situation, one may still perform the re-centering and stretching steps of the RPA method for matching the statistical distributions, turning the transfer learning procedure into a completely unsupervised one.

We have assessed the superiority of the RPA method on several publicly available datasets and have used a heterogeneous panel of statistical tools to results. Also, we have included in our study other recent contributions from the literature, leading to a comprehensive comparison of the performance of state-of-the-art methods. We hope that the breath of the analysis performed here will be useful as a reference for future works related to Transfer Learning on the SPD manifold. In order to foster reproducible research, complete Python code for the results in this paper is available at https://github.com/plcrodrigues/RPA.

Future perspectives for this work shall include an online implementation of our method, where usual drifts in statistics from signals on the same session would be corrected via distribution matching. An important challenge for such procedure would be to detect when changes in the statistics occur as well as when the number of new trials is already large enough so that no information from data points drawn from previous statistical distributions are needed. Another interesting line of work would be to go further in the analysis of Section VI-D by extending the methods proposed in [12], [START_REF] Jayaram | Transfer Learning in Brain-Computer Interfaces[END_REF], and [START_REF] Waytowich | Spectral transfer learning using information geometry for a user-independent brain-computer interface[END_REF] with a statistical matching step based on the RPA. Finally, another interesting topic to investigate would be to consider a Transfer Learning approach for matching data from subjects having different numbers of electrodes.

VIII. ACKNOWLEDGEMENT

This work is partly supported by the ERC Grant CHESS 2012-ERC-AdG-320684. We thank Florent Bouchard for his help on the optimization procedure in Equation [START_REF] Jayaram | Moabb: Trustworthy algorithm benchmarking for bcis[END_REF].

Algorithm 1 :

 1 Transfer Learning via RPA Input: S, T and Tu as defined in (1) and (

u 6

 6 Form the training dataset for the MDM classifier with D train = S (rct) ∪ T (rot) and get the accuracy of classification on the data points from the test dataset Dtest = T (rot) u

Fig. 2 .

 2 Fig.2. Accuracy of the classification of unlabeled data points from the target dataset for different methods of transfer learning. The curve shows how the accuracy for each method evolves when the number of data points in T increases. The generation of the data points is explained in Section V-A.

Fig. 3 .

 3 Fig.3. Accuracies of the cross-subject classification for three different Transfer Learning procedures on the PhysionetMI database. The rows and columns of each subplot were reordered using the seriation procedure explained in the text. The colormap shows white for accuracies of 0.5 or less and black when it is 1.0. The compared methods are described in Section V-C and we consider the cases when there are one and ten labeled matrices in T .

Fig. 4 .

 4 Fig. 4. Results of the statistical tests on each pair of pipelines for all possible values of N on each dataset as indicated in Table II (for instance, on Cho2017 we have N 1 = 1, N 2 = 5, N 3 = 10, and N 4 = 25). The color/pattern of the squares indicate whether there's no statistical difference between two methods (white squares), if the Left method is superior to the Right one (L and R in the legend) (dark gray squares) or the contrary (squares with crossed patterns).All conclusions are with p < 0.05 corrected via Holm's adjustement[START_REF] Holm | A simple sequentially rejective multiple test procedure[END_REF]. For instance, we see that for dataset Cho2017, the method RPA is inferior to RCT when N = 1, but RPA becomes superior when N = 25. Furthermore, for this same dataset, there's no statistical difference in the comparison between PRL and RCT for any N .

Fig. 5 .

 5 Fig. 5. Scatter plots comparing the accuracies of the cross-subject classification on the MunichMI dataset for the RPA and Calibration methods. We consider two sizes for T . The percentage numbers indicate the proportion of dots above or below the diagonal line.

  Riemannian Procrustes Analysis : Transfer Learning for Brain-Computer Interfaces Pedro L. C. Rodrigues, Member, IEEE, Christian Jutten, Fellow, IEEE, and Marco Congedo, Member, IEEE

TABLE I MAIN

 I FEATURES DESCRIBING EACH DATASET USED IN THIS WORK.

	dataset	paradigm	subjects classes trials per class reference
	PhysionetMI	MI	109	2	22	[38]
	Cho2017	MI	50	2	100	[39]
	SSVEP	SSVEP	12	3	8	[37]
	P300	P300	24	2	72 and 360	[36]
	BNCI2014001	MI	9	4	72	[40]
	BNCI2014002	MI	15	2	80	[41]
	BNCI2015001	MI	13	2	100	[42]
	MunichMI	MI	11	2	150	[43]

TABLE II MEAN

 II VALUES OF THE CROSS-SUBJECT AUC (AREA UNDER THE ROC CURVE) FOR FIVE PIPELINES (ALL DESCRIBED IN SECTION V-C) ON EIGHT DIFFERENT DATASETS. PARAMETER N IS THE NUMBER OF TRAINING POINTS AVAILABLE ON THE target DATASET ON EACH SITUATION. THE BEST METHOD IN EACH INSTANCE IS WRITTEN IN BOLD.

	PhysionetMI	N 1 0.54 DCT 5 0.55 10 0.56 15 0.57 DCT 1 0.64 2 0.67 4 0.72 6 0.74	RCT 0.61 0.65 0.67 0.68 RCT 0.67 0.71 0.76 0.78	MEAN AUC PRL 0.61 0.65 0.67 0.68 MEAN AUC PRL 0.66 0.71 0.76 0.78	OPT 0.59 0.60 0.60 0.60 OPT 0.59 0.59 0.59 0.57	RPA 0.56 0.63 0.66 0.67 RPA 0.70 0.75 0.80 0.82	Cho2017 P300	N 1 0.54 DCT 5 0.55 10 0.55 25 0.57 N DCT 6 0.57 12 0.62 32 0.71 48 0.74	0.59 0.61 0.62 0.64 RCT 0.56 0.64 0.74 0.76	MEAN AUC PRL 0.58 0.61 0.62 0.64 PRL 0.56 0.64 0.74 0.76 MEAN AUC	OPT 0.57 0.57 0.57 0.58 OPT 0.58 0.61 0.67 0.69		RPA 0.54 0.59 0.62 0.66 RPA 0.55 0.63 0.73 0.75	BNCI2015001	N 1 0.58 DCT 6 0.59 18 0.61 36 0.64 N DCT 1 0.51 5 0.51 10 0.52 25 0.54		RCT 0.69 0.73 0.76 0.78 RCT 0.56 0.57 0.59 0.62	PRL 0.69 0.73 0.76 0.77 MEAN AUC MEAN AUC PRL 0.55 0.57 0.58 0.62	OPT 0.65 0.65 0.65 0.66 OPT 0.54 0.55 0.56 0.56	RPA 0.62 0.71 0.76 0.79 RPA 0.57 0.60 0.63 0.65		BNCI2014002	N 1 5 10 20 N 1 0.55 DCT 0.55 0.56 0.57 0.58 DCT 25 0.58 50 0.60 75 0.62	RCT 0.68 0.71 0.72 0.73 RCT 0.63 0.69 0.71 0.72	MEAN AUC PRL 0.67 0.70 0.71 0.72 PRL 0.63 0.69 0.71 0.72 MEAN AUC	OPT 0.64 0.64 0.65 0.65 OPT 0.61 0.62 0.62 0.62	RPA 0.63 0.70 0.72 0.73 RPA 0.55 0.68 0.72 0.73
							N 1
		RCT	DCT
		PRL	DCT
		OPT	DCT
		RPA	DCT
		RPA	RCT
		RPA	PRL
		PRL	RCT
		OPT	RCT
		OPT	PRL
		OPT	RPA
				PhysionetMI	Cho2017	P300-Erwan	SSVEP	BNCI2014001	BNCI2014002	BNCI2015001	MunichMI	PhysionetMI	Cho2017	P300-Erwan	SSVEP	BNCI2014001	BNCI2014002	BNCI2015001	MunichMI	PhysionetMI	Cho2017	P300-Erwan	SSVEP	BNCI2014001	BNCI2014002	BNCI2015001	MunichMI	PhysionetMI	Cho2017	P300-Erwan	SSVEP	BNCI2014001	BNCI2014002	BNCI2015001	MunichMI


Keywords: cocaine, crack, HTR2B gene, impulsivity, ADHD

Objectives: Cocaine dependence has a strong heritability component. This study was aimed to investigate the putative association between the serotonin 2B receptor gene (HTR2B), patients with crack use disorders and impulsivity. Methods: A French Caribbean male population of patients with crack use disorders (n=80) was compared to healthy male controls (n=60). Comorbid ADHD and impulsivity were assessed. Five Single Nucleotide Polymorphisms (SNPs) in the HTR2B gene on chromosome 2 were selected: rs643700, rs6736017, rs1549339, rs17586428 and rs3806545 These five SNPs were chosen to cover most of the linkage disequilibrium (LD) blocks in HTR2B. The French translation of the Baratt Impulsivity Scale BIS-11 was used to evaluate impulsivity. Comorbid ADHD was diagnosed using the Wender Utah Rating Scale-25 item for the Attention Deficit-Hyperactivity Disorder. Results: We have found a positive association between the rs6736017 polymorphism and crack use disorders in a French Afro-Caribbean male population. Conclusion: The risk effect of HTR2B rs6736017 appeared to be specific to individuals who are crack users rather than being driven by impulsivity alone or ADHD alone.

A positive association between a polymorphism in the HTR2B gene and cocaine-crack in a French Caribbean population

Jerome Lacoste1 , Sandrine Lamy 2 , Nicolas Ramoz 3 , Nicolas Ballon 4 , Louis Jehel 5 , Luc Maroteaux 6 , Florence Thibaut 7 .

Introduction

Cocaine addiction remains a serious public health concern worldwide. The annual prevalence rate of cocaine use (defined as the number of people aged 15-64 years who have consumed the drug at least once during the last 12 months) is 0.35 worldwide with important regional differences between North America (1.77), Oceania (1.54), Western and Central Europe (1.08) and Asia (0.04) [START_REF] Peacock | Global statistics on alcohol, tobacco and illicit drug use: 2017 status report[END_REF]. When compared to controls with no history of substance abuse, cocaine and/or crack users exhibited higher levels of impulsiveness (for review [START_REF] Czermainski | Assessment of inhibitory control in crack and/or cocaine users: a systematic review[END_REF]. In addition, attention deficit-hyperactivity disorders (ADHD) and substance use disorders often coexist in adults, the risk of developing cocaine abuse is two times higher in individuals with (vs. without) ADHD; the prevalence of ADHD is also high among patients with cocaine dependence (20-25%) (for review [START_REF] Miguel | Cognition and impulsivity in adults with attention deficit hyperactivity disorder with and without cocaine and/or crack dependence[END_REF]. Using the BIS-11 scale, ADHD patients who asked for treatment of cocaine use disorders, compared with those without ADHD, had higher motor impulsivity (Pérez de los [START_REF] De | Features and prevalence of patients with probable adult attention deficit hyperactivity disorder who request treatment for cocaine use disorders[END_REF].

Cocaine dependence has a strong heritability component. Siblings of cocaine-dependent subjects had an estimated relative risk of 1.8 of developing cocaine dependence compared with those without cocaine-dependent siblings [START_REF] Bierut | Familial transmission of substance dependence: alcohol, marijuana, cocaine, and habitual smoking: a report from the Collaborative Study on the Genetics of Alcoholism[END_REF]. A higher rate of concordance was also observed among monozygotic twins relative to dizygotic twins [START_REF] Kendler | Specificity of genetic and environmental risk factors for use and abuse/dependence of cannabis, cocaine, hallucinogens, sedatives, stimulants, and opiates in male twins[END_REF]. A relatively common stop codon (Q20*) was identified in the serotonin 2B receptor gene (HTR2B) in a Finnish founder population and it was associated with impulsivity [START_REF] Bevilacqua | population-specific HTR2B stop codon predisposes to severe impulsivity[END_REF]. Moreover, [START_REF] Tikkanen | Impulsive alcohol-related risk-behavior and emotional dysregulation among individuals with a serotonin 2B receptor stop codon[END_REF] have reported that the presence of the HTR2B Q20* predicted impulsive and aggressive behaviors particularly under the influence of alcohol.

Although crack is derived from cocaine, it has distinct patterns of consumption and potential for addiction [START_REF] Gossop | Cocaine: patterns os use, route of administration and severity of dependence[END_REF][START_REF] Gossop | Concurrent use of alcohol and cocaine : differences in patterns of use and problems among users of crack cocaine and cocaine powder[END_REF]). Yet, there is a lack of studies investigating the genetic factors involved in patients with crack use disorders, whose use is gradually increasing. This study was aimed to investigate the putative association between the serotonin 2B receptor gene (HTR2B), crack-cocaine use disorders and impulsivity in a French Caribbean crack-cocaine male population of patients with crack use disorders as compared to healthy male controls.

Subjects and Methods

Sixty healthy male controls and eighty French Afro-Caribbean male patients with crack use disorders were included in this study. Participants were diagnosed patients with crack use disorders using DSM-IV-R criteria. Healthy controls were consecutively recruited at the blood transfusion center of Martinique, without lifetime personal and family history of psychiatric disorder, including substance abuse. The clinical evaluations were performed by a psychiatrist. All subjects (patients and controls) were male, African Caribbean, and had at least three African-Caribbean grandparents. Healthy controls and patients with crack use disorders were matched for age.. All subjects agreed to participate in the study and signed an informed consent. Medical and drug use histories, socio-demographic characteristics were reported using self-report and medical files. Personal and family psychiatric history of patients and controls were evaluated by a psychiatrist or a psychologist, blind to the genotyping results and trained to the Diagnostic Interview of Genetic Studies (DIGS). [START_REF] Nurnberger | Diagnostic interview for genetic studies. Rationale, unique features, and training. NIMH Genetics Initiative[END_REF]. The French translation of the Baratt Impulsivity Scale BIS-11 [START_REF] Patton | Factor structure of the Barratt Impulsiveness Scale[END_REF]Bayle et al., 2000) was used to evaluate impulsivity. A score above 71 indicates a high level of impulsivity [START_REF] Stanford | Fifty years of the Barratt Impulsiveness Scale: An update and review[END_REF]. Comorbid ADHD was diagnosed using the Wender Utah Rating Scale-25 item for the Attention Deficit-Hyperactivity Disorder with a cut off value at 46 [START_REF] Ward | The Wender Utah Rating Scale: An aid in the retrospective diagnosis of childhood Attention Deficit Hyperactivity Disorder[END_REF]French translation Baylé et al., 2003). Sensation seeking (using the Zuckerman Sensation-Seeking Scale form V-40 item [START_REF] Zuckerman | Sensation seeking in England and America : cross-cultural, age and sex comparisons[END_REF]French translation Carton et al., 1990) and comorbid psychiatric disorders (using the MINI [START_REF] Lecrubier | The Mini International Neuropsychiatric Interview (MINI). A short diagnostic structured interview: reliability and validity according to the CIDI[END_REF] were also assessed (Table 1). DNAs were collected from blood or saliva samples. A 4 ml venous blood or 2 ml saliva sample (Oragene DNA OG-500 self-collection kit) was collected from each subject. Genomic DNA was extracted according to standard methods. All SNP in samples were genotyped using the Taqman assays (Thermo Fisher). DNAs were genotyped for HTR2B polymorphisms (Tables 2 and3). The HTR2B gene contains 4 exons and 3 introns. As described in Figures 1 and2 and Tables 2 and3, five SNPs in the HTR2B gene on chromosome 2 were selected based on the 1000 Genomes project (https://www.ncbi.nlm.nih.gov/variation/tools/1000genomes/) including rs643700, rs6736017, rs1549339, rs17586428 and rs3806545. All five SNPs met the criteria of a call rate of >95%; minor allele frequency of >0.05; and Hardy-Weinberg equilibrium of p > 0.05. The allele frequencies in our sample were very similar to those of the same ethnicity in the HapMap dataset. These five SNPs were chosen to cover most of the linkage disequilibrium (LD) blocks in HTR2B (Figures 2 and3). The distribution of the SNPs encompassing HTR2B gene were in Hardy-Weinberg equilibrium.

Statistical analyses

Data analyses were conducted using chi square, ANOVA, Fisher's exact tests and a general linear model. Fisher's exact test interactions between ADHD or impulsivity (qualitatively) and the association between crack status and polymorphisms were tested. A general linear model was used to test interactions between ADHD or impulsivity (quantitatively) and the association between crack status and polymorphisms. Significance was set at p < 0.05.

Results

No differences were observed between controls and patients with crack use disorders for age, sensation seeking or BIS attentional scores. Using the WURS-25, and a cut off score of 46, the prevalence of ADHD among patients with crack use disorders (71.4%, n=40) was significantly higher compared with controls (27.3%, n=12) (p<0.001). Using the BIS-11, compared to controls, patients with crack use disorders had a higher level of total, motor and nonplanning impulsivity (p<0.001) (Table 1). Similarly, using a cut off score of 71 (BIS total > 71), patients with crack use disorders had a significantly higher prevalence of total impulsivity (21.7% n=13) as compared to controls (3.4%, n=2) (p<0.004). Therefore statistical analyses were adjusted for these variables. A statistically significant association (p=0.029) was observed between the rs6736017 polymorphism allele T and patients with crack use disorders (Table 2). The level of significance was even higher when 12 patients with crack use disorders with a comorbidity of schizophrenia or bipolar disorder were excluded from the sample, p=0.015, Odds Ratio 2.99, confident interval of 95% [1.19-7.49] (Table 3). A corrected p value of 0.075 was carried out for multiple tests, showing a trend for association between the rs6736017*T and crack use disorders. A tendency was observed for rs1549339 polymorphism with patients with crack use disorders (Table 2) but it disappeared when psychiatric comorbidity (schizophrenia or bipolar disorder) was taken into account (Table 3). There was neither interaction between ADHD, crack status and any polymorphism studied (Supplementary Table 1), nor between a high score of impulsiveness, crack status and any polymorphism studied (Supplementary Table 2). Using a general linear model to study the interactions between the scores obtained using the BIS-11 and the WURS-25, clinical status × HTR2B genotype interactions were not significant suggesting that the relationship between impulsivity or ADHD and HTR2B genotype in patients with crack use disorders was not different from that observed in controls (Supplementary Tables 3 &4). There was no clinical status × HTR2B genotype interaction for novelty seeking.

Discussion

We have observed a positive association between the rs6736017 polymorphism (a missense variant located in exon 4 of the HTR2B gene) and crack use disorders in a French Afro-Caribbean male population. The risk effect of HTR2B rs6736017 appeared to be specific to patients with crack use disorders rather than being driven by impulsivity or ADHD alone; however these results need replication in larger sample of patients as well as in populations of different genetic backgrounds. The level of impulsiveness and the prevalence of ADHD were actually significantly higher in our population of patients with crack use disorders as compared to controls, but a lack of power may have contributed to the lack of interaction observed. One unsolved question is the putative functional consequence of the rs6736017 polymorphism introducing a Met421/Val substitution at the C-terminal cytoplasmic tail of the 5-HT2B receptor (Supplementary Figure 3), which was not found in any sequenced HTR2B in primates. [START_REF] Bevilacqua | population-specific HTR2B stop codon predisposes to severe impulsivity[END_REF] observed that a stop codon in exon 2 (HTR2B Q20*) led to an interruption in the expression of the 5-HT 2B receptor and found that this mutation was associated to impulsive behavior. Among these HTR2B Q20* carriers, temperamental traits resembled a passive-dependent personality profile, and the presence of the HTR2B Q20* predicted impulsive and aggressive behaviors particularly under the influence of alcohol [START_REF] Tikkanen | Impulsive alcohol-related risk-behavior and emotional dysregulation among individuals with a serotonin 2B receptor stop codon[END_REF]. These impulsive/aggressive traits have also been found in mice lacking 5-HT 2B receptor gene (Htr2b -/-) [START_REF] Bevilacqua | population-specific HTR2B stop codon predisposes to severe impulsivity[END_REF]. The HTR2B gene may function by modulating levels of both serotonin and dopamine in brain regions involved in impulsivity [START_REF] Doly | Serotonin 5-HT2B receptors are required for 3,4 methylenedioxymethamphetamine-induced hyperlocomotion and 5HT release in vivo and in vitro[END_REF]. In agreement witht the fact that, to date, the HTR2B Q20* was exclusively reported in a Finnish founder population (0.012% of the population) [START_REF] Bevilacqua | population-specific HTR2B stop codon predisposes to severe impulsivity[END_REF], we have not found this stop codon in our population. Furthermore, a recent genome-wide association study (GWAS) of physical aggression occurring under the influence of cannabis in African Americans implicated the HTR2B-rs17440378 TT genotype, which was not previously considered a key cannabis target [START_REF] Montalvo-Ortiz | Translational studies support a role for serotonin 2B receptors in aggression-related cannabis response[END_REF]. Interestingly, treating the 5-HT 2B receptor knockout mice with THC resulted in increased aggressive behavior in the resident-intruder paradigm, while the opposite effect was observed in wild type mice. Thus, HTR2B variation in expression moderates the effects of cannabis on aggression and social interactions [START_REF] Montalvo-Ortiz | Translational studies support a role for serotonin 2B receptors in aggression-related cannabis response[END_REF]. Interestingly, null expression or blockade of 5-HT 2B receptors diminishes the reinforcing effects of psychoactive drugs by modulating serotonin and dopamine signaling. The direct effects of HTR2B-rs17440378 (an intronic polymorphism) on protein function, including regulation is unknown. In general, crack users have a higher level of impulsivity as compared to cannabis users, which means that it may be more difficult to disentangle the respective roles of impulsivity and crack use disorders effects in the association with SNPs, especially in a small sample. Finally, a major haplotype A-A-G-A (with a frequency of 48%) for rs6437000-rs10194776-rs16827801-rs1549339 (four SNPs located in introns 2 and 3) showed a significant positive association with behavioral-approach-system (BAS) fun seeking in a Han Chinese population [START_REF] Zhu | Association between the HTR2B gene and the personality trait of fun seeking[END_REF]. BAS fun seeking was also correlated with impulsive behaviors such as drug abuse. Among these four significant SNPs, rs6437000 and rs1549339 were reported to have high linkage disequilibrium with the rare HTR2B stop codon Q20*. According to [START_REF] Bevilacqua | population-specific HTR2B stop codon predisposes to severe impulsivity[END_REF], HTR2B Q20* was found on a single haplotype background A-G-A-(G/A)-A-A for rs6437000-rs4973377-rs1549339-(stop codon)-rs17586428-rs3806545. In contrast, [START_REF] Tsuchimine | No association between a polymorphism in the serotonin receptor 2B (HTR2B) gene and personality traits in healthy Japanese subjects[END_REF] have reported that the HTR2B rs10194776 polymorphism was not associated with personality traits, including novelty seeking and impulsivity in healthy Japanese subjects. The 5-HT 2B receptors are expressed in the adult human brain, especially in the frontal cortex, although at a low level. The recent work by [START_REF] Doly | Serotonin 2B receptors in mesoaccumbens dopamine pathway regulate cocaine responses[END_REF] has shown a stronger locomotor response to cocaine, with a doubling the distance traveled by 5-HT 2B receptor mutant mice. The 5-HT 2B receptors involved have been located on dopaminergic neurons, who innervate the nucleus accumbens (NA), a brain area involved in the reward circuit. Their absence decreases the activity of these latter neurons. Paradoxically, the lack of 5-HT 2B receptors induces a decrease of dopamine secretion in the NA [START_REF] Doly | Serotonin 2B receptors in mesoaccumbens dopamine pathway regulate cocaine responses[END_REF]. These results indicate that, by binding on 5-HT 2B receptors, serotonin can control the amount of dopamine released by dopaminergic neurons that express the receptor. When these receptors are lacking, this secretion is abnormally low in the NA. Several imaging studies in humans have shown that individuals addicted to cocaine present paradoxically a weaker activation of the NA and a decreased secretion of dopamine in this region [START_REF] Volkow | Stimulant-induced dopamine increases are markedly blunted in active cocaine abusers[END_REF]. Thus, the absence of the 5-HT 2B receptor induces a neuronal situation similar to that observed in drug addicts. This similarity between these mutant mice and cocaine addicts suggests that individuals bearing a genetic mutation that invalidates the 5-HT 2B receptor or associated to a dysfunction of this receptor could have a greater propensity to become addicted to this drug. Similarly, the chronic 5-HT 2B -receptor inhibition makes mice behaving as animals already exposed to cocaine with higher cocaine-induced locomotion. Associated changes in dopaminergic neuron reactivity that includes a local hypodopaminergy, paradoxically and ultimately results in increased cocaine psychostimulant responses without significant change in motivation to the drug whatever the level of impulsivity or hyperactivity observed.

Strengths and limitations of our study

Our results need replication in larger populations, especially, to conclude that there is an association between genetic variants of HTR2B gene and crack use disorders, independently to ADHD and impulsivity. We calculated a power analysis using sample size calculator and a risk alpha of 0.05. Considering our sample (80 cases and 60 controls), the percentages of rs6736017*T (0.94 and 0.86 respectively), the estimated power was of 36%. To detect an effect of this magnitude with a power of 80%, one would need to recruit 250 cases and 188 controls. Furthermore, for the interaction studies (with ADHD and impulsivity), the estimated power was of 22% and 32% respectively. Nevertherless, we observed a positive association between crack use disorders and a polymorphism located on the HTR2B gene in our sample of 80 male French Afro-Caribbean patients. It should be noted, however, that the SNPs used in this study have different minor allele frequency (MAF) in different ethnic populations based on HapMap Data (http://www.hapmap.org) and that replication studies on populations of different genetic backgrounds are needed. The selection of SNPs was based on a high frequency of the heterozygosity and of the minor allele in Caucasian populations since no African Caribbean population was available at the time of our study Fundings LM's work has been supported by funds from the Centre National de la Recherche Scientifique, the Institut National de la Santé et de la Recherche Médicale, the Sorbonne Universités, and by grants from the Fondation pour la Recherche sur le Cerveau, the Fondation pour la Recherche Médicale "Equipe FRM DEQ2014039529", the French Ministry

Table 1 .

 1 Characteristics of patients and controls.

	Clinical characteristics	Controls Mean+/S.D. (n=60)	Patients Mean+/S.D. (n=80)	d.f.	F	P value
	Age	45.6+9.2	43.2+7.5	1/117	2.31	0.13
	WURS 25	40.0+17.7	60.6+26.6	1/98	19.49	<0.001*
	WURS total	52.2+23	83.0+32.9	1/98	27.87	<0.001*
	BIS 11 motor	6.5+4.5	22.2+5.8	1/117 266.34	<0.001*
	BIS 11 attentional	17.2+7.3	17.3+3.5	1/117	0.004	0.95
	BIS 11 nonplanning	13.2+7.6	24.4+5.5	1/117	84.48	<0.001*
	BIS 11 total	37+16.9	63.9+11.9	1/117 101.09	<0.001*
	SSS total	15.5+5.5	16.6+5.6	1/119	1.09	0.30
	An ANOVA was performed to compare patients and controls clinical characteristics.
	d.f.: degree of freedom; N: number: S.D.: Standard Deviation
	* Significant

Table 2 .

 2 Allele frequencies and associations of five candidate SNPs in the HTR2B gene and clinical status (all patients with crack use disorders (n=80) or controls (n=60)).

	Polymorphism Allele	Case/Control	Case/Control,	Chi square P value
			counts	frequencies
	rs6736017	T	148:10/ 103:17	0.94/ 0.86	4.78	0.029*
	rs6437000	C	48:11/31:89	0.30/0.26	0.59	0.44
	rs1549339	G	64:90/36:80	0.42/0.31	3.14	0.07
	rs17586428	A	111:47/77:41	0.70/0.65	0.78	0.38
	rs3806545	A	107:49/75:43	0.69/0.64	0.76	0.38
	*Significant
	rs6736017: genetic association OR=2.44 CI95%[1.06-5.55]
	rs1549339: trends of association

Table 3 .

 3 Allele frequencies and associations of five candidate SNPs in the HTR2B gene and clinical status (patients with crack use disorders (n=68) or controls (n=60)).

	Polymorphism Allele	Case/Control	Case/Controls,	Chi	P value
			counts	frequencies	square
	rs6736017	T	127:7/ 103:17	0.95/ 0.86	5.92	0.015*
	rs6437000	C	37:99/31:89	0.27/0.26	0.06	0.80
	rs1549339	G	51:81/36:80	0.39/0.31	1.57	0.21
	rs17586428	A	91:43/77:41	0.68/0.65	0.2	0.65
	rs3806545	A	87:45/75:43	0.66/0.64	0.15	0.70
	*Significant
	rs6736017: genetic association OR=2.99 CI95%[1.20-7.50]

Table 2 .

 2 No association was observed between SNPs and impulsivity in patients with crack use disorders compared to controls Crack use disorders with a BIS > 50 (medium + high level of impulsivity) N: Number ; SNP: single nucleotide polymorphism ; %: percentage * All cases of patients with crack use disorders (Total sample: N=80) ** Patients with crack use disorders without comorbid schizophrenia or bipolar disorder (Total sample: N=68)

Jerome Lacoste MD PhD CHU Martinique, Department of Psychiatry and Addictology, Fort-de-France, Martinique;

of Research (Agence Nationale pour la Recherche ANR-17-CE16-0008 and the Investissements d'Avenir programme ANR-11-IDEX-0004-02). LM's team is part of the École des Neurosciences de Paris Ile-de-France network and of the Bio-Psy Labex and as such this work was supported by French state funds managed by the ANR within the Investissements d'Avenir programme under reference ANR-11-IDEX-0004-02. JL's work (CRACK-ANT study) was sponsored by the University Hospital of Martinique for regulatory and ethic submission and was partially supported by a grant from the Clinical Research and Innovation Unit (DRCI) of the University Hospital of Martinique (2010).

Acknowledgements

The authors thank P Mazeau for editing this manuscript.

Statement of Interest

Florence Thibaut is Editor-in-Chief of Dialogues in Clinical Neuroscience (the journal receives a grant from Servier). Other authors have no conflict of interest to declare.

Legends