Cédric Allain
email: cedric.allain@inria.fr

Alexandre Gramfort
email: alexandre.gramfort@inria.fr

Thomas Moreau
email: thomas.moreau@inria.fr

DRIPP: DRIVEN POINT PROCESSES TO MODEL STIM-ULI INDUCED PATTERNS IN M/EEG SIGNALS

The quantitative analysis of non-invasive electrophysiology signals from electroencephalography (EEG) and magnetoencephalography (MEG) boils down to the identification of temporal patterns such as evoked responses, transient bursts of neural oscillations but also blinks or heartbeats for data cleaning. Several works have shown that these patterns can be extracted efficiently in an unsupervised way, e.g., using Convolutional Dictionary Learning. This leads to an event-based description of the data. Given these events, a natural question is to estimate how their occurrences are modulated by certain cognitive tasks and experimental manipulations. To address it, we propose a point process approach. While point processes have been used in neuroscience in the past, in particular for single cell recordings (spike trains), techniques such as Convolutional Dictionary Learning make them amenable to human studies based on EEG/MEG signals. We develop a novel statistical point process model -called driven temporal point processes (DriPP) -where the intensity function of the point process model is linked to a set of point processes corresponding to stimulation events. We derive a fast and principled expectation-maximization (EM) algorithm to estimate the parameters of this model. Simulations reveal that model parameters can be identified from long enough signals. Results on standard MEG datasets demonstrate that our methodology reveals event-related neural responses -both evoked and inducedand isolates non-task-specific temporal patterns.

INTRODUCTION

Statistical analysis of human neural recordings is at the core of modern neuroscience research. Thanks to non-invasive recording technologies such as electroencephalography (EEG) and magnetoencephalography (MEG), or invasive techniques such as electrocorticography (ECoG) and stereotactic EEG (sEEG), the ambition is to obtain a detailed quantitative description of neural signals at the millisecond timescale when human subjects perform different cognitive tasks [START_REF] Baillet | Magnetoencephalography for brain electrophysiology and imaging[END_REF]. During neuroscience experiments, human subjects are exposed to several external stimuli, and we are interested in knowing how these stimuli influence neural activity.

After pre-processing steps, such as filtering or Independent Component Analysis (ICA; [START_REF] Winkler | On the influence of high-pass filtering on ICA-based artifact reduction in EEG-ERP[END_REF] to remove artifacts, common techniques rely on epoch averaging -to highlight evoked responses -or time-frequency analysis to quantify power changes in certain frequency bands [START_REF] Cohen | Analyzing Neural Time Series Data: Theory and Practice[END_REF] -for induced responses. While these approaches have led to numerous neuroscience findings, it has also been criticized. Indeed, averaging tends to blur out the responses due to small jitters in time-locked responses, and the Fourier analysis of different frequency bands tends to neglect the harmonic structure of the signal, leading to the so-called "Fourier fallacy" [START_REF] Jasper | Charting the sea of brain waves[END_REF][START_REF] Jones | When brain rhythms aren't 'rhythmic': implication for their mechanisms and meaning[END_REF].

In so doing, one may conclude to a spurious correlation between components that have actually the same origin. Moreover, artifact removal using ICA requires a tedious step of selecting the correct components.

Driven by these drawbacks, a recent trend of work aims to go beyond these classical tools by isolating prototypical waveforms related to the stimuli in the signal [START_REF] Cole | Brain Oscillations and the Importance of Waveform Shape[END_REF][START_REF] Dupré La Tour | Multivariate convolutional sparse coding for electromagnetic brain signals[END_REF][START_REF] Donoghue | Parameterizing neural power spectra into periodic and aperiodic components[END_REF]. The core idea consists in decomposing neural signals as combinations of time-invariant patterns, which typically correspond to transient bursts of neural activity [START_REF] Sherman | Neural mechanisms of transient neocortical beta rhythms: Converging evidence from humans, computational modeling, monkeys, and mice[END_REF], or artifacts such as eye blinks or heartbeats. In machine learning, various unsupervised algorithms have been historically proposed to efficiently identify patterns and their locations from multivariate temporal signals or images [START_REF] Lewicki | Coding time-varying signals using sparse, shift-invariant representations[END_REF][START_REF] Jost | MoTIF: an efficient algorithm for learning translation invariant dictionaries[END_REF][START_REF] Heide | Fast and flexible convolutional sparse coding[END_REF][START_REF] Bristow | Fast convolutional sparse coding[END_REF]Wohlberg, 2016b), with applications such as audio classification [START_REF] Grosse | Shift-invariant sparse coding for audio classification[END_REF] or image inpainting (Wohlberg, 2016a). For neural signals in particular, several methods have been proposed to tackle this task, such as the sliding window matching (SWM; [START_REF] Gips | Discovering recurring patterns in electrophysiological recordings[END_REF], the learning of recurrent waveforms [START_REF] Brockmeier | Learning recurrent waveforms within EEGs[END_REF], adaptive waveform learning (AWL; [START_REF] Hitziger | Adaptive waveform learning: A framework for modeling variability in neurophysiological signals[END_REF] or convolutional dictionary learning (CDL; [START_REF] Jas | Learning the morphology of brain signals using alpha-stable convolutional sparse coding[END_REF][START_REF] Dupré La Tour | Multivariate convolutional sparse coding for electromagnetic brain signals[END_REF]. Equipped with such algorithms, the multivariate neural signals are then represented by a set of spatio-temporal patterns, called atoms, with their respective onsets, called activations. Out of all these methods, CDL has emerged as a convenient and efficient tool to extract patterns, in particular due to its ability to easily include physical priors for the patterns to recover. For example, for M/EEG data, Dupré la Tour et al. ( 2018) have proposed a CDL method which extracts atoms that appertain to electrical dipoles in the brain by imposing a rank-1 structure. While these methods output characteristic patterns and an event-based representation of the temporal dynamics, it is often tedious and requires a certain domain knowledge to quantify how stimuli affect the atoms' activations. Knowing such effects allows determining whether an atom is triggered by a specific type of stimulus, and if so, to quantify by how much, and with what latency.

As activations are random signals that consist of discrete events, a natural statistical framework is the one of temporal point processes (PP). PP have received a surge of interest in machine learning [START_REF] Bompaire | Machine learning based on Hawkes processes and stochastic optimization[END_REF]Shchur et al., 2020;Mei et al., 2020) with diverse applications in fields such as healthcare [START_REF] Lasko | Efficient inference of gaussian-process-modulated renewal processes with application to medical event data[END_REF][START_REF] Lian | A multitask point process predictive model[END_REF] or modelling of communities on social networks [START_REF] Long | NetCodec: Community Detection from Individual Activities[END_REF]. In neuroscience, PP have also been studied in the past, in particular to model single cell recordings and neural spike trains [START_REF] Truccolo | A point process framework for relating neural spiking activity to spiking history, neural ensemble, and extrinsic covariate effects[END_REF][START_REF] Okatan | Analyzing functional connectivity using a network likelihood model of ensemble neural spiking activity[END_REF][START_REF] Kim | A granger causality measure for point process models of ensemble neural spiking activity[END_REF][START_REF] Rad | Information rates and optimal decoding in large neural populations[END_REF], sometimes coupled with spatial statistics [START_REF] Pillow | Spatio-temporal correlations and visual signalling in a complete neuronal population[END_REF] or network models [START_REF] Galves | Modeling networks of spiking neurons as interacting processes with memory of variable length[END_REF]. However, existing models do not directly address our question, namely, the characterization of the influence of a deterministic PP -the stimuli onsets -on a stochastic one -the neural activations derived from M/EEG recordings.

In this paper, we propose a novel method -called driven point process (DriPP) -to model the activation probability for CDL. This method is inspired from Hawkes processes (HP; [START_REF] Hawkes | Point spectra of some mutually exciting point processes[END_REF], and models the intensity function of a stochastic process conditioned on the realization of a set of PP, called drivers, parametrized using truncated Gaussian kernels to better model latency effects in neural responses. The resulting process can capture the surge of activations associated to external events, thus providing a direct statistical characterization of how much a stimulus impacts the neural response, as well as the mean and standard deviation of the response's latency. We derive an efficient expectation-maximization (EM) based inference algorithm and show on synthetic data that it reliably estimates the model parameters, even in the context of MEG/EEG experiments with tens to hundreds of events at most. Finally, the evaluation of DriPP on the output of CDL for standard MEG datasets shows that it reveals neural responses linked to stimuli that can be mapped precisely both in time and in brain space. Our methodology offers a unified approach to decide if some waveforms extracted with CDL are unrelated to a cognitive task, such as artifacts or spontaneous brain activity, or if they are provoked by a stimulus -no matter if they are 'evoked' or 'induced' as more commonly described in the neuroscience literature [START_REF] Tallon-Baudry | Stimulus specificity of phase-locked and non-phase-locked 40 hz visual responses in human[END_REF]. While these different effects are commonly extracted using different analysis pipelines, DriPP simply reveals them as stimuli-induced neural responses using a single unified method, that does not require any manual tuning or selection.

DRIVEN TEMPORAL POINT PROCESS (DRIPP)

A temporal point process (PP) is a stochastic process whose realization consists of discrete events {t i } occurring in continuous time, t i ∈ R + [START_REF] Daley | An introduction to the theory of point processes[END_REF]. In case where the probability that an event occurs at time t only depends on the past events F t := {t i , t i < t}, PP are usually characterized through the conditional intensity function λ : R + → R + :

λ (t|F t ) := lim dt→0 P (N t+dt -N t = 1|F t ) dt , (1)
where N t := i≥1 1 ti≤t is the counting process associated to the PP. This function corresponds to the expected infinitesimal rate at which events are occurring at time t given the arrival times of past events prior to t [START_REF] Daley | An introduction to the theory of point processes[END_REF].

The proposed model DriPP is adapted from the Hawkes process (HP; Hawkes 1971), as the occurrence of a past event in the driver increases the likelihood of occurrence of activation events in the near future. However, here we suppose that the stochastic point process in our model of neural activations does not have the self-excitatory behavior characteristic of HP. Instead, the sources of activation in the DriPP model are either the drivers or some spontaneous background activity, but not its own previous activations. More specifically, in DriPP, the intensity function at time t between a stochastic process k -whose set of events is denoted A k -and a non-empty set of drivers P -whose events are denoted

T p := {t (p) 1 , . . . , t (p)
np }, p ∈ P -is composed of a baseline intensity µ k ≥ 0 and triggering kernels κ k,p : R + → R:

λ k,P (t) = µ k + p∈P i,t (p) i ≤t α k,p κ k,p t -t (p) i , (2)
where α k,p ≥ 0 is a coefficient which controls the relative importance of the driver p on the occurrence of events on the stochastic process k. Note that when the driver processes are known, the intensity function is deterministic, and thus corresponds to the intensity of an inhomogeneous Poisson process [START_REF] Daley | An introduction to the theory of point processes[END_REF]. The coefficient α k,p is set to be non-negative so that we only model excitatory effects, as events on the driver only increase the likelihood of occurrence of new events on the stochastic process. Inhibition effects are assumed non-existent. Figure 1 illustrates how events T p on the driver influence the intensity function after a short latency period.

A critical parametrization of this model is the choice of the triggering kernels κ k,p . To best model the latency, we use a parametric truncated normal distribution of mean m k,p ∈ R and standard deviation

σ k,p > 0, with support [a ; b] ⊂ R + , b > a. Namely, κ k,p (x) := κ (x; m k,p , σ k,p , a, b) = 1 σ k,p φ x-m k,p σ k,p Φ b-m k,p σ k,p -Φ a-m k,p σ k,p 1 a≤x≤b , (3)
where φ (resp. Φ) is the probability density function (resp. cumulative distribution function) of the standard normal distribution. This parametrization differs from the usual exponential kernel usually considered in HP, that captures responses with low latency. Note that the truncation values a, b ∈ R + are supposed independent of both the stochastic process and the drivers, hence they are similar for all kernel p ∈ P. Indeed, in the context of this paper, those values delimit the time interval during which a neuronal response might occur following an external stimulus. In other words, the interval [a ; b] denotes the range of possible latency values. In the following, we denote by T := T (k) the duration of the process k.

PARAMETERS INFERENCE WITH AN EM-BASED ALGORITHM

We propose to infer the model parameters Θ k,P = (µ k , α k,P , m k,P , σ k,P ), where we denote in bold the vector version of the parameter, i.e., x k,P = (x k,p ) p∈P , via maximum-likelihood using an EM-based algorithm [START_REF] Lewis | A nonparametric EM algorithm for multiscale Hawkes processes[END_REF][START_REF] Xu | Learning granger causality for Hawkes processes[END_REF]. The pseudocode of the algorithm is presented in Algorithm 1. The expectation-maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] is an iterative algorithm that allows to find the maximum likelihood estimates (MLE) of parameters in a probabilistic model when the latter depends on non-observable latent variables. First, from (2), we derive the negative log-likelihood of the model (see details in Appendix A.1):

L k,P (Θ k,P ) = µ k T + p∈P α k,p n p - t∈A k log   µk + p∈P i,t (p) i ≤t α k,p κ k,p t -t (p) i    . (4)
Expectation step For a given estimate, the E-step aims at computing the events' assignation, i.e., the probability that an event comes from either the kernel or the baseline intensity. At iteration n,

P (n) k (t) = µ (n) k λ (n) k,P (t) and ∀p ∈ P, P (n) p (t) = α (n) k,p i,t (p)
i ≤t κ

(n) k,p t -t (p) i λ (n) k,P (t) , (5)
where θ (n) denotes the value of the parameter θ at step n of the algorithm, and similarly, if f is a function of parameter θ, f (n) (x; θ) := f x; θ (n) . Note that ∀t ∈ [0 ; T ] , P p (t), we here derive succinctly the update for parameters µ and α:

∂ L k,P ∂µ (n) k Θ (n) k,P = 0 ⇔ T - t∈A k 1 λ (n) k,P (t) = 0 ⇔ T - t∈A k P (n) k (t) µ (n) k (t) = 0 ⇔ µ (n+1) k = 1 T t∈A k P (n) k (t) (6) ∂ L k,P ∂α (n) k,p Θ (n) k,P = 0 ⇔ n p - t∈A k P (n) p (t) α (n) k,p = 0 ⇔ α (n+1) k,p = π R + 1 n p t∈A k P (n) p (t) (7)
where π E (•) denotes the projection onto the set E. These two updates amount to maximizing the probabilities that the events assigned to the driver or the baseline stay assigned to the same generation process. Then, we give the update equations for m and σ, which corresponds to parametric estimates of each truncated Gaussian kernel parameter with events assigned to the kernel. Detailed computations are provided in Appendix A.2.

m (n+1) k,p = α (n) k,p t∈A k i,t (p) i ≤t t-t (p) i κ (n) k,p t-t (p) i λ (n) k,P (t) t∈A k P (n) p (t) -σ (n) k,p 2 C m m (n) k,p , σ (n) k,p , a, b C m (n) k,p , σ (n) k,p , a, b (8)
σ (n+1) k,p = π [ε ;+∞)      C m (n) k,p , σ (n) k,p , a, b C σ m (n) k,p , σ (n) k,p , a, b α (n) k,p t∈A k i,t (p) i ≤t t-t (p) i (t)-m (n) k,p 2 κ (n) k,p t-t (p) i λ (n) k,P (t) t∈A k P (n) p (t)      1/3 (9) where, C (m, σ, a, b) := b a exp -1 2 (u-m) 2 σ 2 du, C m (m, σ, a, b) := ∂ C ∂m (m, σ, a, b) and C σ (m, σ, a, b) := ∂ C ∂σ (m, σ, a, b).
Here ε > 0 is predetermined to ensure that σ remains strictly positive. In practice, we set ε such that we avoid the overfitting that can occur when the kernel's mass is too concentrated. Note that once the initial values of the parameters are determined, the EM algorithm is entirely deterministic.

Algorithm 1: EM-based algorithm input :A k , T p , a, b, T, N output :The estimated values for parameters µ, α, m and σ

Initialize µ (0) , α (0) , m (0) , σ (0) ; for i = 0, . . . , N -1 do if α (i) = 0 R #P then µ (i+1) = µ (MLE) ; break; end Define λ (i) ; Compute µ (i+1) , α (i+1) , m (i+1) , σ (i+1) ; end return µ (i+1) , α (i+1) , m (i+1) , σ (i+1)
Also, when the estimate of parameter m is too far from the kernel's support [a ; b], we are in a pathological case where EM is diverging due to indeterminacy between setting α = 0 and pushing m to infinity due to the discrete nature of our events. Thus, we consider that the stochastic process is not linked to the considered driver, and fall back to the MLE estimator. The algorithm is therefore stopped and we set α = 0 R #P .

It is worth noting that if ∀p ∈ P, α k,p = 0, then the intensity is reduced to its baseline, thus the negative log-likelihood is

L k,p (Θ k,p ) = µ k T -#A k log µ k ,
where #A denotes the cardinality of the set A. Thus, we can terminate the EM algorithm by directly computing the MLE for µ k , namely: µ

(MLE) k = #A k /T .

Initialization strategy

We propose a "smart start" initialization strategy, where parameters are initialized based on their role in the model. It reads:

µ (0) k = #A k -# p∈P D k,p T -λ p∈P t ∈Tp [t + a ; t + b] (10) α (0) k,p = #D k,p λ t ∈Tp [t + a ; t + b] -µ (0) k , ∀p ∈ P (11) m (0) k,p = 1 #D k,p d∈D k,p d and σ (0) k,p = 1 #D k,p d∈D k,p d -m (0) 2 , ∀p ∈ P , (12)
where λ (•) denotes the Lebesgue measure, and where

D k,p := {t -t (p) * (t), t ∈ A k } ∩ [a ; b]
is the set of all empirical delays possibly linked to the driver p, with t (p) * (t) := max{t , t ∈ T p , t ≤ t} denoting the timestamp of the last event on the driver p that occurred before time t. Here, the baseline intensity µ (0) is set to the average number of process' events that occur outside any kernel support, i.e., the events that are guaranteed to be exogenous or spontaneous. Similarly, the kernel intensity α (0) is computed as the increase in the average number of activations over the kernel support, compared to µ (0) . The initial guess for m (0) and σ (0) are obtained with their parametric estimates, considering that all event on the kernel support are assigned to the considered driver.

EXPERIMENTS

We evaluated our model on several experiments, using both synthetic and empirical MEG data. We used Python (Python Software Foundation, 2019) and its scientific libraries [START_REF] Virtanen | SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python[END_REF][START_REF] Hunter | Matplotlib: A 2d graphics environment[END_REF][START_REF] Harris | Array programming with NumPy[END_REF]. We relied on alphacsc for CDL with rank-1 constraints on MEG [START_REF] Dupré La Tour | Multivariate convolutional sparse coding for electromagnetic brain signals[END_REF] and we used MNE [START_REF] Gramfort | MEG and EEG data analysis with MNE-Python[END_REF] to load and manipulate the MEG datasets. Computations were run on CPU Intel(R) Xeon(R) E5-2699, with 44 physical cores.

Ground truth Estimated

Figure 2: True and estimated intensity functions following a driving event at time zero for two different kernels, on synthetic data. Left: "wide" kernel with σ = 0.2. Right: "sharp" kernel with σ = 0.05. Parameters used are T = 10000, P/S = 0.6. On synthetic data, the EM algorithm successfully retrieves the true values of parameters, for both shapes of kernels.

EVALUATION OF THE EM CONVERGENCE ON SYNTHETIC DATA

For a given number of drivers and a set of corresponding parameters Θ, we first generate the drivers' processes and then simulate the stochastic process for a pre-determined duration T . Each driver's timestamps are simulated as follows: given an interstimuli interval (ISI), a set of S = T ISI equidistant timestamps is generated -where • denotes the floor function. Then P timestamps are uniformly sampled without replacement from this set. In all our experiments, we fixed the ISI to 1 s for the "wide" kernel, and to 1.4 s for the "sharp" one. Finally, a one-dimensional non-homogeneous Poisson process is simulated following Lewis' thinning algorithm [START_REF] Lewis | Simulation of nonhomogeneous Poisson processes by thinning[END_REF], given the predefined intensity function λ and the drivers' timestamps.

Figure 2 illustrates the intensity function recovery with two drivers considered together: the first one has a "wide" kernel with standard deviation σ = 0.2 s, and the second one has a "sharp" kernel with σ = 0.05 s. Both kernels have support [0.03 s ; 0.8 s] and mean m = 0.4 s, the coefficients α are both set to 0.8 and the baseline intensity parameter µ to 0.8. We report 8 estimated intensities obtained from independent simulations of the processes -using T = 10000 s and P/S = 0.6 -that we plot over each one of the driver's kernel's support. The EM algorithm is run for 50 iterations using the "smart start" initialization strategy described in Section 2. Note that here, the randomness only comes from the data generation, as the EM algorithm uses a deterministic initialization. Figures demonstrate that the EM algorithm is able to successfully recover the parameters for both shapes of kernels.

To provide a quantitative evaluation of the parameters' recovery, we compute, for each driver p ∈ P, the ∞ norm between the intensity λ * computed with the true parameters Θ * p and the estimated intensity λ p with parameters Θ p :

λ * -λ p ∞ := max t∈[0 ;T ] µ * + α * p κ * p (t) -μ -αp κp (t) . ( 13
)
The rationale for using the ∞ norm is to ensure that errors during baseline and within the kernel support are given equal importance. Figure 3 presents the parameter recovery for the same scenario with varying P/S and T . To easily compare the EM performances on the two shapes of kernels, Figure 3 (resp. Figure A.1) reports the mean (resp. standard deviation) of the relative ∞ norm -that is the ∞ divided by the maximum of the true intensity λ * -computed for each of the driver over 30 repetitions with different realizations of the process. The results show that the more data are available, either due to a longer process duration (increase in T ) or due to a higher event density (increase in P/S), the better are the parameter estimates. The convergence appears to be almost linear in these two cases. Moreover, the average computation time for an EM algorithm in Figure 3 The accuracy of the EM estimates increases with longer and denser processes.

EVOKED AND INDUCED EFFECTS CHARACTERIZATION IN MEG DATA

Datasets Experiments on MEG data were run on two datasets from MNE Python package [START_REF] Gramfort | MNE software for processing MEG and EEG data[END_REF]2013): the sample dataset and the somatosensory (somato) dataset1 . These datasets were selected as they elicit two distinct types of event-related neural activations: evoked responses which are time locked to the onset of the driver process, and induced responses which exhibit random jitters.

Complementary experiments were performed on the larger Cam-CAN dataset [START_REF] Shafto | The Cambridge Centre for Ageing and Neuroscience (Cam-CAN) study protocol: a cross-sectional, lifespan, multidisciplinary examination of healthy cognitive ageing[END_REF] 2 . Presentation of the dataset, data pre-processing and obtained results on 3 subjects are presented in Appendix A.7. The presented results are self-determined as they exhibit, for each subject, the atoms that have the higher ratio α/µ. For all studied datasets, full results are presented in supplementary materials.

The sample dataset contains M/EEG recordings of a human subject presented with audio and visual stimuli. In this experiment, checkerboard patterns are presented to the subject in the left and right visual field, interspersed by tones to the left or right ear. The experiment lasts about 4.6 min and approximately 70 stimuli per type are presented to the subject. The interval between the stimuli is on average of 750 ms, all types combined, with a minimum of 593 ms. Occasionally, a smiley face is presented at the center of the visual field. The subject was asked to press a button with the right index finger as soon as possible after the appearance of the face. In the following, we are only interested in the four main stimuli types: auditory left, auditory right, visual left, and visual right. For the somato dataset, a human subject is scanned with MEG during 15 min, while 111 stimulations of his left median nerve were made. The minimum ISI is 7 s.

Experimental setting For both datasets, only the 204 gradiometer channels are analyzed. The signals are pre-processed using high-pass filtering at 2 Hz to remove slow drifts in the data, and are resampled to 150 Hz to limit the atom size in the CDL. CDL is computed using alphacsc (Dupré la Tour et al., 2018) with the GreedyCDL method. For the sample dataset, 40 atoms of duration 1 s each are extracted, and for the somato dataset, 20 atoms of duration 0.53 s are estimated. The extracted atoms' activations are binarized using a threshold of 6 × 10 -11 (resp. 1 × 10 -10 ) for sample (resp. somato), and the times of the events are shifted to make them correspond to the peak amplitude time of the atom. Then, for every atom, the intensity function is estimated using the EM-based algorithm with 400 iterations and the "smart start" initialization strategy. Kernels' truncation values are hyper-parameters for the EM and thus must be pre-determined. The upper truncation value b is chosen smaller than the minimum ISI. Here, we used in addition some previous domain knowledge to set coherent values for each dataset. Hence, for the sample (resp. somato) dataset, kernel support is fixed at [0.03 s ; 0.5 s] (resp. [0 s ; 2 s]). See Appendix A.4 for an analysis on how these hyperparameters influence on the obtained results presented below. are specifically handpicked to exhibit the usual artifacts, and the last two are selected as they have the two bigger ratios α/µ for their respective learned intensity functions. Even though the intensity is learned with the two stimuli conjointly, we plot the two corresponding "intensities at the kernel separately, i.e., ∀p ∈ P, ∀t ∈ [0 ; 0.5], we plot λ k,p (t), k = 0, 1, 2, 6.

Evoked responses in sample dataset

Spatial and temporal representations of atom 0 (resp. atom 1) indicate that it corresponds to the heartbeat (resp. the eye blink) artifact. These two atoms are thus expected not to be linked to any stimuli. This is confirmed by the shape of the intensities estimated with DriPP that is mostly flat, which indicates that the activation of these two atoms are independent of auditory and visual stimuli. Note that these two artifacts can also be recovered by an Independent Component Analysis (ICA), as shown in Figure A.5. Indeed, the cosine similarity between the spatial maps of the eye blink (resp. the heartbeat) artifact extracted with CDL and its corresponding component in ICA analysis is 99.58 % (resp. 99.78 %), as presented in Figure A.6. In contrast, by looking at the spatial and temporal patterns of atom 2 (resp. atom 6), it can be associated with an auditory (resp. visual) evoked response. Given the spatial topography of atom 2, we conclude to a bilateral auditory response and the peak transient temporal pattern suggests an evoked response that is confirmed by the estimated intensity function that contains a narrow peak around 85 ms post-stimulus. This is the M100 response -here the auditory one -well known in the MEG literature (its equivalent in EEG is the N100) [START_REF] Näätänen | The N1 wave of the human electric and magnetic response to sound: A review and an analysis of the component structure[END_REF]. The M100 is indeed a peak observed in the evoked response between 80 and 120 milliseconds after the onset of a stimulus in an adult population. Regarding atom 6, topography is right lateralized in the occipital region, suggesting a visual evoked response. This is confirmed by the intensity function estimated that reports a relationship between this atom and the visual stimuli. Here also, the intensity peak is narrow, which is characteristic of an evoked response. This reflects a right lateralized response along the right ventral visual stream in this subject. This may be connected to the P200, a peak of the electrical potential between 150 and 275 ms after a visual onset. Moreover, the intensities estimated with DriPP for the unrelated tasks are completely flat. We have α = 0, which indicates that atoms' activations are exogenous or spontaneous relatively to unrelated stimuli. For comparison, we present in Appendix A.5 similar results obtained with dedicated MEG data analysis tools, such as evoked responses and time-frequency plots. Figure 5: Spatial and temporal patterns of 3 atoms from somato dataset, and their respective estimated intensity functions following a somatosensory stimulus (cue at time = 0 s). The eye-blink artifact (atom 0) is not linked to the stimulus, and neither is the α-wave (atom 7). A somatosensory stimulus will induce a neural response similar to atom 2, with a mean latency of 1 s.

Induced response in somato dataset Results on the somato dataset are presented in Figure 5. Similar to the results on sample, spatial and temporal patterns of 3 handpicked atoms are plotted alongside the intensity functions obtained with DriPP. Thanks to their spatial and temporal patterns, and with some domain knowledge, it is possible to categorize these 3 atoms: atom 2 corresponds to a µ-wave located in the secondary somatosensory region (S2), atom 7 corresponds to an α-wave originating in the occipital visual areas, whereas atom 0 corresponds to the eye-blink artifact. As α-waves are spontaneous brain activity, they are not phase-locked to the stimuli. It is thus expected that atom 7 is not linked to the task, as confirmed by its estimated intensity function where α = 0.

For atom 2 -that corresponds to a µ-wave -, its respective intensity is nonflat with a broad peak close to 1 s, which characterizes an induced response. Moreover, similar to results on the sample dataset, we recover the eye-blink artifact that also has a flat intensity function. This allows us to be confident in the interpretation of the obtained results. Some other µ-wave atoms -atoms 1 and 4are presented in Figure A.10 in Appendix A.6. They have an estimated intensity similar to atom 2, i.e., non-flat with a broad peak close to 1 s. The usual time/frequency analysis reported in Figure A.9 exhibits the induced response of the µ-wave.

DISCUSSION

This work proposed a point process (PP) based approach specially designed to model how external stimuli can influence the occurrences of recurring spatio-temporal patterns, called atoms, extracted from M/EEG recordings using convolutional dictionary learning (CDL). The key advantage of the developed method is that by estimating few parameters (one baseline parameter and 3 parameters per considered driver), it provides a direct statistical characterization of when and how each stimulus is responsible for the occurrences of neural responses. Importantly, it can achieve this with relatively limited data which is well adapted to MEG/EEG experiments that last only a few minutes, hence leading to tens or hundreds of events at most. This work proposed an EM algorithm derived for a novel kernel function: the truncated Gaussian, which differs from the usual parametrization in PP models that capture immediate responses, e.g., with exponential kernels. As opposed to competing methods that can involve manual selection of task-related neural sources, DriPP offers a unified approach to extract waveforms and automatically select the ones that are likely to be triggered by the considered stimuli. Note however that DriPP has been developed based on a point process framework, which is event-based. When working with continuous stimuli, other techniques must be considered (e.g., STRF, spatio-temporal response functions; [START_REF] Drennan | Cortical tracking of complex sound envelopes: modeling the changes in response with intensity[END_REF]. Future work will explore the modeling of inhibitions effects (α < 0) frequently observed in neuroscience data. Yet, this could potentially degrade the interpretability of the model parameters, and it also requires a new definition of the intensity function as one needs to ensure it stays non-negative at all times.

A APPENDIX A.1 DETAILS OF THE NEGATIVE LOG-LIKELIHOOD COMPUTATION

We defined the negative log-likelihood for parameters Θ k,P = (µ k , α k,P , m k,P , σ k,P ) as follows:

L k,P (Θ k,P ) = T 0 λ k,P (t)dt - t∈A k log λ k,P (t) . ( 14
)
We show here in details that T 0 λ k,p (t)dt = µ k T + p∈P α k,p n p . Without loss of generality, we can first assume that ∀p ∈ P, t

(p) np + b ≤ T . Hence, ∀p ∈ P, ∀i = 1, . . . , n p , T 0 κ k,p t -t (p) i dt = 1 .
Thus, we have that

T 0 λ k,p (t) dt = T 0   µk + p∈P i,t (p) i <t α k,p κ k,p t -t (p) i    dt = µ k T + p∈P α k,p    t (p) i ∈Tp T 0 κ k,p t -t (p) i dt    = µ k T + p∈P α k,p n p .
Finally, we have that the negative log-likelihood writes

L k,P (Θ k,P ) = µ k T + p∈P α k,p n p - t∈A k log λ k,P (t) .

A.2 DETAILS OF EM-BASED ALGORITHM COMPUTATIONS

First, recall that the negative log-likelihood for parameters Θ k,P = (µ k , α k,P , m k,P , σ k,P ) is as follows:

L k,P (Θ k,P ) = µ k T + p∈P α k,p n p - t∈A k log   µk + p∈P i,t (p) i ≤t α k,p κ k,p t -t (p) i ; m k,p , σ k,p    .
Recall also that we defined P

P (n) k (t) = µ (n) k λ (n) k,p (t) and ∀p ∈ P, P (n) p (t) = α (n) k,p i,t (p) i ≤t κ (n) k,p t -t (p) i (t) λ (n) k,P (t) , (15)
where θ (n) denotes the value of the parameter θ at step n of the algorithm, and similarly, if f is a function of parameter θ, f (n)

(x; θ) := f (x; θ (n) ).
Below are details of calculation for parameters update for m and σ. We first rewrite the kernel function to have it under a form that simplifies further computations.

κ(x; m, σ, a, b) = 1 σ φ x-m σ Φ b-m σ -Φ a-m σ 1 a≤x≤b = exp -1 2 (x-m) 2 σ 2 C (m, σ, a, b) 1 a≤x≤b where C (m, σ, a, b) := b a exp - 1 2 (u -m) 2 σ 2 du . (16)
Hence, we can precompute its derivatives with respect to m and σ:

∂ ∂m κ(x; m, σ, a, b) = x -m σ 2 - C m (m, σ, a, b) C (m, σ, a, b) κ(x; m, σ, a, b) , (17)
and

∂ ∂σ κ(x; m, σ, a, b) = (x -m) 2 σ 3 - C σ (m, σ, a, b) C (m, σ, a, b) κ(x; m, σ, a, b) . ( 18
)
where

C m (m, σ, a, b) := ∂ C ∂m (m, σ, a, b) and C σ (m, σ, a, b) := ∂ C ∂σ (m, σ, a, b).
Update equation for m k,p

∂ L k,P ∂m k,p (Θ k,P ) = - t∈A k i,t (p) i ≤t t -t (p) i -m k,p σ 2 k,p - C m (m k,p , σ k,p , a, b) C (m k,p , σ k,p , a, b) α k,p κ k,p t -t (p) i λ k,P (t) = m k,p σ 2 k,p + C m (m k,p , σ k,p , a, b) C (m k,p , σ k,p , a, b) t∈A k P p (t) - α k,p σ 2 k,p t∈A k i,t (p) i ≤t t -t (p) i κ k,p t -t (p) i λ k,P (t) (19)
Hence, by canceling the previous derivative,

m (n+1) k,p = α (n) k,p t∈A k i,t (p) i ≤t t-t (p) i κ (n) k,p t-t (p) i λ (n)
k,P (t)

t∈A k P (n) p (t) -σ (n) k,p 2 C m m (n) k,p , σ (n) k,p , a, b C m (n) k,p , σ (n) k,p , a, b (20)
Update equation for σ k,p

∂ L k,P ∂σ k,p (Θ k,P ) = - t∈A k i,t (p) i ≤t    t -t (p) i -m k,p 2 σ 3 k,p - C σ (m k,p , σ k,p , a, b) C (m k,p , σ k,p , a, b)    α k,p κ k,p t -t (p) i λ k,P (t) = C σ (m k,p , σ k,p , a, b) C (m k,p , σ k,p , a, b) t∈A k P p (t) - α k,p σ 3 k,p t∈A k i,t (p) i ≤t t -t (p) i (t) -m k,p 2 κ k,p t -t (p) i λ k,P (t) (21)
Hence, by canceling the previous derivative,

σ (n+1) k,p =      C m (n) k,p , σ (n) k,p , a, b C σ m (n) k,p , σ (n) k,p , a, b α (n) k,p t∈A k i,t (p) i ≤t t-t (p) i (t)-m (n) k,p 2 κ (n) k,p t-t (p) i λ (n)
k,P (t)

t∈A k P (n) p (t)      1/3 (22)
Finally, to ensure that the σ coefficient stays strictly positive in order to avoid computational errors, we add a projection step onto [ε ; +∞), with ε > 0: σ In this section, we present figures that are complementary to the ones obtained on synthetic data, presented in subsection 4.1. In Figure A.1, we report the standard deviation associated with Figure 3, i.e., the standard deviation of the relative infinite norm for different process duration T and different proportion of events kept P/S, for the two kernel shapes.

(n+1) k,p = π [ε ;+∞) σ (n+1) k,p

A.4 IMPACT OF MODEL HYPERPARAMETER

In this section, we dwell on the analysis of how setting hyperparameter values may impact the obtained results, with the aim of determining whether it is possible to set these parameters using a general rule of thumb, without degrading previous results. More specifically, we will look at the impact of two hyperparameters on the results we obtained on the MNA sample dataset: the threshold value -applied to the atoms' activation values to binarized them, currently set at 6 × 10 -11 -, and the kernel support, currently set at [0.03 s ; 0.5 s] using previous and domain knowledge. To do so, we conducted two experiments on sample where each varies one hyperparameter. We plot the intensity function learned by DriPP for the same four atoms -namely, the two artifacts (heartbeat and eye blink) and the auditory and visual responses -, separately for the two stimuli (auditory and visual), similarly to Figure 4.

For the first experiment, presented in Figure A.3, we varied the threshold, expressed as a percentile, between 0 and 80. A threshold of 20 means that we only keep activations whose values are above the 20 % percentile computed over all strictly positive activations, i.e., the smaller the threshold is, the more activations are kept. The value of the threshold used in all other experiments is the 60 % percentile. One can observe that for the two artifacts (atoms 0 and 1), when the threshold gets smaller, the learned intensity functions get flatter, indicating that the stimulus has no influence on the atom activation. However, for the two others atoms, the effect of a smaller threshold is the opposite, as the intensity functions have a higher peak, indicating a bigger value of the α parameter, and thus strengthening the link stimulus-atom. Thus, the threshold value could be set to a small percentile and therefore computed without manual intervention, without degrading the current results.

For the second experiment, we now focus on the kernel truncation values a and b. Results are presented in Figure A.4. We set a = 0, as we did for somato and Cam-CAN datasets, and we vary b between 0.5 -the current value -and 10, a large value compared to the ISI. One can observe that for b = 0.5 or b = 1, the results are either unchanged (atoms 1, 2, and 6) or better (atom 0, as the intensity function is totally flat for this artifact), indicating that setting a = 0 and b close to the average ISI of 0.750 s does not hinder the results. However, when b is too large, the results degrade quickly, up to the point where all learned intensities are flat lines, indicating that our model does not find any link between the stimuli and the atoms. This is due to the fact that this hyperparameter is of great importance in the initialization step, as the greater b is, the more atom's activations are considered being on a kernel support. Thus, setting the upper truncation value to a value close to the average ISI seems to give reliable results.

A.5 USUAL M/EEG DATA ANALYSIS

We present in this section some results obtained using usual M/EEG data analysis, such as Independent Component Analysis (ICA), epoch averaging, or time/frequency analysis. First, on MNE sample dataset, we proceed to an ICA to manually identify usual artifacts. To do so, similarly as the CDL pre-processing, the raw signals are filtered (high-pass filter at 2 Hz), and 40 independent components are fitted. The two components 1 and 3, that we manually identify as corresponding to the eye blink and heartbeat artifacts, are presented in Figure A.5. In Figure A.6, we associate for each of the CDL atoms presented in Figure 4 the ICA component that has the maximum cosine similarity. One can observe that the artifact atoms and components are highly similar, suggesting that CDL and ICA have equal performance on the artifact detection. Note that this high similarity is only based on the spatial pattern. Indeed, ICA does not provide temporal waveforms for the atoms as well as their temporal onsets, contrarily to CDL.

However, for the auditory and visual response, the result is different. For the auditory one (atom 1), there is not really an ICA equivalent, as it is most correlated with the eye-blink ICA component.

Regarding the visual atom (atom 6), there is an ICA component that presents a high similarity. While the two related components correspond to neural sources on the occipital cortex, the atom 6 obtained with CSC is more right lateralized, suggesting a source in the right ventral visual stream. Note that unlike ICA, which recover full time courses for each source, CDL also provides the onset of the patterns, which we later use for automated identification of event related components. Finally, this demonstrates that CDL is a strong competitor to ICA for artifact identification, while simultaneously enabling to reveal evoked or induced neural responses in an automated way.

As mentioned, the ICA is commonly used to remove manually identified artifacts to reconstruct the original signals free of those artifacts. However, there are two drawbacks of this method, the first one being that some domain-related knowledge is needed in order to correctly identify the artifacts. The second drawbacks happen when the signal is reconstructed after the removal of certain components. Indeed, such a reconstruction will lead to a loss of information across all channels, as the artifacts are shared by all the sensors. Thanks to the Convolutional dictionary learning (CDL) that extracts the different artifacts directly from the raw data, our method does not suffer from these drawbacks. Finally, note that these methods that are commonly used in the M/EEG data analysis comprise a thoughtfully data pre-processing as well as a manual intervention requiring domain knowledge to identify both artifacts and evoked and induced responses. As the proposed method in this paper is composed of a unified pipeline, it is a significant gain. Indeed, DriPP is able to automatically isolate artifacts without prior removal of artifacts using ICA, as well as capture the diversity of latencies corresponding to induced responses.

Regarding the statistical significance of the link between a stimulus and a neural pattern, one can consider performing a statistical test, where one wishes to reject the null hypothesis of independence. We present in this section extra results obtained on the real dataset somato, that are complementary of Figure 5. The Cam-CAN dataset contains data of M/EEG recordings of 643 human subjects submitted to audio and visual stimuli. In this experiment, 120 bimodal audio/visual trials and eight unimodal trials -included to discourage strategic responding to one modality (four visual only and four auditory only) -are presented to each subject. For each bimodal trial, participants see two checkerboards presented to the left and right of a central fixation (34 milliseconds duration) and simultaneously hear a 300 milliseconds binaural tone at one of three frequencies (300, 600, or 1200 Hz, equal numbers of trials pseudorandomly ordered). For unimodal trials, participants either only hear a tone or see the checkerboards. For each trial, participants respond by pressing a button with their right index finger if they hear or see any stimuli [START_REF] Shafto | The Cambridge Centre for Ageing and Neuroscience (Cam-CAN) study protocol: a cross-sectional, lifespan, multidisciplinary examination of healthy cognitive ageing[END_REF]. For each subject, the experiment lasts less than 4 min.

The signals are pre-processed using low-pass filtering at 125 Hz to remove slow drifts in the data, and are resampled to 150 Hz to limit the atom size in the CDL. CDL is computed using alphacsc (Dupré la Tour et al., 2018) with the GreedyCDL method. Twenty atoms of duration 0.5 s are extracted from the signals. The extracted atoms'activations are binarized and, similarly as previous experiments, the events time are shifted to make them correspond to the peak amplitude time in the atom. Then, for every atom, the intensity function is estimated using the EM-based algorithm with 200 iterations and the "smart start" initialization strategy. Kernels' truncation values are set at [0 s ; 0.9 s]. Two drivers per atom are considered. The first driver contains timestamps of all bimodal trials (all frequencies combined) with in addition the 4 auditory unimodal stimuli (denoted catch0).

The second driver is similar to the first one, but instead of the auditory unimodal stimuli, it contains the 4 visual unimodal stimuli (denoted catch1).

This experiment is performed for 3 subjects, for each one we plot the 5 atoms that have the highest ratio α/µ, so that we automatically exhibit atoms that are highly linked to the presented stimuli.

Similarly as results presented in Section 4, we plot the spatial and temporal representation of each atom, as well as the two learned intensity functions that we plot "at kernel". Results are presented on First, we can observe that, as for experiments on sample and somato datasets, we recover the heartbeat artifact (atom 0 in Figure A.12) as well as the eye-blink artifact (atoms 0 in Figure A.11 and in Figure A.13). Because of the paradigm of the experiment, and in particular the instructions given to the subjects -that is, to press a button after every stimulus, the majority of them being visual -, it is not surprising if the eye blink artifact is slightly linked to the stimuli. Indeed, it is often observed in such experiments -where there is no designated time to blink -that the subject "allows" themself to blink after the visual stimulus.

As the majority of the presented stimuli are a combination of visual and auditory, the CDL model struggles to separate the two corresponding neural responses. Hence, that is why it can be observed that most of the first atoms reported exhibit a mixture between auditory and visual response in their topographies (first row). However, one can observe that for some such atoms, DriPP learned intensity is able to indicate what is the main contributing stimulus in the apparition of the atom. For instance, for atom 10 in

  Maximization step Once this assignation has been computed, one needs to update the parameters of the model using MLE. To obtain the update equations, we fix the probabilities P cancel the negative log-likelihood derivatives with respect to each parameter. For given values of probabilities P (n) k (t) and P

  α = 0.8, m = 0.4, σ = 0.05

Figure 3 :

 3 Figure3: Mean of the relative infinite norm as a function of process duration T and the percentage of events kept P/S, for two kernel shapes on synthetic data: wide kernel (left) and sharp kernel (right). The accuracy of the EM estimates increases with longer and denser processes.

Figure 4 :

 4 Figure4: Spatial and temporal patterns of 4 atoms from sample dataset, and their respective estimated intensity functions following a stimulus (cue at time = 0 s), for auditory and visual stimuli. The heartbeat and eye-blink artifacts are not linked to any stimuli. An auditory stimulus will induce a neural response similar to atom 2, with a mean latency of 85 ms.

Figure A. 1 :

 1 Figure A.1: Standard deviation of the relative infinite norm as a function of the process duration T and the percentage of events kept P/S, for two kernel shapes on synthetic data: wide kernel (left) and sharp kernel (right). The variance of the EM estimates continually decreases with longer and denser processes.

Figure A. 2 Figure A. 2 :

 22 Figure A.2 reports some details on the computation time of the experiment done on synthetic data which produced the Figure 3 and Figure A.1. For each value of T , the mean computation time is computed over 90 experiments (3 values of P/S times 30 random seeds). Results are presented in Figure A.2.

  Still on MNE sample dataset, we compute from raw data the epoch average following an auditory stimulus (stimuli in the left ear and in the right ear are combined) and plot on Figure A.7 the obtained

Figure A. 4 :

 4 Figure A.4: Influence of the kernel truncation upper bound b (with a = 0) on the obtained results on MNE sample dataset, for the 4 main atoms, for auditory and visual stimulus. By setting b too high, all intensity functions are completely flat, indicating a total loss of information.

Figure A. 5 :Figure

 5 Figure A.5: The first two components alongside with electrooculography (EOG) signal (top) and their corresponding topomaps (bottom), on MNE sample dataset. ICA000 can be associated to the heartbeat artifact, whereas ICA003 can be associated to the eye-blink artifact.

  Figure A.10 shows 3 atoms that correspond all to a µ-wave located in the secondary somatosensory region (S2), with three different shapes of kernels in their estimated intensity functions.A.7 EXPERIMENTS ON CAM-CAN DATASET

  Figure A.10: Spatial and temporal patterns of 3 µ-wave atoms from somato dataset, alongside with their respective estimated intensity functions. Below, in order to provide further information, are the corresponding brain locations for each atom obtained with dipole fitting.

FigureFigureFigure

  Figure A.11: Spatial and temporal patterns of the 5 atoms from Cam-CAN dataset subject CC620264, a 76.33-year-old female, and their respective estimated intensity functions following an audiovisual stimulus (cue at time = 0 s). Atoms are ordered by their bigger ratio α/µ.

Table 1 :

 1 On MNE sample dataset, statistical univariate Student t-test, H 0 : independence between an atom and the stimulus (auditory or visual). The atom ids correspond to the ones presented in Figure4.We have performed a statistical test using a student t-test to check if the mean activation probability on [a, b] segments is the same as the mean estimated on the baseline[0, T ] \ ti [t i + a, t i + b].We present the results of the performed test in Table1. We can clearly see that both artifacts (heartbeat and eye-blink) are not linked to stimuli while neural responses are linked to related stimuli (p-values are very small). We would like to stress that while this test is interesting, it is much more dependent on the selection of the support interval[a, b] than the proposed method. Moreover, this test does not allow to assess the latency of the responses, and whether the atom is an induced response or an evoked one.

	atom id atom type p-value auditory p-value visual
	0 1 2 6	heartbeat eye-blink 5.85 × 10 -1 2.25 × 10 -1 auditory 2.31 × 10 -97 visual 7.78 × 10 -1	1.19 × 10 -1 8.48 × 10 -1 6.31 × 10 -1 5.12 × 10 -50

-A.6 EXPERIMENTS ON MNE SOMATO DATASET

Both available at https://mne.tools/stable/overview/datasets_index.html

Available at https://www.cam-can.org/index.php?content=dataset

On somato dataset, in order to exhibit the induced response to the somatosensory stimulus on the left median nerve of the subject, we perform a time/frequency analysis, presented on Figure A.9. In order to perform this analysis, a complex process including the use of Morlet wavelets is performed.


Olivier Bernardi

AND Bertrand Duplantier
email: bertrand.duplantier@cea.fr

Philippe Nadeau
email: philippe.nadeau@univie.ac.at

A BIJECTION BETWEEN WELL-LABELLED POSITIVE PATHS AND MATCHINGS

A well-labelled positive path of size n is a pair (p, σ) made of a word p = p 1 p 2 . . . p n-1 on the alphabet {-1, 0, +1} such that P j i=1 p i ≥ 0 for all j = 1 . . . n-1, together with a permutation σ = σ 1 σ 2 . . . σn of {1, . . . , n} such that p i = -1 implies σ i < σ i+1 , while p i = 1 implies σ i > σ i+1 . We establish a bijection between well-labelled positive paths of size n and matchings (i.e. fixed-point free involutions) on {1, . . . , 2n}. This proves that the number of well-labelled positive paths is (2n -1)!! ≡ (2n -1) • (2n -3) • • • 3 • 1. By specialising our bijection, we also prove that the number of permutations of size n such that each prefix has no more ascents than descents is [(n -1)!!] 2 if n is even and n!! (n -2)!! otherwise.

It is shown in [1] that well-labelled positive paths of size n are in bijection with a collection of n-dimensional subpolytopes partitioning the polytope Πn consisting of all points (x 1 , . . . , xn) ∈ [-1, 1] n such that P j i=1 x i ≥ 0 for all j = 1 . . . n. Given that the volume of each subpolytope is 1/n!, our results prove combinatorially that the volume of Πn is

.

Introduction

A well-labelled path of size n is a pair (p, σ) made of a word p = p 1 p 2 . . . p n-1 on the alphabet {-1, 0, +1}, together with a permutation σ = σ 1 σ 2 . . . σ n of [n] ≡ {1, . . . , n} such that p i = -1 implies σ i < σ i+1 , while p i = 1 implies σ i > σ i+1 . We shall represent a path (p, σ) as a lattice path on Z 2 starting at (0,0) and made of steps (1, p i ) for i = 1 . . . n-1 together with the label σ i on the ith lattice point of the path for i = 1 . . . n. For instance, two well-labelled paths of size 10 are represented in Figure 1. A well-labelled path (p, σ) is said Motzkin (resp. In this paper, we define a bijection between well-labelled Motzkin paths of size n + 1 and matchings (fixed-point free involutions) on [2n]. We then define a closely related bijection between well-labelled positive paths of size n, and matchings on [2n]. This proves that these sets of paths are counted by (2n -1)!! ≡ (2n -1) • (2n -3) • • • 3 • 1. Both bijections also allow for a refined enumeration of well-labelled paths (p, σ) according to the number of horizontal steps (the number of letters 0 in p). More precisely, we show that the number of well-labelled positive paths of size n with k horizontal steps is k!(nk -1)!!(nk -3)!! if nk is even and 0 otherwise. Observe that well-labelled paths (p, σ) without horizontal steps are completely determined by the permutation σ. Indeed, in this case the word p encodes the up-down sequence of the permutation σ. Hence, by specialising our results to paths with no horizontal steps (i.e. k = 0), we enumerate permutations whose up-down sequence belong to a certain family. For instance, we prove that the number of permutations of size n such that each prefix has no more ascents than descents is [(n -1)!!] 2 if n is even and n!! (n -2)!! otherwise. We also prove that the number of permutations of size n having one more ascent than descent but such that each prefix has no more ascents than descents is (n -1)!! (n -3)!! if n is even and 0 otherwise. These enumerative results contrast with those in [START_REF] Carlitz | Permutations with prescribed pattern[END_REF][START_REF] De Bruijn | Permutations with given ups and downs[END_REF][START_REF] Foulkes | Enumeration of permutations with prescribed up-down and inversion sequences[END_REF][START_REF] Niven | A combinatorial problem of finite sequences[END_REF] by the fact that we consider here a family of admissible up-down sequences rather than a single sequence.

n k n-1 k k! [(n -k -1)!!]
Well-labelled positive paths appeared recently in a problem concerning the evaluation of the volume of the n-dimensional polytope Π n made of the points (x 1 , . . . , x n ) in [-1, 1] n such that j i=1 x i ≥ 0 for all j = 1 . . . n. Indeed, it was shown in [START_REF] Bernardi | The partition function of a freely-jointed chain in a half-space[END_REF] that the set of well-labelled positive paths of size n is in bijection with a set of n-dimensional subpolytopes forming a partition of Π n and this was our original motivation for studying well-labelled paths. Given that the volume of each subpolytope is 1/n!, our results prove combinatorially that the volume of Π n is (2n-1)!! n! . The paper is organised as follows. In Section 2, we define a recursive decomposition of well-labelled positive and Motzkin paths. We then translate these decompositions in terms of generating functions. For Motzkin paths, solving the generating function equation shows that the number of well-labelled Motzkin paths of size n + 1 is (2n -1)!!. From this, a simple induction shows that the number of well-labelled positive paths of size n is also (2n -1)!!. In Section 3, we use the recursive decomposition of paths in order to define bijections between well-labelled positive paths, well-labelled Motzkin paths and matchings. One step of these bijections uses a construction of Chen [START_REF] Chen | A general bijective algorithm for trees[END_REF] between labelled binary trees and matchings. Lastly in Section 4, we use our bijections to count well-labelled positive and Motzkin paths according to their number of horizontal steps. Specialising this results to the paths with no horizontal steps, we enumerate permutations whose up-down sequence belong to certain families mentioned above.

Decomposition of well-labelled paths

In this section, we define a recursive decomposition of the class A of well-labelled Motzkin paths and the class B of well-labelled positive paths. We then translate these equations in terms of generating functions and obtain our first counting results.

We denote respectively by A n and B n the sets of paths of size n in A and B. We denote respectively by a n and b n the cardinality of A n and B n and by

A(z) = n≥0 a n n! z n and B(z) = n≥0 b n n! z n
the corresponding exponential generating functions. Observe that a 0 = a 1 = 0 and b 0 = 0. The following notation will be useful for relabelling objects: given a set I of n integers, we denote by λ I the order preserving bijection from [n] to I (and by λ -1 I the inverse bijection).

2.1. Decomposition of well-labelled Motzkin paths. We first define a recursive decomposition of the class A of well-labelled Motzkin paths. For i ∈ {-1, 0, 1}, we denote by A i (resp. A i n ) the set of paths (p, σ) in A (resp. A n ) such that p 1 = i. Observe that A -1 is made of a single element α 2 of size 2. The decomposition

A = {α 2 } ⊎ A 0 ⊎ A 1
is illustrated by Figure 2 and the following proposition reveals its recursive nature.

Proposition 1. For any positive integer n,

• the set A 0 n is in bijection with the set [n] × A n-1 , • the set A 1
n is in bijection with the set C n made of all unordered pairs {(I ′ , P ′ ), (I ′′ , P ′′ )} such that I ′ ⊆ [n], I ′′ = [n] \ I ′ and P ′ , P ′′ are well-labelled Motzkin paths of respective size |I ′ | and

|I ′′ |. = α 2 A 0 A 1 A ⊎ ⊎ Figure 2.
Recursive decomposition of well-labelled Motzkin paths.

Proof. • For any path (p, σ) in A 0 n one obtains a pair (k, (p

′ , σ ′ )) in [n]×A n-1 by setting k = σ 1 , p ′ = p 2 . . . p n-1 and σ ′ = σ ′ 1 . . . σ ′ n-1 where σ ′ i = λ -1 [n]\{k} (σ i+1 ) for i = 1 . . . n
= p k+i , • σ ′ = σ ′ 1 . . . σ ′ k and σ ′′ = σ ′′ 1 . . . σ ′′ n-k , where σ ′ i = λ -1 I ′ (σ k+1-i ) and σ ′′ i = λ -1 I ′′ (σ k+i )
. The mapping φ is clearly a bijection between the sets A 1 n and C ′ n , which concludes the proof. Corollary 2. The generating function A(z) of well-labelled Motzkin paths satisfies

A(z) = z 2 2 + zA(z) + A(z) 2 2 . (1)
Proof. For i ∈ {-1, 0, 1}, we denote by a i n the cardinality of A i n and by

A i (z) = n≥0 a i n n! z n the corresponding generating function. The partition A = {α 2 } ⊎ A 0 ⊎ A 1 gives A(z) = z 2 2 + A 0 (z) + A 1 (z).
Moreover, the bijection between

A 0 n and [n]×A n-1 gives a 0 n = n a n-1 , hence A 0 (z) = zA(z) while the correspondence between A 1 n and C n gives a 1 n = 1 2 n k=0 n k a k a n-k , hence A 1 (z) = A(z) 2 2 .
By solving Equation (1) (and using the fact that a 0 = 0), one gets

A(z) = 1 -z - √ 1 -2z. (2)
One can extract the coefficient a n either directly from this expression of A(z) or by applying Lagrange inversion formula to the series

C(z) = A(z)/z. Indeed, Equation (1) gives C(z) = z (1+C(z)) 2 2
, hence

a n+1 = (n + 1)![z n ]C(z) = (n + 1)! n [x n-1 ] (1 + x) 2 2 n = (2n)! 2 n n! = (2n -1)!!.
We will now determine the number b n of well-labelled positive paths of size n. This can be done by exploiting a bijection between B n × [n + 1] × {0, 1} and B n+1 ⊎ A n+1 obtained by adding one step to a positive path. The bijection is as follows: given a well-labelled positive path (p, σ) of size n, an integer k in [n + 1] and an integer b in {0, 1}, we define the labelled path (p

′ , σ ′ ) = ψ((p, σ), k, b) by setting • σ ′ = σ ′ 1 . . . σ ′ n+1 , where σ ′ i = λ [n+1]\{k} (σ i ) for i = 1 . . . n and σ ′ n+1 = k, • p ′ = p ′ 1 . . . p ′ n , where p ′ i = p i for i = 1 . . . n-1 and p ′ n is equal to b -1 if σ ′ n+1 > σ ′
n and equal to b otherwise. Observe that the path (p ′ , σ ′ ) is well-labelled (by the choice of the step p ′ n ) and is either positive or Motzkin (since (p, σ) is positive). Moreover, the mapping ψ is a bijection between

B n × [n + 1] × {0, 1} and B n+1 ⊎ A n+1 showing that 2(n + 1)b n = b n+1 + a n+1 for all n ≥ 0. ( 3
)
Since a n+1 = (2n -1)!! a simple induction shows that b n = (2n -1)!! and proves the following. Proposition 3. The number a n+1 of well-labelled Motzkin paths of size n + 1 and the number b n of well-labelled positive paths of size n are both equal to (2n -1)!!.

2.2. Decomposition of well-labelled positive paths. We now define a recursive decomposition of the class B of well-labelled positive paths. We denote by β 1 the well-labelled path of size 1 and for i ∈ {0, 1}, we denote by B i the set of paths (p, σ) in B of size at least 2 satisfying p 1 = i. For a path (p, σ) of size n in B 1 , we consider the greatest integer k ≤ n such that j-1 i=1 p i ≥ 1 for all j = 2 . . . k-1 and k-1 i=1 p i = 1. We denote by B ′ the subset of paths in B 1 such that k = n and we denote B ′′ = B 1 \ B ′ the complement. We also denote by B 0 n , B ′ n and B ′′ n respectively the paths of size n in B 0 , B ′ and B ′′ . The partition

B = {β 1 } ⊎ B 0 ⊎ B ′ ⊎ B ′′
is illustrated by Figure 3 and the following proposition reveals its recursive nature. Proof. • For any path (p, σ) in B 0 n one obtains a pair (k, (p

′ , σ ′ )) in [n] × B n-1 by setting k = σ 1 , p ′ = p 2 . . . p n-1 and σ ′ = σ ′ 1 . . . σ ′ n-1 where σ ′ i = λ -1 [n]\{k} (σ i+1 ) for i = 1 . . . n
I ′′ = {σ k+1 , . . . , σ n }, • p ′ = p ′ 1 . . . p ′ k-1 and p ′′ = p ′′ 1 . . . p ′′ n-k-1 , where p ′ i = -p k-i and p ′′ i = p k+i , • σ ′ = σ ′ 1 . . . σ ′ k and σ ′′ = σ ′′ 1 . . . σ ′′ n-k , where σ ′ i = λ -1 I ′ (σ k+1-i ) and σ ′′ i = λ -1 I ′′ (σ k+i
). We first want to prove that ((I ′ , (p ′ , σ ′ )), (I ′′ , (p ′′ , σ ′′ ))) is in D n . It is clear that (p ′ , σ ′ ) and (p ′′ , σ ′′ ) are well-labelled paths and moreover, (p ′ , σ ′ ) is a Motzkin path. It remains to prove that (p ′′ , σ ′′ ) is a positive path. Observe that the step p k is non-zero otherwise it contradicts the maximality of k. If p k = -1, then (p ′′ , σ ′′ ) is clearly positive because (p, σ) is positive; and if p k = +1, then (p ′′ , σ ′′ ) is positive otherwise it would contradict the maximality of k. Hence, φ is a mapping from A n to D n . The bijectivity of φ is easy to check after observing that the step p k can be recovered: it is equal to 1 if λ I ′ (σ ′ 1 ) < λ I ′′ (σ ′′ 1 ) and to -1 otherwise.

Proposition 4 will allow to define a bijection between positive paths and matchings in the next section. It also leads to the following relation between the generating functions A(z) and B(z):

(4) B(z) = z + zB(z) + A(z) + A(z)B(z),
which, by (2), gives

B(z) = 1 √ 1 -2z -1.
This result could also have been derived from the observation that b n = a n+1 implies B(z) = A ′ (z).

Bijections with matchings

Proposition 3 suggests that the classes of paths A n+1 and B n are both in bijection with matchings. The goal of this section is to describe such bijections. For this, we will introduce intermediate objects called labelled binary trees.

3.1.

Bijections between well-labelled paths and labelled binary trees. A labelled binary tree of size n is a rooted tree with n leaves having n different labels in [n] and such that each (unlabelled) internal vertex has exactly two unordered children. We call marked labelled binary tree a labelled binary tree in which one of the (internal or external) vertices is marked. A binary tree and a marked binary tree are represented in Figure 4. We denote by T the set of labelled binary trees of size at least 2 and we denote by R the set of marked labelled binary trees. We will now show that the recursive descriptions of the classes T and R parallel those of the classes A and B and obtain bijections between T and A and between R and B.

We use the following notation for relabelling trees: if λ is a bijection between two sets of integers I, J and τ is a binary tree whose leaves have labels in I, then λ(τ ) denotes the tree obtained from τ by replacing each leaf labelled i ∈ I, i = 1 . . . n by a leaf labelled λ(i) ∈ J. Bijection between Motzkin paths and labelled trees. We denote by τ 2 the unique labelled binary tree of size 2 and by T 0 (resp. T 1 ) the set of labelled binary trees made of the root-vertex and two subtrees, one of which (resp. none of which) is a leaf. For all integers n ≥ 2, we denote respectively by T n , T 0 n and T 1 n the set of trees of size n in T , T 0 and T 1 . We now explicit the recursive nature of the partition T = {τ 2 } ⊎ T 0 ⊎ T 1 .

Proposition 5. For any integer n > 2,

• the set T 0 n is in bijection with the set

[n] × T n-1 . • the set T 1
n is in bijection with the set U n of unordered pairs {(I ′ , τ ′ ), (I ′′ , τ ′′ )} such that

I ′ ⊆ [n], I ′′ = [n] \ I ′ and τ ′ , τ ′′ are labelled binary trees in T of respective size |I ′ | and |I ′′ |. Proof. Let n > 2.
• Given an integer k ∈ [n] and a tree τ in T n-1 , one defines a tree τ ′ in T 0 n as the tree made of a root-vertex, one subtree which is a leaf labelled k and one subtree equal to λ [n]\{k} (τ ) (which is not a leaf). This is a bijection between [n] × T n-1 and T 0 n . • Given a pair {(I ′ , τ ′ ), (I ′′ , τ ′′ )} in U n , one defines a tree τ in T 1 n as the tree made of a rootvertex, a subtree equal to λ I ′ (τ ′ ) and a subtree equal to λ I ′′ (τ ′′ ). This is a bijection between U n and T 1 n . Definition of bijection Φ. Comparing Propositions 1 and 5, it is clear that one can define a recursive bijection Φ between the class A of well-labelled Motzkin paths and the class T of labelled binary trees. We now summarise the recursive construction of the image of a welllabelled Motzkin path (p, σ) by the bijection Φ:

(i) If (p, σ) = α 2 , then Φ(p, σ) = τ 2 .
(ii) If (p, σ) has size n > 2 and p 1 = 0, then we set p ′ = p 2 . . . p n-1 and σ

′ = σ ′ 1 . . . σ ′ n-1
where σ ′ i = λ -1

[n]\{σ1} (σ i+1 ) for i = 1 . . . k-1. With these notations, we define Φ(p, σ) as the tree made of a root-vertex, the subtree made of a leaf labelled σ 1 and the subtree λ [n]\{σ1} (Φ(p ′ , σ ′ )). (iii) If (p, σ) has size n > 2 and p 1 = 1, then we consider the least integer k such that k i=1 p i = 0 and we set (as in the proof of Proposition 1):

• I ′ = {σ 1 , . . . , σ k } and I ′′ = {σ k+1 , . . . , σ n }, • p ′ = p ′ 1 . . . p ′ k-1 and p ′′ = p ′′ 1 . . . p ′′ n-k-1 , where p ′ i = -p k-i and p ′′ i = p k+i , • σ ′ = σ ′ 1 . . . σ ′ k and σ ′′ = σ ′′ 1 . . . σ ′′ n-k , where σ ′ i = λ -1 I ′ (σ k+1-i ) and σ ′′ i = λ -1 I ′′ (σ k+i )
. With these notations, we define Φ(p, σ) as the tree made of a root-vertex, the subtree λ I ′ (Φ(p ′ , σ ′ )) and the subtree λ I ′′ (Φ(p ′′ , σ ′′ )). For instance, the image of the Motzkin path represented in Figure 1(a) by the mapping Φ is represented in Figure 4(a). From the definition of Φ and Propositions 1 and 5, we have the following theorem: Theorem 6. For any positive integer n > 1, the mapping Φ is a bijection between well-labelled Motzkin paths of size n and labelled binary trees with n leaves.

Bijection between positive paths and marked trees. We now define a bijection Φ ′ between well-labelled positive paths and marked labelled binary trees. Before defining the bijection Φ ′ , let us explain briefly what led us to consider marked labelled binary trees. As seen in Section 2, the recursive decomposition of positive paths leads to consider blocks corresponding to either positive paths or Motzkin paths. The recursive relation is captured by Equation (4) which can be written B(z) = Ã(z) + Ã(z)B(z), where Ã(z) = z + A(z) is the series of unmarked binary trees (of size n ≥ 1), or equivalently, binary trees marked at their root-vertex. This relation suggests that one can interpret B(z) as counting marked binary trees.

We denote by ρ 1 the marked labelled tree of size 1. We denote by R ′ ⊂ R the set of marked trees of size at least 2 such that the marked vertex is the root. Clearly, this set is in bijection with the set T of unmarked trees. We denote by R 0 (resp. R ′′ ) the set of marked trees made of a non-marked root and two subtrees, one of which (resp. none of which) is a non-marked leaf. For all integer n > 1, we denote respectively by R n , R 0 n , R ′ n and R ′′ n the set of marked trees of size n in R, R 0 , R ′ and R ′′ . We now explicit the recursive nature of the partition

R = {ρ 1 } ⊎ R 0 ⊎ R ′ ⊎ R ′′ .
Proposition 7. For all integer n > 1,

• the set R 0 n is in bijection with the set

[n] × R n-1 , • the set R ′ n is in bijection with T n , • the set R ′′
n is in bijection with the set U n of ordered pairs ((I ′ , τ ′ ), (I ′′ , τ ′′ )) such that

I ′ ⊆ [n], I ′′ = [n] \ I ′ , τ ′ is
a non-marked labelled binary tree of size |I ′ | and τ ′′ is a marked labelled binary tree of size |I ′′ |.

The proof of Proposition 7 is similar to the proof of Proposition 5 and is omitted.

Definition of bijection Φ ′ . Comparing Propositions 4 and 7, it is clear that one can define a recursive bijection Φ ′ between the class B of well-labelled paths and the class R of marked labelled binary trees. We now summarise the recursive construction of the image of a welllabelled positive path (p, σ) by the bijection Φ ′ :

(i) If (p, σ) = β 1 then Φ ′ (p, σ) is the marked tree ρ 1 . (ii) If (p, σ) has size n > 1 and p 1 = 0, then we define p ′ = p 2 . . . p n-1 and σ ′ = σ ′ 1 . . . σ ′ n-1 , where σ ′ i = λ -1 [n]\{σ1} (σ i+1 ) for i = 1 . . . n-1.
With these notations, we define Φ ′ (p, σ) as the tree made of a non-marked root-vertex, the subtree made of a non-marked leaf labelled σ 1 and the marked subtree λ [n]\{σ1} (Φ ′ (p ′ , σ ′ )). (iii) If (p, σ) has size n > 1 and p 1 = 1, then we consider the greatest integer k ≤ n such that j-1 i=1 p i ≥ 1 for all j = 2 . . . k-1 and k-1 i=1 p i = 1, and we set (as in the proof of Proposition 4):

• I ′ = {σ 1 , . . . , σ k } and I ′′ = {σ k+1 , . . . , σ n }, • p ′ = p ′ 1 . . . p ′ k-1 and p ′′ = p ′′ 1 . . . p ′′ n-k-1 , where p ′ i = -p k-i and p ′′ i = p k+i , • σ ′ = σ ′ 1 . . . σ ′ k and σ ′′ = σ ′′ 1 . . . σ ′′ n-k , where σ ′ i = λ -1 I ′ (σ k+1-i ) and σ ′′ i = λ -1 I ′′ (σ k+i ). If k = n (that is, (p, σ) is in B ′ ),
we define Φ ′ (p, σ) as the marked tree obtained by marking the root-vertex of the unmarked tree Φ(p ′ , σ ′ ) (note that I ′′ , p ′′ and σ ′′ are empty in this case). Otherwise (that is, if k < n), we define Φ ′ (p, σ) as the marked tree made of a non-marked root-vertex, the non-marked subtree λ I ′ (Φ(p ′ , σ ′ )) and the marked subtree λ I ′′ (Φ ′ (p ′′ , σ ′′ )). For instance, the image of the positive path represented in Figure 1(b) by the mapping Φ ′ is represented in Figure 4(b). From the definition of Φ ′ and Propositions 4 and 7, we have the following theorem: Theorem 8. The function Φ ′ is a bijection between well-labelled positive paths of size n and marked labelled binary trees with n leaves.

Remark. The final height of a positive path (p, σ) of size n is n-1 i=1 p i . It is not hard to prove inductively that a positive path (p, σ) has an even final height if and only if the mark of the image tree, ρ = Φ ′ (p, σ), is on a leaf. Indeed, if (p, σ) = β 1 the final height is 0, and the mark is on a leaf of the tree ρ = ρ 1 ; while if (p, σ) is in B ′ the final height is 1, and the mark is on an internal vertex (the root-vertex) of ρ. In the other cases ((p, σ) ∈ B 0 ⊎ B ′′ ), the mark is in a subtree ρ ′ of ρ corresponding to a path (p ′ , σ ′ ) having a final height of the same parity as (p, σ).

3.2.

Bijections between labelled binary trees and matchings. We will now present a bijection Ψ due to Chen [START_REF] Chen | A general bijective algorithm for trees[END_REF] between labelled binary trees of size n and matchings on [2n-2]. We follow the exposition from [7, p.16] for defining the bijection Ψ and then define a similar bijection Ψ ′ between marked binary trees of size n and matchings on [2n]. The mappings Ψ and Ψ ′ are represented in Figure 5. The first step of this bijection is to attribute a label to each internal node of the binary tree. 0 0 1 1 00 00 11 11 00 00 11 11 0 0 1 1 0 0 1 1 00 00 11 11 00 00 11 11 00 00 11 11 0 0 1 1 0 0 1 1 00 00 11 11 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 00 00 11 11 0 0 1 1 0 0 1 Definition of bijection Ψ [START_REF] Chen | A general bijective algorithm for trees[END_REF]. Let τ be a labelled binary tree with n leaves labelled on [n]. One defines an induced labelling of the n-2 internal non-root vertices of τ by the following procedure. While there are unlabelled non-root vertices, we consider those among them that have both of their children labelled. There is at least one such vertex and we consider the vertex v which has the child with least label; we then attribute to v the least label in [2n-2] \ [n] which has not yet been attributed. For instance, starting from the tree in Figure 4(a) one obtains the labels indicated in the fully labelled tree represented in Figure 5(a). One then obtains the matching π = Ψ(τ ) on [2n -2] by letting π(i) = j for all pairs i, j ∈ [2n -2] which are the labels of siblings in the fully labelled tree. In [7, p.16] the mapping Ψ is proved to be a bijection between labelled binary trees of size n and matchings on [2n-2].

Definition of bijection Ψ ′ . Let τ ′ be a marked labelled binary tree of size n and let τ be the corresponding unmarked tree. We consider the matching π = Ψ(τ ) on [2n-2] and define a matching π ′ = Ψ ′ (τ ′ ) on [2n] as follows. If the marked vertex v of τ ′ is the root, then π ′ (i) = π(i) for all i in [2n-2] and π ′ (2n -1) = 2n. If v is not the root, then we consider its label k and the label l of its sibling. In this case π ′ (i) = π(i) for all i = k, l in [2n -2], π ′ (k) = 2n and π ′ (l) = 2n -1. It is clear that the mapping Ψ ′ is a bijection between marked labelled binary tree of size n and matchings on [2n].

Combining the bijections Φ, respectively Φ ′ , with the bijections Ψ, respectively Ψ ′ , gives the following bijective proof of Proposition 3.

Enumerative corollaries

We will now study the number of horizontal steps in well-labelled paths through the bijections Φ, Φ ′ , Ψ, Ψ ′ and deduce some enumerative corollaries in terms of the up-down sequences of permutations. Recall that a horizontal step of a well-labelled path (p, σ) is a letter 0 in p. We say that a leaf in a labelled binary tree is single if its sibling is an internal node.

Theorem 10. For all integers n, k, the mappings Φ and Ψ induce successive bijections between • well-labelled Motzkin paths of size n with k horizontal steps,

• labelled binary trees with n leaves, k of which are single leaves,

• matchings on [2n -2] having k pairs (i, j) such that i ∈ {1, . . . , n} and j ∈ {n+1, . . . , 2n-2}.

For example, the Motzkin path of size n = 10 in Figure 1(a) has 2 horizontal steps, the corresponding labelled binary tree represented in Figure 4(a) has 2 single leaves, and the corresponding matching represented in Figure 5(a) has 2 pairs (i, j) such that i ∈ {1, . . . , n} and j ∈ {n+1, . . . , 2n-2}.

Proof. • The correspondence between the number of horizontal steps of a Motzkin path (p, σ) and the number of single leaves in the binary tree Φ(p, σ) follows from a simple induction on the size of (p, σ). Indeed, one creates a single leaf in the recursive construction of Φ(p, σ) exactly when case (ii) thereof (corresponding to a horizontal step of (p, σ)) occurs.

• The correspondence between the number of single leaves in the binary tree τ and the number of pairs (i, j) such that i ∈ {1, . . . , n} and j ∈ {n + 1, . . . , 2n -2} in the matching Ψ(τ ) is an immediate consequence of the fact that the labels of external vertices are in {1, . . . , n} while the labels of internal vertices are in {n+1, . . . , 2n-2}.

Corollary 11. The number of well-labelled Motzkin paths of size n having k horizontal steps is

a n,k = n k n -2 k k! (n -k -1)!! (n -k -3)!!
if nk is even and 0 otherwise.

Proof. By Theorem 10, the number a n,k counts matchings on [2n -2] with exactly k pairs (i, j) such that i ∈ {1, . . . , n} and j ∈ {n+1, . . . , 2n-2}. To enumerate such matchings, first choose these k pairs: there are n k possibilities of choosing the integers i in {1, . . . , n}, there are n-2 k possibilities for choosing the integers j in {n+1, . . . , 2n-2} and then k! possibilities to define the pairing between the chosen integers in {1, . . . , n} and the chosen integers in {n + 1, . . . , 2n -2}.

After that, it remains to choose a pairing of the nk unmatched integers in {1, . . . , n} and a pairing of the nk -2 unmatched integers in {n+1, . . . , 2n-2}. Such matchings exist only if nk is even and there are (nk -1)!!(nk -3)!! of them in this case.

We now examine horizontal steps in positive paths. We say that a leaf in a marked labelled binary tree is quasi-single if it is not marked and its sibling is either marked or internal.

Theorem 12. For all integers n, k, the mappings Φ ′ and Ψ ′ induce successive bijections between • well-labelled positive paths of size n with k horizontal steps,

• marked labelled binary trees of with n leaves, k of which are quasi-single leaves,

• matchings on [2n] having k pairs (i, j) with i ∈ {1, . . . , n} and j ∈ {n + 1, . . . , 2n -1}.

  0 for all j = 1 . . . n-2 and n-1 i=1 p i = -1 (resp.

p

  i ≥ 0 for all j = 1 . . . n-1).

Figure 1 .

 1 Figure 1. (a) A well-labelled Motzkin path. (b) A well-labelled positive path.

  nk)!! (nk -2)!! otherwise. Similarly, the number of well-labelled Motzkin paths of size n with k horizontal steps is n k n-2 k

Proposition 4 .Figure 3 .

 43 Figure 3. Recursive decomposition of well-labelled positive paths.

Figure 4 .

 4 Figure 4. (a) A labelled binary tree. (b) A marked labelled binary tree

Figure 5 .

 5 Figure 5. Bijections between labelled binary trees and matchings.

Theorem 9 .

 9 The composition Ψ • Φ is a bijection between well-labelled Motzkin paths of size n and matchings on [2n-2]. Similarly, Ψ ′ • Φ ′ is a bijection between well-labelled positive paths of size n and matchings on [2n].

  -1. This is clearly a bijection between A 0 n and [n] × A n-1 . • Observe that the set C n is trivially in bijection with the set C ′ n of ordered pairs ((I ′ , P ′ ), (I ′′ , P ′′ )) such that λ I ′ (σ ′ 1 ) < λ I ′′ (σ ′′ 1 ). Let (p, σ) be a path in A 1 n and let k be the least integer such that

	k i=1 p i = 0 (observe that p k = -1). We define an element φ(p, σ) = ((I ′ , (p ′ , σ ′ )), (I ′′ , (p ′′ , σ ′′ ))) n by setting of C ′ • I ′ = {σ 1 , . . . , σ k } and I ′′ = {σ k+1 , . . . , σ n }, • p ′ = p ′ 1 . . . p ′ k-1 and p ′′ = p ′′ 1 . . . p ′′ n-k-1 , where p ′ i = -p k-i and p ′′ i

Acknowledgement: We are very thankful to Sylvie Corteel for fruitful discussions and for providing us with the conjectural formula for the numbers b n,k .

Proof. • The correspondence between the number of horizontal steps of a positive path (p, σ) and the number of quasi-single leaves in the marked tree Φ ′ (p, σ) follows from a simple induction on the size of (p, σ). Indeed, one creates a quasi-single leaf in the recursive construction of Φ ′ (p, σ) exactly when case (ii) of the definition of either Φ or Φ ′ occurs.

• We now consider a marked labelled binary tree τ ′ of size n > 1. Let v be the marked vertex and let τ be the non-marked tree obtained by forgetting the mark. The number k of single leaves in τ and the number k ′ of quasi-single leaves in τ ′ are related by

and its sibling are both leaves. Similarly, the definition of Ψ ′ gives a relation between the number l ′ of pairs (i, j) of the matching Ψ ′ (τ ′ ) such that i ∈ {1, . . . , n} and j ∈ {n + 1, . . . , 2n -1} and the number l of pairs (i, j) of the matching Ψ(τ ) such that i ∈ {1, . . . , n} and j ∈ {n + 1, . . . , 2n -2}:

• l ′ = l if the label of v is larger than n,

• l ′ = l -1 if the label of v is not larger than n and the label of its sibling is larger than n,

• l ′ = l + 1 if the label of v and its sibling are not larger than n. Theorem 10 gives k = l, hence the previous discussion gives k ′ = l ′ and concludes the proof.

Corollary 13. The number of well-labelled positive paths of size n having k horizontal steps is

The proof of Corollary 13 is very similar to the Corollary 11 and is omitted. We now study the consequence of these results in terms of the up-down sequences of permutations.

An ascent of a permutation σ = σ 1 σ 2 . . . σ n is an index i < n such that σ i < σ i+1 ; a descent is an index i < n such that σ i > σ i+1 . The enumeration of permutations with a given sequence of ascents and descents, called up-down sequences (or shape) was investigated for instance in [START_REF] Carlitz | Permutations with prescribed pattern[END_REF][START_REF] De Bruijn | Permutations with given ups and downs[END_REF][START_REF] Foulkes | Enumeration of permutations with prescribed up-down and inversion sequences[END_REF][START_REF] Niven | A combinatorial problem of finite sequences[END_REF]. Here we will count permutations of size n such that their up-down sequences belong to a certain family, while previous works focused on the enumeration according to a single updown sequence.

We say that a permutation σ has a positive up-down sequence if for all j ≤ n the number of ascents less than j is no more than the number of descents less than j. We say that σ has a Dyck up-down sequence if it has one more ascent than descents but for all j < n the number of ascents less than j is no more than the number of descents less than j. Observe that a well-labelled path (p, σ) with no horizontal steps is completely determined by the permutation σ (indeed, the word p is determined by the up-down sequence of σ). Moreover, the well-labelled path (p, σ) is positive (resp. Motzkin) if and only if the permutation σ has a positive (resp. Dyck) up-down sequence. Thus, the following theorem immediately follows by looking at the specialisation k = 0 in Corollaries 11 and 13. Theorem 14. For any integer n, the number of permutations of size n having a positive updown sequence is [(n -1)!!] 2 if n is even and n!! (n -2)!! otherwise. The number of permutations of size n having a Dyck up-down sequence is (n -1)!! (n -3)!! if n is even and 0 otherwise.

We would be happy to see a more direct bijective proof of these specialisations.