Di Ge
email: ge.di@glaizer.com

Jérôme Idier (Member, IEEE Eric Le Carpentier

J Idier

D Farina
email: dario.farina@bccn.uni-goettingen.de

Spike Sorting by Stochastic Simulation

Keywords: Bayesian model, MMSE estimation, Markov

The decomposition of multiunit signals consists of the restoration of spike trains and action potentials in neural or muscular recordings. Because of the complexity of automatic decomposition, semiautomatic procedures are sometimes chosen. The main difficulty in automatic decomposition is the resolution of temporally overlapped potentials. In a previous study [1], we proposed a Bayesian model coupled with a maximum a posteriori (MAP) estimator for fully automatic decomposition of multiunit recordings and we showed applications to intramuscular EMG signals. In this study, we propose a more complex signal model that includes the variability in amplitude of each unit potential. Moreover, we propose the Markov Chain Monte Carlo (MCMC) simulation and a Bayesian minimum mean square error (MMSE) estimator

by averaging on samples that converge in distribution to the joint posterior law. We prove the convergence property of this approach mathematically and we test the method representatively on intramuscular multiunit recordings. The results showed that its average accuracy in spike identification is greater than 90% for intramuscular signals with up to 8 concurrently active units. In addition to intramuscular signals, the method can be applied for spike sorting of other types of multiunit recordings.

The decomposition problem of multiunit recordings, often referred to as spike sorting, is usually solved with template matching approaches [START_REF] Lewicki | A review of methods for spike sorting: the detection and classification of neural action potentials[END_REF]. In some of these approaches, the action potentials that are overlapped in time are not classified and are considered outliers [START_REF] Quian Quiroga | Unsupervised spike detection and sorting with wavelets and superparamagnetic clustering[END_REF]. When the number of sources is small, this limitation may provide an acceptable amount of information on the sources under study. However, when there are many sources active concurrently, overlapped action potentials constitute the majority of the potentials in the recorded signals. Therefore, in some template matching methods, the superpositions are resolved by iterative subtraction of all possible template combinations from unidentified waveforms [START_REF] Lewicki | Bayesian modeling and classification of neural signals[END_REF], [START_REF] Bankman | Optimal detection, classification, and superposition resolution in neural waveform recordings[END_REF]. As an alternative, neural network classifiers are applied to resolve the superposition problem by introducing overlapped spikes into the training data [START_REF] Chandra | Detection, classification, and superposition resolution of action potentials in multiunit singlechannel recordings by an on-line real-time neural network[END_REF].

The main challenge underlying the resolution of superimposed spikes is that the global optimization problem is non-deterministic polynomial-type (NP) hard, i.e., it cannot be solved by polynomial complexity algorithms [START_REF] Ge | Unsupervised Bayesian Decomposition of Multi-Unit EMG Recordings using Tabu Search[END_REF]. Therefore, the existing methods either perform on the restrained search spaces [START_REF] Lewicki | Bayesian modeling and classification of neural signals[END_REF], [START_REF] Atiya | Recognition of multiunit neural signals[END_REF], which reduces the complexity, or are based on recursive algorithms [START_REF] Mcgill | Optimal resolution of superimposed action potentials[END_REF], [START_REF] Mcgill | EMGLAB: An interactive EMG decomposition program[END_REF] with specific trial strategies and residual threshold estimations. For example, Atiya [START_REF] Atiya | Recognition of multiunit neural signals[END_REF] proposed a robust approach for decomposing overlaps of action potentials in neural recordings by comparing all possible combinations of up to two action potentials (restrained search space). Within the family of spike sorting algorithms, the dynamic programming method [START_REF] Lewicki | Bayesian modeling and classification of neural signals[END_REF], which uses the fast exploration technique of a k-d tree, is also limited by the memory space necessary to generate such data structure, resulting in practice in an equally restrained search space of up to two overlapping action potentials. This constraint is not justified in several applications. Similar issues arise when determining a representative training set of overlapped spikes and then using neural networks to identify overlapped sources [START_REF] Chandra | Detection, classification, and superposition resolution of action potentials in multiunit singlechannel recordings by an on-line real-time neural network[END_REF].

Other approaches are not limited in the number of overlapping sources but require an interaction with an operator. For example, a recent algorithm for decomposing intramuscular electromyographic (EMG) signals into the constituent motor unit spike trains recursively matches the templates and reevaluates the residual errors, until it is able to prove that it has found the global optimum [START_REF] Mcgill | Optimal resolution of superimposed action potentials[END_REF]. To avoid too long recursions in cases of difficult spike superpositions, the algorithm stops and requires an expert intervention after a number of trials without reaching the lower bound value.

Other methods make use of independent component analysis (ICA), of clustering and ICA [START_REF] Takahashi | Automatic sorting for multineuronal activity recorded with tetrodes in the presence of overlapping spikes[END_REF], [START_REF] Takahashi | Real-time and automatic sorting of multi-neuronal activity for sub-millisecond interactions in vivo[END_REF], [START_REF] Takahashi | Classification of neuronal activities from tetrode recordings using independent component analysis[END_REF] to avoid the limitation on the number of units concurrently active and the need for the intervention of an operator. However, these methods are applied to multi-channel recordings. Similarly, blind source separation (BSS) approaches such as the Convolution Kernel Compensation (CKC) used for surface EMG decomposition [START_REF] Holobar | Correlation-based decomposition of surface electromyograms at low contraction forces[END_REF], [START_REF] Holobar | Estimating motor unit discharge patterns from highdensity surface electromyogram[END_REF], are usually applied to multi-channel recordings for which at least as many channels as sources are needed.

A Bayesian model coupled with a maximum a posteriori (MAP) estimator was recently proposed for fully automatic decomposition of multiunit recordings in singlechannel signals [START_REF] Ge | Unsupervised Bayesian Decomposition of Multi-Unit EMG Recordings using Tabu Search[END_REF]. The method was tested on intramuscular EMG recordings to identify individual motor unit spike trains and proved to have similar performance as obtained with a semi-automatic decomposition by expert operators [START_REF] Ge | Unsupervised Bayesian Decomposition of Multi-Unit EMG Recordings using Tabu Search[END_REF]. Furthermore, the number of concurrent action potentials was not limited although the search space augmented exponentially.

One assumption of the Bayesian model in [START_REF] Ge | Unsupervised Bayesian Decomposition of Multi-Unit EMG Recordings using Tabu Search[END_REF] is that action potentials discharged by a unit have the same amplitude and shape for the entire recording. However, although the shape of action potentials may not change substantially, the amplitude may be variable in some conditions. For example, the amplitude of single motor unit action potentials in intramuscular EMG may be influenced by the rate of discharge because of the velocityrecovery function of muscle fibers [START_REF] Farina | Effect of muscle-fiber velocity recovery function on motor unit action potential properties in voluntary contractions[END_REF]. Moreover, small displacements of the recording electrodes with respect to the sources may influence action potential amplitude.

In this study, we address the spike sorting problem in single-channel recordings by proposing a more complex signal model than in [START_REF] Ge | Unsupervised Bayesian Decomposition of Multi-Unit EMG Recordings using Tabu Search[END_REF]. This model includes the variability in magnitudes of each unit potential to provide a better fit with the experimental signals in a leastsquared-residual sense. Since action potentials may also change in shape, which is not modeled, this new model is not expected to provide a perfect fit. However, it is hypothesized that it would provide better overall performance compared with the invariant amplitude and shape model. Moreover, instead of the MAP estimator [START_REF] Ge | Unsupervised Bayesian Decomposition of Multi-Unit EMG Recordings using Tabu Search[END_REF], we propose a new approach for the solution of the spike sorting problem. The approach is based on the Markov Chain Monte Carlo (MCMC) simulation tool to build a minimum mean square error (MMSE) estimator for continuous parameters and a marginal MAP estimator for discrete parameters (discharge instants) using samples that converge in distribution to the joint posterior law.

The study also provides the mathematical proof of the convergence property of this approach. As in [START_REF] Ge | Unsupervised Bayesian Decomposition of Multi-Unit EMG Recordings using Tabu Search[END_REF], the proposed method is representatively tested on intramuscular multiunit recordings, although it can also be directly applied to other multiunit signals. The test on intramuscular EMG signals will allow a direct performance comparison with the MAP estimator previously proposed [START_REF] Ge | Unsupervised Bayesian Decomposition of Multi-Unit EMG Recordings using Tabu Search[END_REF]. The main part of the validation is performed on experimental data, however we also present results on EMG signals simulated with very irregular spike trains as a proof of its robustness.

II. MODEL WITH VARIABLE IMPULSE MAGNITUDES

The generation of a multiunit recording is schematically represented in Fig. 1. Each impulse train is convolved by a characteristic action potential whose shape is assumed to vary minimally during the recording. The direct problem model can thus be expressed mathematically as a convolution product:

z = I i=1 h i * s i + ǫ (1)
where z is the recorded signal of length N , s i , i = 1, . . . , I and h i , i = 1, . . . , I are the impulse trains (the discharge patterns) and their linear filters (action potentials) respectively. Note that in the present study, the impulses in s i can have different amplitudes reflecting the amplitude variation of the different action potentials in the train. (•) i denotes the set {•, i = 1 . . . , I}, e.g., (s i ) i = {s i , i = 1 . . . , I} gathers the discharge patterns of all impulse trains. Assuming that ǫ is an independent, identically distributed (i.i.d.) Gaussian noise with unknown variance σ 2 ǫ , the likelihood of the data, given the source parameters (s i , h i ) i , σ 2 ǫ can be written:

P (z | (s i , h i ) i , σ 2 ǫ ) = 1 (2πσ 2 ǫ ) N 2 e - | z-i s i * h i | 2 2σ 2 ǫ
(2)

A. Prior laws on parameters

The modeling of the impulse trains (s i ) i is based on the following assumptions and notations: A1) (s i ) i are supposed mutually independent discrete time impulse trains; A2) the coordinates of non-zero elements (discrete discharge instants) in each s i are contained in a vector noted as x i ; A3) each s i conditional to x i and σ 2 si is supposed to be a truncated Gaussian:

s i | x i , σ 2 si ∼ N (1 xi , σ 2 si diag(1 xi )) • 1 si≥0 (3)
where 1 xi is binary vector with ones at positions x i , and diag(1 xi ) is the diagonal matrix whose diagonal is 1 xi . The sparsity in each spike train is modeled by the vector x i , whereas the variability in the amplitude of spikes generated by the same source is taken into account via the conditional Gaussian law (Eq.( 3)). We note that negative magnitudes are not relevant in spike sorting, thus the law of the magnitudes should be a Gaussian truncated to the positive hyperoctant, i.e., 1 si≥0 imposes a non-negative constraint on each element of s i .

In the literature, the Bernoulli-Gaussian (BG) model has been adopted for the restoration of a sparse and magnitude variable train [START_REF] Kormylo | Maximum-likelihood seismic deconvolution[END_REF], [START_REF] Cheng | Simultaneous wavelet estimation and deconvolution of reflection seismic signals[END_REF], [START_REF] Champagnat | Unsupervised deconvolution of sparse spike trains using stochastic approximation[END_REF], [START_REF] Bourguignon | Bernoulli-Gaussian spectral analysis of unevenly spaced astrophysical data[END_REF]. The BG model first defines the discharge instants x i using the Bernoulli law:

P (x i ) = λ ni (1 -λ) N -ni
where n i = dim(x i ) denotes the unknown number of discharges, and λ ∈ (0, 1) the Bernoulli parameter; and then defines s i conditionally to x i as a Gaussian variable, in the same manner as in Eq. (3). Let us remark that a Gaussian random variable with zero-variance corresponds to a variable of constant value (the Gaussian mean m), for which the distribution can also be defined as a Dirac.

On the other hand, the Bernoulli model needs also adaptations when used for neural spike trains, in order to account for physiological constraints. In this study, the model of the impulse instants (x i ) i is based on the following assumptions: A4) the inter-spike interval (ISI) T ij = x i,j+1x i,j , or the temporal distance between two consecutive impulses of the same source, is assumed to be larger than a threshold value T R . This value represents the refractory period needed to discharge the spikes and thus depends on the specific application; e.g., for intramuscular EMG decomposition, it is the absolute refractory period of muscle fibers;

A5) among all the admissible solutions satisfying the refractory period condition, spike trains that are more regularly spaced are favored. Following [START_REF] Ge | Unsupervised Bayesian Decomposition of Multi-Unit EMG Recordings using Tabu Search[END_REF], this is achieved by a Gaussian-shaped distribution on the variables T ij -T R . The resulting law of x i (in a discrete, regularly sampled time framework) for each unit can be written:

P x i | m i , σ 2 i ∝ 1 C (x i ) f -(ni-1) mi,σ 2 i exp    - 1 2σ 2 i ni-1 j=1 (x i,j+1 -x i,j -m i -T R ) 2    (4)
where

f mi,σ 2 i = +∞ k=-mi exp - k 2 2σ 2 i ,
and C = {x i , x i,j+1x i,j > T R for all j} is the set of admissible x i satisfying the refractory period condition, 1 C is the indicator function of C:

1 C (x i ) = ni-1 j=1 1 {xi,j+1-xi,j>TR} .
Contrarily to the BG case, a spike train s i is no more an i.i.d. sequence, because the discharge instants x i are not independent. In the BG deconvolution framework, a comparable model called modified Bernoulli process has been proposed in [START_REF] Kail | A blind Monte Carlo detection-estimation method for optical coherence tomography[END_REF] to impose the minimum distance constraint on the impulse train. The integration constant f mi,σ 2 i depends on both m i and σ 2 i . In practice, the following approximation can be adopted:

f mi,σ 2 i ≈ +∞ -mi exp - k 2 2σ 2 i dk ≈ +∞ -∞ exp - k 2 2σ 2 i dk = 2πσ 2 i , (5)
under the assumptions that • σ i is large enough to ignore the discretization error; e.g., for σ i = 10ms and a sampling frequency of 10kHz, the 3σ Gaussian lobe is discretized by 600 samples; • σ i /m i is small enough to ignore the integral truncation at -m i , typically if σ i /m i < 1/3, then the truncation error is controlled at 0.135%. Finally, in order to define a proper probability, Eq. ( 4) should be modified to incorporate a prior law on the first discharge instant x i1 . Here we suppose it is uniform for the sake of simplicity.

Using the independence assumption A1, we obtain:

P (x i ) i | (m i , σ 2 i ) i = I i=1 P x i | m i , σ 2 i . (6)

B. Joint posterior law

In the Bayesian framework, the posterior distribution

P (Θ | z) for Θ = (x i , s i , m i , σ 2 i , h i , σ 2 si ) i , σ 2 ǫ
can be expressed:

P (Θ | z) ∝ P (z | (s i , h i ) i , σ 2 ǫ )P (σ 2 ǫ ) i P (s i | x i , σ 2 si )P (σ 2 si )P (x i | m i , σ 2 i )P (m i )P (σ 2 i )P (h i ) (7)
Fig. 2 illustrates the hierarchical Bayesian model with all intervening parameters, among which unknown parameters (in Θ) are circled. The decomposition task consists of deriving an estimator Θ from Eq. ( 7). We note that Θ contains both continuous and discrete parameters and that the combinatorial nature of (x i ) i , precludes the exhaustive exploration method even in one segment [START_REF] Atiya | Recognition of multiunit neural signals[END_REF].

Eqs. ( 2), ( 3) and ( 4) provide the likelihood term and the prior laws on (s i , x i ) i that enter the joint posterior [START_REF] Atiya | Recognition of multiunit neural signals[END_REF]. For the remaining terms, we choose conjugate priors with non-informative hyper-parameters [START_REF] Ge | Unsupervised Bayesian Decomposition of Multi-Unit EMG Recordings using Tabu Search[END_REF] for a discussion on this choice):

(α i , β i ) i , (α s , β s , µ 0 , σ 2 0 , σ 2 h ) (see
m i ∼ N (µ 0 , σ 2 0 ), h i ∼ N (h (0) i , σ 2 h I), σ 2 i ∼ IG(α i , β i ), σ 2 ǫ ∼ IG(α ǫ , β ǫ ), σ 2 si ∼ IG(α s , β s ),
where IG stands for the inverse Gamma distribution and I is an identity matrix of appropriate size. (h

i ) i denote approximate spike shapes that are determined by the preprocessing step described in the following section.

h i s i posterior likelihood σ 2 ǫ m i σ 2 i m 0 , σ 2 0 α i , β i x i α s , β s z h (0) i , σ 2 h σ 2 si a s , b s

III. PREPROCESSING

The resolution of superimposed spikes and complete decomposition is preceded by a preprocessing phase. The preprocessing consists of: 1) filtering the signal to enhance spike train activities, 2) segmenting the filtered signal into temporal intervals containing spikes [START_REF] Mcgill | Automatic decomposition of the clinical electromyogram[END_REF], [START_REF] Mcgill | EMGLAB: An interactive EMG decomposition program[END_REF], 3) classifying isolated individual spikes (those not overlapped temporally with others) to determine the number of units I and the approximate spike shape for each class h (0) i . The preprocessing steps described above have robust solutions described in the literature. The thresholding method proposed in [START_REF] Mcgill | Automatic decomposition of the clinical electromyogram[END_REF] was adopted in this study for the segmentation. The classification of isolated action potentials can be performed with several techniques. For example, the canonically registered discrete Fourier transform (CRDFT) [START_REF] Mcgill | High-resolution alignment of sampled waveforms[END_REF] or non-parametric Bayesian estimation (NPB) [START_REF] Wood | A nonparametric bayesian alternative to spike sorting[END_REF], have been proven to be effective methods for this purpose [START_REF] Ge | Déconvolution impulsionnelle multi-source. Application aux signaux électromyographiques[END_REF]. These solutions for segmentation and classification of isolated spikes have excellent performance and do not need substantial improvements. Therefore, the methods described in [START_REF] Mcgill | Automatic decomposition of the clinical electromyogram[END_REF] and [START_REF] Mcgill | High-resolution alignment of sampled waveforms[END_REF] have been used in this study for preprocessing. On the contrary, the decomposition of temporally overlapped spikes poses more challenges and is the main contribution of this study, as described in the following.

IV. DECOMPOSITION OF OVERLAPPED POTENTIALS

USING MCMC Markov chain Monte Carlo (MCMC) methods are a class of stochastic simulation algorithms used to construct a Markov chain whose stationary distribution is invariant and converges to the desired distribution after a number of iterations (burn-in period). The estimator is then calculated from valid samples (i.e., those after the burn-in period) and its quality improves as a function of the sample population according to the Monte Carlo principle. The Metropolis-Hastings (MH) and Gibbs algorithms are among the most classic algorithms in the MCMC family, while Reversible Jump MCMC (RJM-CMC) [START_REF] Green | Reversible jump MCMC computation and Bayesian model determinaion[END_REF] further extends the application field by including variable dimension problems. These algorithms constitute widely used numerical tools in the field of Bayesian statistics and computational physics [START_REF] Robert | Monte Carlo Statistical Methods[END_REF].

In the multiunit spike sorting context, the joint posterior law in Eq. ( 7) is considered as the distribution of interest. The Markov chain is generated using a Gibbs sampler, by re-sampling iteratively each parameter in Θ according to its posterior conditional law derived from Eq. ( 7) while fixing the other parameters. The choice of conjugate prior laws on

{(m i , σ 2 i , h i , σ 2 ǫ ) i , σ 2 ǫ } facilitates
the re-sampling on continuous parameters (steps ( 6)- [START_REF] Takahashi | Automatic sorting for multineuronal activity recorded with tetrodes in the presence of overlapping spikes[END_REF] of Tab. I). The conditional laws to sample steps ( 6)-( 10) are detailed in [START_REF] Ge | Unsupervised Bayesian Decomposition of Multi-Unit EMG Recordings using Tabu Search[END_REF]. In particular, the conditional law for m i and σ 2 i are, respectively, Gaussian and Inverse Gamma if the two approximations (5) are satisfied.

The difficulty lies in the sampling of the discrete parameters

(x i ) (k) i 1 . We denote C k as the set of firing instants (x i ) (k) i
of all units in the k-th segment that satisfy the refractory period condition with respect to

(x i ) (-k) i
. Thus, each element in C k yields a non-zero joint posterior probability according to Eq. ( 4). Let us roughly evaluate the cardinality of C k to stress the computational burden of sampling

(x i ) (k) i .
Even in a very favorable case for which at most 4 units are active in a segment of duration of 8ms sampled at 10kHz, we would have 81 4 < |C k | < 2 500 . The inferior bound corresponds to a subspace of C k containing at most one impulse for each unit in the segment. For segment lengths that vary between 8 and 60ms and number of units that varies between 4 and 8, the cardinal of |C k | is at least 4.30 × 10 7 and 1.68 × 10 22 in each case. Thus, sampling the conditional probabilities of (x i )

(k)
i by evaluating all probabilities in the set C k for each segment yields unrealistic computational load.

To solve this problem, we propose a Metropolis-Hastings algorithm summarized in Tab. I. The algorithm explores a subspace of C k at each iteration Its validity is shown by verifying that the corresponding Markov chain remains irreducible: all configurations of non-zero probabilities (∈ C k ) are explored with non-zero probabilities regardless of the initial state. This proof is provided in appendix. Consequently, the MCMC algorithm in Tab. I is a valid stochastic simulation algorithm, of Metropolis-Hastings within Gibbs type.

repeat for k = 1 . . . K do % Integrate (si)i analytically Sample (xi) (k) i ∼ P ((xi) (k) i | Θ \ {(si)i, (xi) (k) i }, z) % MH end for Sample (si)i ∼ P ((si)i | Θ \ (si)i, z) Sample (hi)i ∼ P ((hi)i | Θ \ (hi)i, z) Sample (mi)i ∼ P ((mi)i | Θ \ (mi)i) Sample (σ 2 i )i ∼ P ((σ 2 i )i | Θ \ (σ 2 i )i) Sample (σ 2 s i )i ∼ P ((σ 2 s i )i | Θ \ (σ 2 s i )i) Sample σ 2 ǫ ∼ P (σ 2 ǫ | Θ \ σ 2 ǫ , z) until number of iteration reached
In what follows, Sections IV-A and IV-B enter into more details about the sampling of (x i ) (k) i (first step in Tab. I): they respectively focus on its marginal conditional law and on the implementation of the Metropolis-Hastings step.

A. Marginal conditional law

We first rewrite the convolution sum in the following matrix form:

I i=1 h i * s i = HS where H = [H 1 , . . . , H I ] is composed of convolution matrices such that h i * s i = H i s i ,
and S = [s 1 ; . . . ; s I ] represents a N I × 1 column vector by vertical concatenation. In the k-th segment, the data generation can be written:

z (k) = H (k) S (k) + ǫ (k)
where z (k) is the signal in the k-th segment;

H (k) = [H (k) 1 , . . . , H (k) I ] and S (k) = [s (k) 1 ; . . . ; s (k)
I ] denote the corresponding submatrix of H and subvector of S.

From Eq. ( 7), we obtain for each segment k the conditional law of (s i , x i ) (k) i | rest. In the following, the term rest stands for the set {Θ, z}, except the concerned parameters.

P (s i , x i ) (k) i | rest ∝ P z (k) | (s (k) i , h i ) i , σ 2 ǫ i P s (k) i | x (k) i , σ 2 si ) P (x (k) i | m i , σ 2 i , x (-k) i ∝ exp - 1 2σ 2 ǫ z (k) -Ga 1 2 |V| 1 2 δ(a 0 ) exp - 1 2 a 1 -1 t V a 1 -1 i P x (k) i | m i , σ 2 i , x (-k) i (8)
where the k-th segment specific matrices are noted without the superscript (k) for simplicity:

                 G = (H (k) 1 ) x (k) 1 , . . . , (H (k) I ) x (k) I V = diag [σ -2 s1 , . . . , σ -2 s1 n (k) 1 , . . . , σ -2 sI , . . . , σ -2 sI n (k) I ] a 1 = (s (k) 1 ) x (k) 1 , . . . , (s (k) I ) x (k) I a 0 = S (k) \ a 1
δ(•) denotes the multi-dimensional Dirac function, and n

(k) i = dim(x (k) i )
is the number of impulses of unit #i in the segment k. The vector 1 is a column vector of size i n (k) i while a 1 and a 0 are, respectively, the nonzero spike magnitudes of all units in the segment and its complement. {a 1 , a 0 } constitutes thus a permutation of spike trains S (k) .

Let us remark that matrices G and V are functions of

(x i ) (k) i , (h i , σ 2 si ) i , but not of (s i ) (k) i .
It is thus possible to integrate out first a 0 , and then a 1 in Eq. ( 8) since it is merely the product of two exponential terms and thus remains Gaussian.

The conditional probability of (x i )

(k) i
after marginalization of (s i ) (k) i from Eq. ( 8) can be written:

P (x i ) (k) i | rest \ (s i ) (k) i ∝ Σ 1 2 exp 1 2 m t Σ -1 m - 1 2 1 t V1 i σ -n (k) i si P x (k) i | m i , σ 2 i , x (-k) i (9)
where

Σ -1 = 1 σ 2 ǫ (G) t G + V m = Σ 1 σ 2 ǫ (G) t z (k) + V1
Σ -1 and m are, respectively, the variance and mean of the Gaussian vector a 1 . Let us remark that the marginal law in Eq. ( 9) involves Σ -1 and m. In the special case where (

i ) i = ∅, Σ -1 and m cannot be defined. To obtain the marginal distribution, we first rewrite Eq. ( 8):

P (s (k) i ) i = 0, (x i ) (k) i = ∅ | rest ∝ exp - z (k) 2 2σ ǫ δ (s (k) i ) i i P x (k) i = ∅ | m i , σ 2 i , x (-k) i
.

Integrating out (s (k)

i ) i and considering the common exponential factor w.r.t. cases where (x i ) (k) i = ∅, its marginal probability then can be written:

P (x i ) (k) i =∅ | rest \ (s i ) (k) i ∝ i P x (k) i = ∅ | m i , σ 2 i , x (-k) i (10)
Eqs. ( 9) and ( 10) describe the marginal conditional probability of (x i )

(k) i
up to a normalization factor. The term P (x

(k) i | m i , σ 2 i , x (-k) i
) can be directly derived from Eq. ( 4) and measures the regularity of the binary sequences (1 xi ) i for each spike train while the Gaussian variability of magnitudes of the spikes is expressed in the first two terms in Eq. [START_REF] Mcgill | EMGLAB: An interactive EMG decomposition program[END_REF].

We recall that though the marginal conditional probability of (x i ) (k) i can be analytically expressed, its combinatorial space makes it a hard problem either to maximize (as proposed in [START_REF] Ge | Unsupervised Bayesian Decomposition of Multi-Unit EMG Recordings using Tabu Search[END_REF] using a Tabu search) or to simulate the distribution using the MCMC approach. The next section proposes an MCMC algorithm that does not sample directly according to the conditional probability, but rather iteratively on local subspaces.

B. Metropolis-Hastings step

We propose here a reversible Metropolis-Hastings step to avoid the evaluation of all probabilities of (x i )

(k) i ∈ C k .
The strategy is to explore a reasonable subspace of C k per iteration while the exploration of the whole space C k is mathematically guaranteed in the long run. We note that this strategy can be adapted to the full model that includes Gaussian variability in spike magnitudes.

Let ω(u) ⊂ C k be the local subspace (to be specified hereafter) to explore for each iteration and thus contains accessible configurations from the current configuration, noted by u

= [1 x (k) 1 , . . . , 1 x (k) I ] of length I • dim(Seg k ).
It is constructed by concatenating spike trains of all units in the segment k. In Tab. II, P (u) denotes the conditional probability as specified by Eq. ( 9) and ( 10) for a particular configuration u, F (u) the sum of probabilities of the configurations in ω(u), i.e., F (u) = y∈ω(u) P (y), and q(u → u + ) the instrumental law or the probability of proposing u + given the current state u. Note that u + ∈ω(u) q(u → u + ) = 1.

TABLE II METROPOLIS-HASTINGS STEP TO SAMPLE (xi)

P (u)K(u → u + ) = P (u + )K(u + → u),
where K(u → u + ) denotes the transition kernel derived from the Metropolis-Hastings step of Tab. II and equals to:

P (u + ) F (u) ρ(u → u + ) + r(u)δ u (u + ) u + ∈ ω(u) 0 otherwise ( 11
)
where r(u) = u + ∈ω(u)

P (u + ) F (u) 1 -ρ(u → u + )
5) the Markov chain is irreducible: it is capable of exploring the entire space C k (see demonstration in Appendix); 6) the Markov chain is also aperiodic since the kernel (11) satisfies:

∀u, P (u) > 0 ⇒ K(u → u) > 0.
by verifying that P (u) F (u) > 0 and ρ(u → u) = 1. It is thus possible to have two consecutive samples that are identical, and such a Markov chain cannot be periodic.

7) the complexity per iteration (measured by |ω(u)|)

remains linear with respect to to the segment length dim(Seg k ) and the number of units I. In conclusion, the Markov chain generated from the proposed algorithm is irreducible, reversible and aperiodic. Therefore the only equilibrium distribution is P (u) (which means from any initial state u, K n (u, •)

D → P (•) [27, Theorem 1]).
The numerical implementation of the spike decomposition algorithm is summarized in Tab. III.

C. Bayesian estimators

With the convergence property of the Markov chain, we thus obtain for each parameter in Θ a population of random samples distributed according to its marginal posterior law P (• | z). For continuous parameters, i.e., Θ \ (x i , s i ) i , a minimum mean square error estimator (MMSE) optimizes the following criterion:

y = min y0 y -y 0 2 P (y | z)dy = min y0 E P (y | z) y -y 0 2 = E P (y | z) (y)
According to the Monte Carlo principle, this estimator is approximated by:

y = E P (y | z) (y) = yP (y | z)dy ≈ N j=N0+1 y j ,
where {y j } j=1,...,N are simulated samples, N 0 and N represent, respectively, the number of burn-in iterations and the total number of iterations. In our tests, we fixed N = 2N 0 = 200. The same estimator, however cannot be applied on the discrete parameters (x i ) i since averaged firing instants do not have any physical interpretation. Note that different samples may contain vectors of (x i ) i of different

]

Evaluate P (u + ), u + ∈ ω(u) in Eq. ( 9) (10)

Fu ← u + ∈ω(u) P (u + ) Propose u * ∼ P (u + )/Fu, u + ∈ ω(u) Fu * ← u + ∈ω(u * ) P (u + ) Accept u * with probability ρ = min {1, Fu/Fu * } % ----Magnitudes sampling------ [s (k) 1 , . . . , s (k) I ] ← 0 if (x (k) i )i = ∅ then G ← (H (k) 1 ) x (k) 1 , . . . , (H (k) I ) x (k) I d ← [σ -2 s 1 ones(n (k) 1 , 1), . . . , σ -2 s I ones(n (k) I , 1)] V ← diag(d) Σ -1 ← 1 σ 2 ǫ G t G + V Q ← chol(Σ -1 ) % Cholesky factor of Σ -1 b ← randn( i n (k) i , 1) % iid Gaussian vector % a1 ∼ N (m, Σ), Eq. (9) a1 ← Q -1 (Q -t (G t z (k) /σ 2 ǫ + V • 1) + b) [(s (k) 1 ) x (k) 1 , . . . , (s (k) I ) x (k) I ] ← a1 end if end for % -Sampling of remaining parameters [1] - h = [h1; . . . ; hI ], h ∼ P ((hi)i | rest) For each mi, mi ∼ P ((mi)i | rest) For each σ 2 i , σ 2 i ∼ P ((σ 2 i )i | rest) For each σ 2 s i , σ 2 s i ∼ P ((σ 2 s i )i | rest) σ 2 ǫ ∼ P (σ 2 ǫ | rest) until number of iteration reached
dimensions since the number of detected spikes are not fixed. An alternative is to estimate binary sequences (1 xi ) i using a marginal component MAP detector, i.e., the marginal majority vote for each instant (component).

However, considering only one component marginally at a time may result in counter-intuitive estimators. Suppose for example that in an interval J , 2 < dim(J ) < T R , the probability that there is one and only one detected event is 1 ( j∈J P (s i,j = 0) = 1) for a particular sample population. Further suppose that no one event is more likely than the sum of the others (in terms of probability), and thus P (s i,j = 0) < 0.5 for all j ∈ J . From these assumptions, it would then follow that s i,j = 0 for all j ∈ J , using the marginal component MAP detector. The result is clearly counterintuitive because the probability that there is no event in J is zero.

Recently, Kail et al. [START_REF] Kail | A blind Monte Carlo detection-estimation method for optical coherence tomography[END_REF] proposed a MAP block detector. The main idea is to perform the majority vote on a block of the binary variables in 1 xi . The number of binary combinations in each block should be much smaller than the number of available samples (the Markov chain length) while the length of the block is big enough to avoid false negatives using the marginal component MAP detector due to spike shifts in the sample population. The same strategy is adopted in this study for the estimation of the discrete parameters (x i ) i .

V. VALIDATION ON EXPERIMENTAL AND SIMULATED

SIGNALS

A. Experimental signals

The method proposed can be applied to several spike sorting applications, such as intraneural or intramuscular recordings. The assumptions of the model include a refractory period, which can be set depending on the application, and the cost associated to the firing irregularity. Note that the latter does not imply that irregular firings cannot be identified but simply that a solution with very irregular firings is less favored than one with more regular firings.

The validation of the method was performed on a set of experimental intramuscular EMG signals, where the spikes represent the activation of spinal motor neurons. These signals were chosen for direct comparison between the proposed method with the Tabu search method proposed in [START_REF] Ge | Unsupervised Bayesian Decomposition of Multi-Unit EMG Recordings using Tabu Search[END_REF]. The same experimental signals as described in [START_REF] Ge | Unsupervised Bayesian Decomposition of Multi-Unit EMG Recordings using Tabu Search[END_REF] were thus used for this study.

Briefly, the experimental signals were recorded from the abductor digiti minimi muscle of five healthy men (age, mean ± SD = 25.3 ± 4.5 yr) with a pair of wire electrodes made of Teflon coated stainless steel (A-M Systems, Carlsborg, WA, USA; diameter 50µm) inserted into the muscle with a 25 G needle. The intramuscular EMG signals were amplified and provided one bipolar recording (Counterpoint EMG, DANTEC Medical, Skovlunde, Denmark) that was band-pass filtered (500 Hz-5 kHz) and sampled at 10 kHz. A reference electrode was placed around the wrist. The decomposition method was evaluated on intervals of 20 s of signals recorded during isometric contractions at 10% of the maximal voluntary contraction (MVC) force. The validation was performed by comparing the results of the proposed methods with those provided as reference results by manual decomposition of an expert operator using the EMGLAB tool. The proposed method was applied in a fully automatic way, without any intervention by the operator.

B. Simulated signals

The experimental signals were recorded in conditions where the firing patterns were expected to be rather regular (isometric contractions). In addition to experimental signals, the method was also applied to a set of simulated signals. This set of simulations was used to test the performance of the algorithm in cases of very irregular spike firings. The simulated signals were generated with the intramuscular EMG model proposed in [START_REF] Farina | A model for the generation of synthetic intramuscular EMG signals to test decomposition algorithms[END_REF] and also previously applied for the validation of the method proposed in [START_REF] Ge | Unsupervised Bayesian Decomposition of Multi-Unit EMG Recordings using Tabu Search[END_REF]. The simulated signals were obtained by an imposed coefficient of ISI variation of 60%, corresponding to very irregular firings. The other parameters of the model were set as in a previous study [START_REF] Ge | Unsupervised Bayesian Decomposition of Multi-Unit EMG Recordings using Tabu Search[END_REF], including the level of noise that corresponded to a SNR of 10dB). In total, five signals were simulated with the same set of parameters. These signals differed from each other because the shape of the action potentials were varied within the library of shapes of the model [START_REF] Farina | A model for the generation of synthetic intramuscular EMG signals to test decomposition algorithms[END_REF] and because the spike trains were generated randomly, according to the imposed statistics.

C. Performance index

The performance index for each spike source was defined as [START_REF] Luca | Decomposition of Surface EMG Signals[END_REF]:

A(i) = N Dis (i) -N FP (i) -N FN (i) N Dis (i) × 100%
where N Dis (i) are spike numbers of the unit #i detected by EMGLAB for the experimental signals or simulated by the model for the synthetic signals. N FN (i) and N FP (i) are respectively the false negatives and the false positives with respect to the reference (EMGLAB or the model). A detected impulse was considered a correct one if within a 1ms-window centered on the instant of reference for the same unit. The global performance of the decomposition was then measured by the following index [START_REF] Luca | Decomposition of Surface EMG Signals[END_REF]:

A = 1 I I i=1 A(i)
For the experimental signals, discharge patterns for which σ i /m i > 0.3 were excluded from the analysis, since the ISI variability of motor unit spike trains is usually lower than this limit in the experimental conditions [START_REF] Moritz | Discharge Rate Variability Influences the Variation in Force Fluctuations Across the Working Range of a Hand Muscle[END_REF], [START_REF] Fuglevand | Models of recruitment and rate coding organization in motor-unit pools[END_REF]. Therefore, among the discharge patterns obtained by the application of the proposed method to experimental signals, only those satisfying the condition σ i /m i < 0.3 were selected for analysis and the global performance criterion A was thus averaged on the A(i) of the validated units. This a-posteriori selection of sources is not necessary in other applications in which the firing instants are less regular. The decomposition results of simulated firings provide a test of the method in conditions of very irregular spike firing. The a-posteriori selection of sources was not applied to the simulated signals and results are reported as average over all motor units. Note that the ISI variability of the simulated spike trains is doubled compared with the upper limit imposed to the experimental spike train decomposition.

D. Results of decomposition

The constant parameters of the Bayesian model (Fig. 2) and the Markov chain initializations are listed in Tab. IV. (h (0) i ) i , σǫ are results of the preprocessing and we fixed σ h = 0.1 max(h) in the tests. αi βi αs βs as bs 100ms (30ms

) 2 1 1 1 1 1 1 mi σ 2 i si xi σ 2 s i hi σǫ 100ms (30ms) 2 0 ∅ (0.15) 2 h (0) i σǫ
Fig. 3 shows an example of raw experimental signal and the decomposition of two adjacent segments with overlapped potentials using the proposed method. In this example, the fully automatic method proposed provided the same result as that obtained by the reference decomposition tool used by an expert operator. The numbers on top of the raw signal represent the source labels as identified by the automatic decomposition with the proposed method (top numbers) and the decomposition with interaction of an expert operator using EMGLAB (lower numbers). The labelling of sources in this example led to perfect matching with the reference result. The raw (solid line) and the reconstructed signal (dashed line) are shown in the upper panel whereas the residual error is plotted in the lower panel.

Fig. 4 shows an example of decomposition of simulated signals to better represent the ability of the method to identify complex superimpositions of action potentials. In this example, four overlapping action potentials are present in two segments. The shape of action potential #4, which is isolated in the second segment, is fully cancelled by overlaps of other action potentials in the first segment. Similarly, action potential #1 was relatively isolated in the first segment and cancelled in the second. Despite these complex overlaps, the automatic method was able to identify all action potentials in the two segments correctly. In this example, the reference decomposition result was provided by the simulation tool and is thus exact. ) and the number of validated spike trains according to the a-posteriori criterion on the regularity. Because the main challenge in the decomposition is dealing with overlapped action potentials [START_REF] Lewicki | A review of methods for spike sorting: the detection and classification of neural action potentials[END_REF], especially when more than two action potentials are superimposed, in Tab.V we also report the percentage of action potentials that are overlapped with others. Therefore, in Tab. V, % overlap and % overlap ≥ 3 indicates the percentage of overlapped action potentials (over the total number of action potentials) and the percentage of action potentials in superposition with at least other two action potentials (superposition of at least 3 action potentials), respectively. The algorithm runtime complexity was measured by computational time in seconds per iteration for the decomposition of each second of the experimental signal z. In order to decompose EMG signals of sufficient length (≈ 20 s to get reliable statistics) and upon several experimental signals, the Markov chain was systematically run for 200 iterations. The estimation of the discharge patterns was then obtained by a majority vote on each temporal window (of length = T R ), whose advantages are discussed in [START_REF] Kail | A blind Monte Carlo detection-estimation method for optical coherence tomography[END_REF].

The decomposition accuracy of the proposed method, as evaluated with the global performance index A, was improved (by approximately 3%) compared to that obtained by the joint maximization by Tabu search [START_REF] Ge | Unsupervised Bayesian Decomposition of Multi-Unit EMG Recordings using Tabu Search[END_REF] (see Tab. V). Most importantly, more spike trains were validated with the proposed method than with the Tabu search. For all tested methods in Tab. V, the detection errors (both false positives and false negatives) occurred less often in segments with single action potentials, while most errors were found in segments with overlaps, as it was expected. Therefore, the accuracy of the method on identifying overlapped action potentials was also computed. Over the entire number of overlapped action potentials, the accuracy was 89% and over all action potentials in superimpositions of at least three action potentials, the accuracy was 85%. Note that by using the MCMC method, the decomposition results would not vary for the same signal even if the latter is analyzed several times with different initialization values. This convergence property (independent of initialization values) derives directly from the convergence in the distribution sense of the MCMC method (see discussion in Section IV).

Tab. V also shows the decomposition results obtained by the proposed MCMC algorithm with the constant spike magnitude model for each source, i.e., without including the variable amplitude model. The same expert results were used in each case and the "#Sources" row reports the number of detected sources. The decomposition performance degraded using this simplified model compared to the model with variable magnitudes: the latter have achieved better results in both extracting more valid trains and better accuracy calculated upon these valid trains. The background noise variance estimation σ 2 ǫ can also be considered a good indicator of decomposition quality, since it measures the difference between the observed signal and its reconstruction from the detected discharge patterns and the estimated spikes. Indeed, σ 2 ǫ is sampled according to its conditional law IG(ᾱ s , βs ), where

α s = α s + N 2 , β s = β s + 1 2 z - i s i * h i 2
and the mean value decreases as zi s i * h i 2 decreases. Tab. VI shows σ 2 ǫ from the preprocessing (shown as reference), after full decomposition with optimization by Tabu search [START_REF] Ge | Unsupervised Bayesian Decomposition of Multi-Unit EMG Recordings using Tabu Search[END_REF] and with the proposed MCMC approach with and without spike magnitude variability. The lower value of residual noise level with the proposed approach (Tab. VI) is in agreement with the improvement of A by the MCMC method and the magnitude variability model that adjusts the energy for each detected spike. The lower error likely represents a better fit of the model w.r.t. the data rather than overfitting the noises since the reduction in residual error is also accompanied by superior accuracy, as shown in Tab. V. Finally, simulated signals with high variability in firing were used to evaluate the proposed decomposition method. One example of decomposition of simulated signals is presented in Fig. 4, as commented above. The average accuracy (over all motor units and over the 5 signals simulated) obtained for these signals was 83 ± 5%. The accuracy calculated for overlapped action potentials in the simulations was 80 ± 3%.

VI. DISCUSSION AND CONCLUSION

We have proposed a novel method for fully automatic spike sorting in multiunit recordings. The main difficulty in the problem of spike sorting is the resolution of superimposed spikes. The isolated spikes can indeed be detected and clustered with robust techniques, that have high performance. On the contrary, the global optimization of overlapped action potentials is a nondeterministic polynomial-type (NP) hard problem, and thus cannot be solved by polynomial complexity algorithms. Therefore, methods for full spike sorting differ mainly for their performance when separating overlapped action potentials. We have previously proposed a method for the solution of overlapped spikes in neural or intramuscular recordings based on a Bayesian model coupled with a MAP estimator and the Tabu search technique [START_REF] Ge | Unsupervised Bayesian Decomposition of Multi-Unit EMG Recordings using Tabu Search[END_REF]. In the present study, we have expanded the signal model to allow for a variability in the spike magnitudes in the same source. Moreover, we have proposed a new method for the solution of the spike sorting problem based on MCMC simulation. The convergence property of this method was proven theoretically (Appendix), whereas the Tabu search algorithm described in [START_REF] Ge | Unsupervised Bayesian Decomposition of Multi-Unit EMG Recordings using Tabu Search[END_REF] is heuristic. The performance of this approach is further improved by including variable magnitudes in the spike model.

The model is based on the assumption of a refractory period in neural firing, which implies to set a minimum time lapse between two consecutive discharges of the same source. Further, the model also favors regular firings with respect to irregular firings, but without constraining the solution space to regular firings. Accordingly, the method was tested on simulated signals with very large firing variability and proved to be relatively accurate also in those conditions. As it was expected, the accuracy was lower in case of irregular than for regular firings but the method could still be successfully applied in a fully automatic way with accuracy greater than 80% for an ISI variability of 60% (simulation results). It has to be noted that the imposed variability in the simulations was so large that the spike train regularity could not be used as a relevant information for reducing the solution space in the automatic method. The model is also based on the assumption of independence of the spike trains. We note however that this assumption is not always met since neural cells often fire more synchronously than by mere chance. However, the spike train independence assumption should not be viewed as a limitation of the approach since it does not imply that dependent spike trains cannot be identified, similarly to the discussion above for regular firings. The assumption only implies that the potential correlation between firing patterns is not taken into account in the prior laws, thus the correlation property is simply not exploited by the model. Accordingly, the method was proved successful in decomposing motor unit activities in a hand muscle, characterized by a high degree of motor unit synchronization [START_REF] Dideriksen | Comparison between the degree of motor unit short-term synchronization and recurrence quantification analysis of the surface EMG in two human muscles[END_REF].

A limitation of the method is that, although it includes changes in amplitude, it does not track changes in the shapes of the action potentials, which often occur over time. Some previous methods (e.g., [START_REF] Lefever | A procedure for decomposing the myoelectric signal into its constituent action potentials-part I: technique, theory, and implementation[END_REF], [START_REF] Studer | An Algorithm for Sequential Signal Estimation and System Identification for EMG Signals[END_REF]) have included this possibility, at least for slow changes over time. Inclusion of slow changes in shape in the current approach could be incorporated in the first preprocessing phase in which the individual motor unit action potential shapes are estimated from segments containing only one potential. However, this extension has not been tested in the current study.

In conclusion, a new fully automatic method for spike sorting of multiunit recordings has been derived and compared to a previously proposed solution on experimental intramuscular recordings. The method was proven to provide good accuracy in signal decomposition with relatively limited computational time even in the case of several spikes overlapped in time.

Fig. 1 .

 1 Fig. 1. Multiunit direct model.

Fig. 2 .

 2 Fig. 2. Directed acyclic graph representing the hierarchical Bayesian model. Unknown parameters are circled in the graph.

1 :

 1 Propose u → u + using the instrumental law:q(u → u + ) = 0 if u + / ∈ ω(u) P (u + )/F (u) otherwise 2: Accept u + with probability ρ(u → u + ) = min{1, F (u)/F (u + )} In this application, we define ω(u) = {v | |u-v| 1 ≤ 2 and |u| 1 -|v| 1 ≤ 1} with | • | 1 the L1-norm. The following conditions are verified: 1) |u -v| 1 = 0 ⇔ u = v, s.t. u ∈ ω(u); 2) |u -v| 1 = 1, adding or removing a spike; 3) |u -v| 1 = 2 and |u| 1 -|v| 1 ≤ 1, shiftingan existing spike in the same spike train or being replaced by a spike in another train; 4) u ∈ ω(v) ⇔ v ∈ ω(u), together with condition 1), assures the reversibility of the chain: for all pairs (u, u + ) ∈ C k ,

Fig. 3 .

 3 Fig.3. Example of EMG decomposition of two adjacent segments using the proposed approrach in comparison with the reference decomposition. In these two segments a total of 7 sources are active. The numbers on top of the raw signal represent the source labels as identified by the automatic decomposition with the proposed method (top numbers) and the decomposition with interaction of an expert operator using EMGLAB (lower numbers). The labelling of sources in this example led to perfect matching with the reference result. The raw (solid line) and the reconstructed signal (dashed line) are shown in the upper panel whereas the residual error is plotted in the lower panel.

Fig. 4 .

 4 Fig.4. Illustration of decomposition result of one segment containing up to four overlapped action potentials and the neighboring segment with three overlapped action potentials. Upper labels represent the source numbers identified by the proposed automatic method while lower labels those from the simulation tool. The residual error is also shown in the lower panel.

TABLE I STOCHASTIC

 I SIMULATION ALGORITHM BY MCMC.

TABLE III IMPLEMENTATION

 III OF THE MCMC ALGORITHM.

	repeat
	for k = 1, . . . , K do % for each segment k
	% ---Metropolis-Hastings step---------
	u ← [1 x (k) 1	, . . . , 1 x (k) I

TABLE V COMPARISON

 V OF DECOMPOSITION RESULTS ON THE 5 EXPERIMENTAL EMG SIGNALS (10%MVC, 20 S RECORDINGS).

	EMG #	1	2	3	4	5
	# Sources	5	5	8	8	4
	% overlap	81	84	91	88	82
	% overlap ≥ 3	63	61	74	69	65
			MCMC on extended model
	# Val	5	4	6	5	2
	A	91.1% 90.6% 90.5% 89.3% 94.4%
	Computing time (s.)	13.2	14.9	15.2	13.8	15.0
		MCMC on uniform magnitude model
	# Val	4	4	4	4	2
	A	87.1% 90.4% 92.3% 89.2% 90.7%
		Optimization using Tabu search [1]
	# Val	3	4	4	4	2
	A	89.4% 85.1% 87.7% 87.2% 92.5%

TABLE VI COMPARISON

 VI OF σ 2 ǫ BY OPTIMIZATION WITH TABU SEARCH AND BY MCMC ON THE 5 EXPERIMENTAL EMG SIGNALS (10% MVC, 20 S RECORDINGS).

	EMG #	1	2	3	4	5
	preprocessing	0.0012 0.0014 0.0012 0.0019 0.0011
	Tabu search	0.0023 0.0029 0.0030 0.0034 0.0026
	MCMC (full model) 0.0016 0.0027 0.0024 0.0030 0.0016
	MCMC (si = 1x i ) 0.0018 0.0027 0.0026 0.0033 0.0019

The superscript (k) and (-k) denote parameters in the k-th segment and parameters in all other segments, respectively.

APPENDIX

Proof: To show the irreducibility of the chain in Tab. I, it is sufficient to show that for all u ∈ C k , such that P (u) > 0, the probability of a particular trajectory composed of {u (0) = 0, . . . , u (i) , . . . u (L) = u}, where L = |u| 1 and |u (i)u (i-1) | 1 = 1, is non-zero.

Such a chain is constructed by adding a spike of a certain source in each iteration. It is then sufficient to show that for all i, K(u (i-1) → u (i) ) > 0, or equivalently,

It is evident that u (i) ∈ ω(u (i-1) ), and it follows that:

The inequality (A-1) is thus guaranteed by the fact that P (u (i) ) > 0 for all i. Since P (u) > 0 implies that