Associate Professor, S F Masri

Steady-

Numbers in brackets designate References at end of paper. Contributed by the Applied Mechanics Division and presented at the Applied Mechanics Summer Conference, University of California, La Jolla, Calif., June 26-28, 1972, of THE AMERICAN SOCIETY OP MECHANICAL ENGINEERS.

Discussion on this paper should be addressed to the Editorial Department, ASME, United Engineering Center, 345 East 47th Street, New York, N. Y. 10017, and will be accepted until April 20, 1973. Discussion received after this date will be returned. Manuscript received by ASME Applied Mechanics Division, September, 1971; final revision, January, 1972. Paper No. 72-APM-39. motion is consistent with that which has been found to predominate in the experimental studies of impact dampers as observed by most investigators in this field.

Steady-State Solution

The equation of motion of the system in the absence of the damper is

[m]x + [c]x + Wx = F(0 (1)
where [m], [c], [k] are the mass, damping, and stiffness matrices, respectively, and F(<) = col {0, 0, . . . , 0, F h (t), 0, . . ., 0}. Assume that the damping mechanism is the proportional-type satisfying the condition

[c] = a[m] + 0[k] (2)
where a and /3 are constants.

Consider the steady-state motion of the system shown in Pig. 1 with the origin of the time axis shifted to coincide with the time of occurrence of an impact at t = to. The result of this shift is to modify the excitation force, giving in the new time scale

F h (t) = F 0 sin (Qt + ao) (3)
where a 0 = ilto is a phase angle to be determined from the steadystate motion.

Using the normal mode approach, equation ( 1) can be transformed into the form

l^M<\q + [ x C\]q + [^\]q = Q«M (4)
where the n-components of q are the normal coordinates, [ x iW\], [ X C\], and [ X .K\] are diagonal matrices corresponding to the

M&i + C4i + K,Qi = Qex,» = <PkiFo sin (tit + a 0 ) (5)
and its solution is

?,« = exp(- - tit)\- UiSiny-Slt + 77; cos -' -tit ) qt>i H I sin --tit ) q ai Ai ( . Vi o, . i)i n \ • ~ -I f • sln -tit + Vi cos -tit I sin Ti Vi \ r i '•» / --r v I sin --tit) cos T t -} + A< sin (tit + 1i V r ( I J i=l,2,...,n (6)
where

-Vf; Vi = Vi -$-,-2 fi = VuFo 2£V,- ^ == tan" 1 -~L ± - 1 -r,-» ?oi = ?i(0) f« = C, ti A t = fi/K< V(i (1 -rW + (2f<r 4 )» g«i = g,(0+)
and the + subscript indicates conditions immediately after the specified time. Letting the initial displacement and velocity at t = 0 + be

x(0) = xo = [<p]qo, x(0 + ) = x" = [<p]q a (7) then x(0 = [Bn(.t))k a + [B, s (*)]xo + [B 23 «]Si + [S 24 (01S 2 + MS"(0 (8) k(t) = [B 31 (*)]X. + [B3 2 («)]xo+ [B"a(0]Si + [iJ 3 4(*)]S2 + [^]S 4 («) (9)
where the undefined matrices and vectors are functions of the system parameters. Letting y{t) be the relative displacement of the damper mass m d with respect to m jt then

y{t) = z(t) -XjOL). ( 10
)
Since the time origin coincides with the time of occurrence of an impact,

2/(0)^=2/0 = z(0) -xj(0) = ±- (11)
During an impact (idealized to be a discontinuous process), the conditions of the system remain unchanged except for the velocities of m d and nij, whose instantaneous changes are calculated using the momentum equation and the coefficient of restitution. Noting that in steady-state motion 2ti , z(0)+ = (x oy + ?/o) 7T

then the velocity vectors before and after an impact are related by

it b = [ s .Bes]x 0 (13)
-Nomenclature" rrit = ith mass in multidegree-of-freedom system hi -ith spring constant in multidegree-of-freedom system a, fi = constants used to specify proportional-type damping where [^.BON] is a constant diagonal matrix whose elements are equal to unity, except for the jth element which is equal to (1 -e -2/*)/(l -e -2ne).

Xj(t) = absolute displacement
In steady-state motion with two symmetric impacts/cycle of excitation on opposite ends of damper container,

x(«)|«=»-= -x(0) = -xo (14) x(«)|w=^_ = -x(0)_ = -x, = -[\B 6 s]x" (15)
Using ( 14) and (15) in conjunction with (8) and (9) xo = S7 sin a 0 + Sg cos a 0 (16)

k a = S9 sin a 0 + Sw cos ao (17)
Equations ( 16) and ( 17), together with (13), result in

a 0 ± = tan" 1 [(hh, ± kJi^/ihzh, =F hju)] ( 18
)
where the h's are functions of S 7 -Sio. With a 0 determined from (18), the rest of the unknowns can be found by back substitution.

Applications

Consider the 10-degree-of-freedom lumped parameter system shown in Fig. 2(c) whose natural frequencies, mode shapes, and damping characteristics approximate those obtained experimentally from full-scale dynamic tests of a modern 10-story building [START_REF] Kuroiwa | Vibration Test of a Multistory Building[END_REF].

To investigate the effects of various system parameters, the models shown in Figs. 2(a) and (6) were also studied. The hypothetical systems in Figs. 2(a) and (6) consist of the first 4 and first 7 stories, respectively, of the system in Fig. 2(c). Fig. 2(d) shows the first three mode shapes and the corresponding frequencies for the lumped parameter representation of the 10story building in Fig. 2(c). Fig. 3 shows a typical solution curve for the 7-story building in Fig. 2(6). The left-hand-side (LHS) ordinate in Fig. 3 is the ratio of the maximum displacement of the ith floor, for the structure using the damper divided by the peak displacement, x vi , of the same floor with no damper being used. This ratio is, therefore, a measure of the effectiveness of the damper (without a damper the solution curve is a straight line at ima*,/ x v% -!)• The abscissa in Fig. 3 is the clearance ratio d/x s t", where x"t" is the static deflection of the top floor with the excitation force and the damper both applied at the top.

For a given set of parameters there are two possible solutions to equation (18) labeled a 0 + and a 0 ~, which correspond to two distinct steady-state solutions. Experiments with an electronic analog computer were conducted to check the theoretical results and to determine which branch of the solution, if any, was asymptotically stable. As indicated by the experimental results shown in Fig. 3, the solution curve associated with a 0 + (which, in general, leads to lower response amplitude than solution curve a 0 ~) was found to be stable throughout most of its range.

The validity of the analytical results were further compared to, and were found to be in excellent agreement with, published analytical and experimental results pertaining to single-degreeof-freedom systems provided with impact dampers [START_REF] Masri | On the Stability of the Impact Damper[END_REF].

Journal of Applied Mechanics

Effects of System Parameters. A detailed investigation of the independent parameters on the response of the system was undertaken, and the results can be summarized as follows:

1 Mass Ratio Mr. The system under discussion is effective even with mass ratios on the order of 1 percent. Note from Fig. 3 that even with a mass ratio n T = 0.005 the response amplitude can be reduced by as much as 80 percent.

For "tall" structures operating in the vicinity of their first mode, similar to the one under discussion, the optimum clearance ratio (in regard to vibration attenuation) satisfies the condition Hjd/xttn » 0.6 and results in a reduced amplitude which satisfies the relationship /x T ^mai,/a;p, » 0.001, (i = 1, 2, . . . , re) for relatively high values of e.

2 Natural Mode Shapes. The response of the 10-story building under discussion vibrating in each of its first three modes is given in Pig. 4. It can be seen that, with a given damper mass, the maximum percentage reduction in the response of the structure is achieved in the first mode. The maximum reduction decreases with mode number from 90 percent in the first, to 85 percent in the second, and 60 percent in the third. One of the reasons for this reduction in the efficiency of the damper is that f " the ratio of critical damping in the ith mode, increases with mode number from ft » 0.01 to f 2 » 0.02 and f 8 « 0.03. Figs. 4(a), (6), and (c) also show the corresponding response of an equivalent single-degree-of-freedom (SDOF) system whose ratio of critical damping is the same as the 10-story building in a given mode. It is clear from Fig. 4 that factors in addition to damping are involved in modifying the performance of the damper. Since the "moment arm" of the damper with respect to the base of the building is not affected by the mode shape, it appears that the main factor that influences the efficiency of the damper is the magnitude of the displacement of m, to which the damper is attached.

3 Damper Location /'. Fig. 5 illustrates some of the possible With the force imposed on mass TO 6 (i.e., k = 5), then, for the parameters in Fig. 6, the maximum reduction in amplitude will be =2 percent if the damper is attached to mi, Fig. 6 percent if the damper is placed at the top floor, Fig. 6(c). This progressive increase in efficiency with the damper "elevation" above the ground floor is to be expected since the moment arm increases with increasing j to a maximum at the top floor (i.e., j = n).

For modes other than the first, increase in moment arm must be modified to account for the relative locations of the damper and the nodal points.

4 Force Location k. Although the application of an exciting force of a given amplitude at different points along the structure will affect the magnitude of the response of the structure, it was found that the response curves for an impact damped structure are virtually independent of the point of application of the exciting force, provided that the results are expressed in terms of d/xpj and Zmaxj/xp,-.

5 Telescopic Effect. Fig. 7 shows the relation between the amplification ratio x mo x t '/xst,-, (i = 1, . . . , n) and the clearance ratio d/x s t" for the three buildings shown in Fig. 2, where each one is vibrating in its first mode and with both force and damper acting on the top floor so that j = k = n.

It is clear that if the damper mass ratio fi T is kept constant, the effectiveness of the damper is proportional to the height of the building. This is due to the fact that as the height of the building is increased, the moment arm of the impact damper is simultaneously increased to amplify the dampening effects. The deviation of the results from the corresponding ones for an equivalent SDOF system increases as the number of stories increase.

6 Damping in the Primary System (a,/3). An increase in the viscous damping of the primary system is detrimental to the efficiency of the impact damper.

7 Excitation Frequency fi. If the structure is excited at a fixed frequency, the design parameter can be optimized according to the aforementioned criteria. However, if the response is to be controlled effectively over a relatively wide frequency band, an increase in /J. T and damping or a decrease in e will improve the performance of the damper.

8 Base Excitation. The present theory can also analyze cases where the excitation is supplied through support displacement rather than force excitation. The base-excited two-degree-offreedom system in Fig. 8(a) is equivalent to the force-excited three-degree-of-freedom system in Fig. 8(6) where mi replaces the support and the excitation amplitude is adjusted so that So/Wh) = 1.

The theoretical predictions and the corresponding experimental measurements obtained by using an actual mechanical model are shown in Fig. 8(c). As mentioned earlier, of the two possible steady-state solutions, the one associated with ott, + is the stable one. This is true in spite of the fact that, for the case in the frequency range O/coi < 0.95, the solution «o + leads to larger amplitudes than a 0 ~. Note also the existence of a peak in the response curve at 0/cOi » 1/(1 + Mr)-Practical Design Considerations. Some of the design problems inherent in the application of practical impact dampers are (a) large accelerations developed during impact; (6) the accompanying noise; (c) design of the concentrated impacting mass (in massive structures, such as buildings, even a damper mass ratio of ~1 percent is still a large weight to have suspended in the structure). Analytical and experimental studies have shown that the foregoing problems can be alleviated to some extent by the use of multiple-unit impact dampers operating in parallel, instead of one single "particle."

Summary and Conclusions

An exact theory is presented for determining the steady-state motion of a damped multidegree-of-freedom lumped parameter system that is provided with an impact damper attached to some arbitrary point in the structure and is subsequently excited by a sinusoidal force, also applied to an arbitrary point.

The results were applied to several multidegree-of-freedom systems, including a model of a modern 10-story building. A detailed investigation of the effects of various parameters of the system on the response was conducted. It was found that regardless of where the excitation is applied (including the base) the efficiency of the damper increases as its moment arm in- creases. Due to the "telescoping effect," it is feasible to use impact dampers, with mass ratios on the order of 1 percent of the mass of a lightly damped structure, to considerably attenuate the response of the structure over a relatively broad suppression band around any frequency of interest. Predictions of the theory were corroborated by experimental studies with an electronic analog computer and with a mechanical model.

Journal of Applied Mechanics

Fig. 1

 1 Fig. 1 Model of system

1 H

 1 of ith mass Xj(t) = absolute velocity of ith mass z(t) = absolute displacement of damper m d y(t) = relative displacement of damper with respect to mass m,j m d = damper mass Ft, = amplitude of sinusoidol exciting force Fo sin tit ti = frequency of sinusoidol excitation Fo sin tit #o = phase angle, initially unknown CO; = ith natural frequency of multidegree-of-freedom system f t = ratio of critical damping corresponding to ith mode fi T = mass ratio of damper to total n mass of structure, m d / 2_i m i » = = mass ratio of damper to mass rtij to which it is attached, mjnij e = coefficient of restitution d = clearance in which damper is free to oscillate j = mass number to which damper is attached k = mass number to which force is applied n = total number of masses in ndegree-of-freedom system Xst; = static deflection of ith mass in an n-degree-of-freedom system when force is applied to mass m k Xpi = peak displacement amplitude of mass m,-in the absence of damper £mox; = maximum displacement amplitude of mass m v with damper operating

Fig. 2

 2 Fig. 2 Example structure: (a) four-degree-of-freedom system; (b) seven-degree-of-freedom system; (c) 10-degree-of-freedom system; (d) first three modes of 10-degree-of-freedom system

Fig. 4

 4 Fig. 4 Effect* of roods shape: / = It = n = 10, HT = 0.05, • Fig. 3 Comparison of analytical and experimental results; / = k = (a) il/on = 1; (b) fi/wj = 1; and (e) 0/«3 = 1 n = 7; n/toi = 1-(a) e = 0.075; (b) e = 0.25; and (e) e = 0.75

Fig. 5

 5 Fig. 5 Different combinations of force and damper location for 10-story building; (a) / = 1; (b) / = 5; and (e) / = 10

  (a), about 50 percent if the damper is at m 5 , Fig-6(6), and more than 90

Fig. 6 Fig. 7 "

 67 Fig. 6 Effects of force and damper location: n = 10, Qfui = 1, HT = 0.05, e = 0.25; (a) / = 1; (b) / = 5; and (e) / = 10

Fig. 8

 8 Fig. 8 Bate excitation: (a) base-excited two-degree-of-freedom model; (b) equivalent force excited model; and (c) frequency response of a base-excited two-degree-of-freedom system

Acknowledgment

This work was supported by the National Science Foundation Grant No. GK-11553. The assistance of Dr. D. S. Margolias in the preparation of this manuscript is appreciated.


Thibaut Guirimand

Charlie Pauvert

Sandra Derozier

Anne-Laure Abraham

Mahendra Mariadassou

Valentin Loux

Pierre Renault

Cheese ecosystems

come

Genes Annotations


L Roques

E Walker

P Franck

S Soubeyrand

E K Klein

Using genetic data to estimate diffusion rates in heterogeneous landscapes

Keywords: Reaction-diffusion, stochastic differential equation, inference, mechanisticstatistical model, allele frequencies, genotype measurements

Having a precise knowledge of the dispersal ability of a population in a heterogeneous environment is of critical importance in agroecology and conservation biology as it can provide management tools to limit the effects of pests or to increase the survival of endangered species. In this paper, we propose a mechanistic-statistical method to estimate space-dependent diffusion parameters of spatially-explicit models based on stochastic differential equations, using genetic data. Dividing the total population into subpopulations corresponding to different habitat patches with known allele frequencies, the expected proportions of individuals from each subpopulation at each position is computed by solving a system of reaction-diffusion equations. Modelling the capture and genotyping of the individuals with a statistical approach, we derive a numerically tractable formula for the likelihood function associated with the diffusion parameters.

In a simulated environment made of three types of regions, each associated with a different diffusion coefficient, we successfully estimate the diffusion parameters with a maximumlikelihood approach. Although higher genetic differentiation among subpopulations leads to more accurate estimations, once a certain level of differentiation has been reached, the finite size of the genotyped population becomes the limiting factor for accurate estimation.

Introduction

Dispersal is one of the main forces driving population redistribution [START_REF] Turchin | Quantitative Analysis of Movement: Measuring and Modeling Population Redistribution in Animals and Plants[END_REF], gene flow [START_REF] Slatkin | Gene flow and the geographic structure of natural populations[END_REF][START_REF] Bohonak | Dispersal, gene flow, and population structure[END_REF]) and genetic diversity [START_REF] Hewitt | The genetic legacy of the quarternary ice ages[END_REF][START_REF] Roques | Allee effect promotes diversity in traveling waves of colonization[END_REF][START_REF] Roques | The effect of competition on the neutral intraspecific diversity of invasive species[END_REF]. Dispersal directly affects population flows between different spatial positions. Having a precise knowledge of these flows is of critical importance in agroecology and conservation biology, as it can provide management tools to limit the effects of pests [START_REF] Gilligan | Sustainable agriculture and plant diseases: an epidemiological perspective[END_REF][START_REF] Papaïx | Influence of cultivated landscape composition on variety resistance: an assessment based on wheat leaf rust epidemics[END_REF] or to increase the survival of endangered species [START_REF] Hanski | Metapopulation Biology: Ecology, Genetics, and Evolution[END_REF] by acting on the landscape structure.

In heterogeneous environments, dispersal is often influenced by local landscape features. This local effect of the landscape on mobility can be captured by most spatially-explicit models at the scale of individual movement, such as in random walk models and stochastic differential equations [START_REF] Roques | Modeling animal movements using stochastic differential equations[END_REF][START_REF] Smouse | Stochastic modelling of animal movement[END_REF] or at the scale of population density in reaction-diffusion models [START_REF] Shigesada | Biological Invasions: Theory and Practice[END_REF][START_REF] Cantrell | Spatial ecology via reaction-diffusion equations[END_REF][START_REF] Ovaskainen | Bayesian methods for analyzing movements in heterogeneous landscapes from mark-recapture data[END_REF][START_REF] Roques | Modèles de réaction-diffusion pour l'écologie spatiale[END_REF], through a space-dependent mobility parameter.

Traditionally, Euclidian distances and least-cost distances have been used, e.g. in models of isolation by distance [START_REF] Wright | Isolation by distance[END_REF][START_REF] Rousset | Genetic differentiation and estimation of gene flow from f-statistics under isolation by distance[END_REF][START_REF] Broquet | Genetic isolation by distance and landscape connectivity in the american marten (martes americana)[END_REF]), to quantify the population flows. However, these approaches have several drawbacks; in particular, least-cost distances are based on a subjective definition of a cost function and assume a single and optimal migration path. Resistance-based approaches are more realistic [START_REF] Mcrae | Isolation by resistance[END_REF][START_REF] Graves | Estimating landscape resistance to dispersal[END_REF], as the resistance distance is computed in a similar manner as in an electrical network, where all possible paths are taken into account. However, the computation of resistance distances is time-consuming and generally does not allow for a precise estimation of the local resistance parameter, e.g. by maximum likelihood, but rather to test a limited set of conjectured resistance values [START_REF] Graves | Current approaches using genetic distances produce poor estimates of landscape resistance to interindividual dispersal[END_REF], leading again to a subjective parametrization of the local resistance values. Additionally, although there exist random walk interpretations of the effective resistance in electrical networks [START_REF] Doyle | Random walks and electric networks[END_REF], these approaches are not based on a mechanistic description of the spatio-temporal dynamics of a population (see Theorems 1 and 2 in [START_REF] Tetali | Random walks and the effective resistance of networks[END_REF], and therefore do not directly quantify population flows.

The method developed in this paper enables direct and fast estimations of a spatiallyheterogeneous parameter D(x) measuring the local mobility of individuals at the space position x, in mechanistic models based on stochastic differential equations. More precisely, we assumed that the individual trajectories followed Itō diffusion processes, corresponding to uncorrelated random walks with spatially-varying speed. This framework is widely used for analyzing movement, see [START_REF] Roques | Modeling animal movements using stochastic differential equations[END_REF]; [START_REF] Smouse | Stochastic modelling of animal movement[END_REF] and the references therein. In this framework, the expected population density at any time and space position can be computed using the corresponding Fokker-Planck partial differential equation [START_REF] Gardiner | Stochastic methods[END_REF] with diffusion parameters D(x). Here, we also took into account death events occurring at exponentially distributed times, leading to a reaction-diffusion description of the population density. Thanks to a well-developed theory for their numerical analysis and efficient softwares (e.g., freefem++, see [START_REF] Hecht | New development in freefem++[END_REF], or Comsol Multiphysics c ), the numerical computation of the solutions of such reaction-diffusion equations is fast and reliable, even in the presence of heterogeneous coefficients. This makes them ideally suited for parameter estimation [START_REF] Soubeyrand | Parameter estimation for reaction-diffusion models of biological invasions[END_REF].

Several types of data can be used for the estimation of dispersal parameters in population models. Most studies bear on abundance data [START_REF] Roques | A statistical-reaction-diffusion approach for analyzing expansion processes[END_REF] or on mark-recapture experiments [START_REF] Turchin | Quantitative Analysis of Movement: Measuring and Modeling Population Redistribution in Animals and Plants[END_REF][START_REF] Ovaskainen | Bayesian methods for analyzing movements in heterogeneous landscapes from mark-recapture data[END_REF], where individuals are marked with different technics, such as color markers or radioactive isotopes [START_REF] Southwood | Ecological methods[END_REF]. As opposed to this marking experiments, passive surveys of the relative abundance of genetic markers that are naturally present in a population lead to spatio-temporal data which can be easier to obtain and more informative (Robledo-Arnuncio 2012). Here, we considered the dispersal of individuals starting from several habitats and which eventually died, defining the end of the dispersal period. The estimation of the parameter D(x) was based on measurements of the genotypes of individuals captured at several positions during the dispersal process, the location Comment citer ce document : Roques, L., Walker, E., Franck, P., Soubeyrand, S., Klein, E. (2016). Using genetic data to estimate diffusion rates in heterogeneous landscapes. Journal of Mathematical Biology, 73 (2), 397-422. DOI : 10.1007/s00285-015-0954-4

of the habitats and the allele frequencies before dispersal in these sites being considered as known. Although numerous methods have been developed in landscape genetics to estimate dispersal from molecular markers [START_REF] Hamrick | Using population genetic analyses to understand seed dispersal patterns[END_REF] we are not aware of any article estimating diffusion parameters from such data.

The method is based on a mechanistic-statistical approach [START_REF] Berliner | Physical-statistical modeling in geophysics[END_REF][START_REF] Wikle | Hierarchical models in environmental science[END_REF][START_REF] Ovaskainen | Bayesian methods for analyzing movements in heterogeneous landscapes from mark-recapture data[END_REF][START_REF] Soubeyrand | Spatio-temporal structure of host-pathogen interactions in a metapopulation[END_REF][START_REF] Roques | A statistical-reaction-diffusion approach for analyzing expansion processes[END_REF][START_REF] Soubeyrand | Parameter estimation for reaction-diffusion models of biological invasions[END_REF] included in the framework of state-space models [START_REF] Patterson | State-space models of individual animal movement[END_REF][START_REF] Durbin | Time series analysis by state space methods[END_REF]. This approach typically combines a mechanistic model describing the dynamics under investigation with a statistical model conditional on the dynamics, describing how the measurements have been collected, bridging the gap between the data and the model for the dynamics. In the mechanistic part of our model, we divided the total population into subpopulations at Hardy-Weinberg equilibrium, each one corresponding to the individuals coming from a different habitat patch. The dynamics of the different subpopulations were then described by a system of reaction-diffusion equations, as in [START_REF] Roques | Allee effect promotes diversity in traveling waves of colonization[END_REF]. Given a diffusion coefficient D(x), this allowed us to compute the expected proportions of individuals from each subpopulation at each space position. Conversely, the genotype data contain information about these proportions; namely, the probability to observe a given genotype at some trapping location depends on the respective contributions of each subpopulation and on the allele frequencies in these subpopulations. Modelling the capture of the individuals with a statistical approach, and using the Hardy-Weinberg equilibrium assumption in each subpopulation as a key ingredient, it was therefore possible to derive a numerically tractable formula for the likelihood function associated with the diffusion parameters D(x), given the genotypes of the captured individuals.

Note: a summary of the notations used throughout this paper is provided in Table 1.

Models

Modelling dispersal and death. We begin with a Lagrangian description of the individual movements. We assumed that the positions of the individuals followed 2-dimensional spaceheterogeneous Itō diffusion processes, corresponding to uncorrelated random walks. This means that the individuals travel at random, with no drift in any particular direction [START_REF] Roques | Modeling animal movements using stochastic differential equations[END_REF][START_REF] Smouse | Stochastic modelling of animal movement[END_REF]. However, the mobility of the individuals can be influenced by their position. The corresponding stochastic differential equation for the position X t ∈ R 2 of an individual at time t can be written: We also assumed that each individual had a life expectancy ν > 0, the death events being independent and identically distributed and modelled by exponential distributions with parameter 1/ν.

dX t = 2D(X t )
Under these assumptions, we can switch to an Eulerian description of the population. The expected population density u(t, x) at time t and position x, starting from an initial distribution u 0 satisfies the following Fokker-Planck reaction-diffusion equation (see e.g. [START_REF] Gardiner | Stochastic methods[END_REF]:

∂u ∂t = ∆(D(x) u) - u ν , t > 0, x ∈ Ω, u(0, x) = u 0 (x), (2.2)
where

∆ = ∂ 2 ∂x 2 1 + ∂ 2 ∂x 2 2
is the 2-dimensional Laplace diffusion operator. The set Ω ⊂ R 2 is the study region. For this equation to be well-posed, some conditions on the boundary ∂Ω of Ω have to be specified. In the computations of this paper, we assumed absorbing conditions (u(t, x) = 0 on ∂Ω). Reflecting conditions (∇(D(x)u(t, x)) • n(x) = 0 on ∂Ω, where n(x) is the outward normal to Ω) could have been assumed as well.

The quantity u 0 is called the pre-dispersal density. It corresponds to the density of individuals at a motionless stage, such as eggs, pupae or larvas in insect populations or immature seeds in plants. The quantity u(t, x) is the (expected) density of dispersers, e.g. insects at the adult stage or dispersing seeds. We assumed a zero pre-dispersal density outside some known disjoint subsets (habitats) Ω h , h = 1, . . . , H, where u 0 is positive and constant (equal to α > 0):

u 0 (x) = H h=1 α 1 x∈Ω h , for all x ∈ Ω, (2.3)
where 1 x∈Ω h is the characteristic function of the set Ω h : it takes the value 1 in Ω h and 0 anywhere else. A possible gradual release of the pre-dispersal populations could also be assumed by considering a slightly modified version of the equation (2.2), see Appendix A.

Modelling the capture of individuals from different sources. We considered the case of nonattractive traps, corresponding to disjoint sets θ τ , τ = 1, . . . , J in the study site. We assumed that the expected number of individuals C τ captured in a trap θ τ was proportional to the cumulated population in θ τ : (2.4) and β τ the capture rate (number of captured individuals per unit of time per unit of area in the trap θ τ ). Note that u(t, x) ≤ e -t/ν max u 0 ; this means that the value of w ∞ (x) can be precisely

C τ = β τ θτ w ∞ (x) dx, with w ∞ (x) = ∞ 0 u(t, x) dt,
Comment citer ce document : Roques, L., Walker, E., Franck, P., Soubeyrand, S., Klein, E. (2016). Using genetic data to estimate diffusion rates in heterogeneous landscapes. Journal of Mathematical Biology, 73 (2), 397-422. DOI : 10.1007/s00285-015-0954-4

approached by computing the above integral over a finite interval. For sufficiently small traps, w ∞ (x) can be considered constant in θ τ , which leads to:

C τ = β τ |θ τ | w ∞ (x τ ), (2.5)
where x τ is the location of the center of the trap θ τ , and |θ τ | is the trap area.

Remark 1 With this approach, the trapping process has no influence on the species dynamics, i.e., the trapped individuals are not removed from the system. To avoid this lack of realism, a sink term -β τ u1 x∈θτ could be added to the right-hand side of (2.2). For the sake of simplicity, we assumed that the traps were small enough to consider that this term could be neglected.

Consider now the density u h (t, x) of individuals coming from a given habitat Ω h . Since all the individuals are supposed to share the same dispersal and death characteristics independently of their origin, the densities u h (t, x) satisfy [START_REF] Roques | Allee effect promotes diversity in traveling waves of colonization[END_REF]:

∂u h ∂t = ∆(D(x) u h ) - u h ν , t > 0, x ∈ Ω, (2.6) and u h (0, x) = u h 0 (x) = α 1 x∈Ω h , for all x ∈ Ω, (2.7)
with the same boundary conditions as the total population u. The dynamics of the different fractions u h of the total population u is therefore described by a system of H decoupled reactiondiffusion equations. Summing up all these equations, it can be checked that for all t, x,

u(t, x) = H h=1 u h (t, x).
The expected number of individuals coming from a habitat Ω h and which are captured in a trap θ τ is then given by:

C h τ = β τ |θ τ | w h ∞ (x τ ), with w h ∞ (x) = ∞ 0 u h (t, x) dt. (2.8)
We assumed that the population was large enough so that the number of captured individuals was larger than a constant G in all traps. This constant corresponds to the number of individuals genotyped in each trap.

Parameters and data

Our goal was to estimate the diffusion parameters D(x) for all x ∈ Ω. Other unknown parameters are the pre-dispersal density in the habitats (α) and the capture rates in the different traps (β τ ). We assumed that the global population before dispersal, u 0 , was organised into several subpopulations each of which was at Hardy-Weinberg equilibrium and linkage equilibrium among loci. For the sake of clarity, we assumed that there were exactly H subpopulations, with densities u h 0 = α 1 x∈Ω h , each one corresponding to a habitat Ω h . More complex assumptions are also possible, see Remark 3.

The positions of the habitats Ω h and of the traps θ τ were known. For each subpopulation h and each locus λ (e.g. microsatellites) out of Λ loci, the pre-dispersal frequencies of A λ alleles were known and designated as:

F hλ = (p hλa ) a=1,...,A λ .
(3.9)

The individuals captured in θ τ were genotyped at the same Λ loci. These individuals were assumed to be diploid; thus, each genotype was described by:

G = (a 1 λ , a 2 λ ) λ=1,...,Λ .
Comment citer ce document : Roques, L., Walker, E., Franck, P., Soubeyrand, S., Klein, E. (2016). Using genetic data to estimate diffusion rates in heterogeneous landscapes. Journal of Mathematical Biology, 73 (2), 397-422. DOI : 10.1007/s00285-015-0954-4

4 Computation of the likelihood

The computation of genotype likelihoods builds on a combination of classical genetic assignment studies [START_REF] Paetkau | Microsatellite analysis of population structure in canadian polar bears[END_REF][START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF] and seed dispersal analyses from trap data (Robledo-Arnuncio and Garcia 2007; [START_REF] Klein | Seed dispersal kernels estimated from genotypes of established seedlings: does density-dependent mortality matter?[END_REF]). Among the individuals captured in a trap θ τ , G individuals have been genotyped. This led to G genotypes G iτ , i = 1, . . . , G.

The conditional probability that an individual i carries alleles (a 1 , a 2 ) ∈ {1, . . . , A λ } 2 at locus λ, given that this individual comes from a habitat Ω h , can be deduced from the allele frequencies in subpopulation h (see Section 3). The two alleles being independent, which follows from the Hardy-Weinberg equilibrium assumption in the subpopulation u h 0 , we get:

P((a 1 , a 2 )|Ω h ) = 2 k λ p hλa 1 p hλa 2 , (4.10)
where k λ = 0 if the individual is homozygous at locus λ (a 1 = a 2 ) and k λ = 1 otherwise. Using the linkage equilibrium assumption among loci, we get the conditional probability of genotype G iτ :

P(G iτ |Ω h ) = 2 ki Λ λ=1 p hλa 1 p hλa 2 , (4.11)
where k i is the number of heterozygous loci in the genotype G iτ .

Remark 2 For the sake of simplicity, the dependence of a 1 and a 2 with respect to the locus λ and the individual i have been dropped in our notations. For instance, in formula (4.11), a 1 and a 2 may designate different alleles, depending on the locus λ and on the individual i.

The law of total probability leads to:

P(G iτ ) = 2 ki H h=1 Λ λ=1
p hλa 1 p hλa 2 P(indiv. i comes from Ω h ). (4.12)

We have seen in Section 2 that the expected number of individuals trapped in θ τ was given by C τ , and the expected number of individuals coming from a habitat Ω h and which are captured in a trap θ τ was given by C h τ . Let us denote by I h τ the number of individuals coming from Ω h and captured in θ τ , and let us set

I τ = S h=1 I h τ ≥ G, (4.13)
the total number of individuals captured in θ τ . Assume that I h τ follows a Poisson distribution:

I h τ ∼ P(C h τ ). (4.14)
Thus, I τ also follows a Poisson distribution with parameter

C τ = S h=1 C h τ .
It can be verified that the conditional distribution of I h τ given I τ satisfies a binomial distribution with parameters I τ and C h τ /C τ . Thus, the conditional expectation of the proportion

I h τ /I τ given I τ is C h τ /C τ , which is independent of I τ .
Finally, this shows that:

E I h τ I τ = C h τ C τ .
The genotyping process corresponds to a sampling without replacement, this means that the number of genotyped individuals coming from habitat Ω h follows a multivariate hypergeometric Comment citer ce document : Roques, L., Walker, E., Franck, P., Soubeyrand, S., Klein, E. (2016). Using genetic data to estimate diffusion rates in heterogeneous landscapes. Journal of Mathematical Biology, 73 (2), 397-422. DOI : 10.1007/s00285-015-0954-4 distribution with parameters I τ , I h τ and G. For large values of I τ , this distribution converges to the multinomial distribution with parameters G and (I 1 τ /I τ , . . . , I H τ /I τ ). Using this multinomial distribution, we can compute the probability that a genotyped individual i trapped in θ τ comes from a habitat h:

P(indiv. i comes from Ω h ) = E I h τ I τ = C h τ C τ . (4.15)
The hypergeometric distribution would lead to the same formula, but the advantage of the multinomial distribution is that it guarantees an independence assumption between the individuals trapped at a same location.

Assuming that the trapping and genotyping processes are independent, and using formulas (4.12) and (4.15), we finally get the likelihood function associated with the unknown parameter D, given the genotypes G iτ :

L(D) = τ =1,...,J i=1,...,G P(G iτ |D) = τ =1,...,J i=1,...,G P(G iτ |C h τ /C τ ) = 2 k τ =1,...,J i=1,...,G H h=1 C h τ Cτ Λ λ=1 p hλa 1 p hλa 2 , (4.16)
where k is the total number of heterozygous loci in the genotyped population. Coming back to the definitions (2.5) and (2.8) of C τ and C h τ , we can compute the ratio:

C h τ C τ = w h ∞ (x τ ) w ∞ (x τ ) . (4.17)
This shows that C h τ /C τ is independent of the capture rates β τ . From the linearity of the equations (2.2) and (2.6), it follows that w h ∞ and w ∞ are proportional to α which means that the ratio C h τ /C τ is also independent of the choice of α. Thus, the likelihood L(D) can be computed with an arbitrary choice of parameters, e.g. α, β τ = 1. Note that if the source intensity α was spatially variable, it would not simplify in the expression C h τ /C τ . In such case, α should be estimated in each source, which is only possible up to a multiplicative constant since C h τ would be proportional to the value of α in each habitat h.

Remark 3 We recall that u 0 was organised into several subpopulations at Hardy-Weinberg equilibrium. Here, we assumed that, before dispersal, these subpopulations coincided with the habitats, which were disjoint. A first easy generalization would be to consider several subpopulations in the same habitat. Assume that there are R subpopulations with pre-dispersal densities

u r 0 (x) = H h=1
µ hr 1 x∈Ω h , and that the allele frequencies F rλ are known in these subpopulations.

In this case, the probability P(G iτ |Ω h ) does not satisfy formula (4.11) but can be computed as follows:

P(G iτ |Ω h ) = R r=1
µ hr P(G iτ |indiv. i comes from subpop. r), (4.18)

which leads to:

P(G iτ |Ω h ) = 2 ki R r=1 µ hr Λ λ=1 p rλa 1 p rλa 2 . (4.19)

Numerical computations

The aim of this section was to validate the maximum likelihood estimator of Section 4 on a simulated data set.

Comment citer ce document : Roques, L., Walker, E., Franck, P., Soubeyrand, S., Klein, E. (2016). Using genetic data to estimate diffusion rates in heterogeneous landscapes. Journal of Mathematical Biology, 73 (2), 397-422. DOI : 10.1007/s00285-015-0954-4

Simulated data set

Landscape. The study site Ω was a unit square (0, 1) × (0, 1) containing H = 6 habitats Ω h , described by balls B R (x h ) of radius R = 0.05 and centers x h ∈ Ω and a rectangular region Q = (q -R/2, q + R/2) × (0, 1) modelling a barrier to dispersal (q = 0.5). The rest of the study site was considered as the matrix (Fig 1). The general idea was to consider diffusion parameters D 1 in the matrix, D 2 = D 1 /2 in the habitats and D 3 = D 1 /10 in the barrier. However, for the well-posedness of the reactiondiffusion equations for u and u h , the coefficient D(x) had to be positive and smooth. We thus defined the heterogeneous diffusion parameters D(x) as:

D(x) = exp d 1 + d 2 H h=1 φ(x -x h ) + d 3 ψ(x) ,
(5.20)

for smooth positive functions φ and ψ, such that φ(x-x h ) was compactly supported in B 2 R (x h ) for any h = 1, . . . , S and ψ was compactly supported in (q -R, q + R) × (0, 1) and max φ = max ψ = 1. The precise shape of φ and ψ is detailed in Appendix B. The numerical values of d 1 , d 2 and d 3 that we used in our computations were: for random walk movements with one direction change per minute this value of the diffusion parameter corresponds to a flying speed of 53m per minute. Under the same assumptions, D = 10 -2 /2 and D = 10 -2 /10 correspond to flying speeds of 37m and 17m per minute, respectively. See [START_REF] Kareiva | Local movement in herbivorous insects: Applying a passive diffusion model to markrecapture field experiments[END_REF] for reference values of diffusion parameters of several insect species.

d 1 = log(0.01), d 2 = -log (
Allele frequencies. The number of loci and the number of alleles at each locus were Λ = 10 and A = 10. Following [START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF], we drew the allele frequencies in a flat Dirichlet distribution, with concentration parameter q in each subpopulation h and for each locus λ :

F sλ ∼ D(q).
(5.22)

Values of q larger than 1 lead to evenly distributed frequencies, whereas smaller values of q lead to distributions which are concentrated on a few components. In our simulations, we adjusted the value of q such that the fixation index F ST among subpopulations (see Appendix C) was close to 0.01, 0.05 or 0.1 (Table 2), with a relative tolerance of 0.1%.

Simulation of the measured data. We solved the reaction-diffusion models for the cumulated population densities w ∞ and w h ∞ , for h = 1, . . . , 6 with the true diffusion parameter values (5.21) and ν = 5 for the life expectancy of the dispersers (see Appendix D). The proportions (C h τ /C τ ) h=1,...,6 have been computed using formulas (2.5) and (2.8). We recall that, from formulas (4.15) and (4.17

), C h τ /C τ = w h ∞ (x τ )/w ∞ (x τ )
is the probability that an individual trapped in θ τ comes from the habitat Ω h . The probability w h ∞ (x)/w ∞ (x) can be computed at each point in Ω; it is depicted in Fig. 2. Movies of the dynamics of the probability w h t (x)/w t (x) that an individual trapped in θ τ between the times 0 and t comes from the habitat Ω h are available as supplementary materials.

In each trap θ τ (see Fig. 1 for the locations of the traps), the numbers of genotyped individuals coming from the habitats Ω h followed a multinomial distribution with parameters G and (C h τ /C τ ) h=1,...,6 . We tested the effect of the number of traps by using J = 10 or 20 traps and of

Comment citer ce document : Roques, L., Walker, E., Franck, P., Soubeyrand, S., Klein, E. (2016). Using genetic data to estimate diffusion rates in heterogeneous landscapes. Journal of Mathematical Biology, 73 (2), 397-422. DOI : 10.1007/s00285-015-0954-4 the total number of genotyped individuals by using J × G = 500 or 2000 genotyped individuals (Table 2).

For each genotyped individual i, the genotype G iτ was randomly drawn according to the allele frequencies F hλ in its habitat of origin. The simulated observations consisted in the genotypes G iτ , i.e., the set:

M = {G iτ , i = 1, . . . , G, τ = 1, . . . , J}.
(5.23)

For each of the 12 sets of parameters (F ST , J, J G), we simulated 70 data sets M. For each data set, the estimator ( d1 , d2 , d3 ) of (d 1 , d 2 , d 3 ) has been obtained by minimizinglog(L(D)) (see formula (4.16)). The minimization was performed using the Matlab constrained gradientbased minimization algorithm fmincon, with the constraint

( d1 , d2 , d3 ) ∈ (d 1 -5, d 1 + 5) × (d 2 -5, d 2 + 5) × (d 3 -5, d 3 + 5).
(5.24)

In the computation oflog(L(D)), the numerical evaluation of the quantities (C h τ /C τ ) h=1,...,6 was based on the finite element method as described in Appendix D. The average computation time for one estimation was about 45min with a dual Core Intel processor, while the computation of the likelihood took about 5sec. The Bayesian method [START_REF] Marin | Bayesian Core[END_REF] would be more computationally intensive, but it could also be used to compute a posterior distribution of D(x) at each position x.

Results

The typical profile of the log of the likelihood function (4.16) is depicted in Fig 3 for one simulation with F ST = 0.05, J = 20 traps and G J = 2000 individuals. The likelihood tends to decay as ( d1 , d2 , d3 ) diverges from the true value (d 1 , d 2 , d 3 ), suggesting an efficient parameter estimation by maximum likelihood, even with the constraint (5.24) was relaxed.

The direct analysis of the quality of the estimator ( d1 , d2 , d3 ) is presented in Table 2, in terms of the F ST value, the number of traps and the number of genotyped individuals. In this section, we focus on the more biologically meaningful estimators D1 = exp d1 , D2 = exp d1 + d2 , D3 = exp d1 + d3 which can be directly compared to the values D 1 , D 2 , D 3 of the diffusion parameter in each of the three regions (matrix far from the other regions, center of the habitats and center of the barrier, respectively).

For each F ST value (0.01, 0.05 and 0.1), with J = 20 traps and J G = 2000 genotyped individuals, we observe in Fig. 4 that the median of D is correctly centered on the true value of D in each of the three regions (see Fig. 1 and Section 5.1 for the definition of the three regions). We also note that there are no outliers among the estimated values of D in the matrix (the largest region in Fig. 1) for all F ST values. Outliers far from the true value of D (about 12

Comment citer ce document : Roques, L., Walker, E., Franck, P., Soubeyrand, S., Klein, E. (2016) times the true value) appear in the other regions for small F ST values (F ST =0.01). For larger F ST values, the number of outliers is reduced and they are closer to the true value (2 outliers in the habitats and the barrier for F ST =0.05 and only 2 outliers in the habitats region for F ST =0.1). In the three regions (matrix/habitats/barrier), the interquartile ranges are reduced as the F ST is increased from 0.01 to 0.05; this decrease of the interquartile ranges by a factor 2 in average indicates a strong effect of the F ST on the accuracy of the estimation. There is no clear difference between the interquartile ranges corresponding to F ST = 0.05 and F ST = 0.1; however, the bias and the standard deviation of ( d1 , d2 , d3 ) are still improved when the F ST is increased from 0.05 to F ST = 0.1 (see Table 2). Table 2 Effect of the F ST , the number of traps and of the total number of genotyped individuals on the quality of the estimator ( d1 , d2 , d3 ).

FS T = 0 .0 1 FS T = 0 .0 5 FS T = 0 .1
To better quantify the power of our genetic system to discriminate between different habitats of origin, depending on the F ST , we computed the average posterior probability that an individual comes from its true habitat, say Ω h * , given its genotype G, that is:

P(i comes from Ω h * |G).
Using Bayes theorem and considering that all of the H habitats are equally likely a priori, we get:

P(i comes from Ω h * |G) = P(G|i comes from Ω h * )/H H h=1 P(G|i comes from Ω h )/H . (5.25)
For 300 F ST values between 0 and 0.1, we simulated 1000 data sets (allele frequencies in the H sites, see Section 5.1), we sampled 1000 genotypes G in one of the sites, and we averaged the quantity P(i comes from Ω h * |G) over the 1000 individuals. We call the obtained quantity the Comment citer ce document : Roques, L., Walker, E., Franck, P., Soubeyrand, S., Klein, E. (2016)

F ST P (Ω s * |G)
H=6, Λ=10, A=10 H=3, Λ=10, A=10 H=12, Λ=10, A=10 Fig. 5 Discrimination power: average posterior probability that an individual comes from its true habitat, given its genotype G, in terms of the F ST . Left: effect of the number of loci Λ and of the number of alleles per locus A. Right: effect of the number of habitats H. In both figures the red curve corresponds to the parameters that have been used in this study.

discrimination power of our genetic system. As expected, increasing the F ST leads to higher discrimination power (Fig. 5). With our parameters (H = 6, Λ = 10, A = 10), the discrimination power is about 0.5 for an F ST of 0.01; it reaches 0.9 for an F ST of 0.05 and 0.99 for an F ST of 0.1. This explains the gap between the quality of the estimators obtained with the F ST values larger than 0.05, compared to the case F ST = 0.01. This could also explain the small difference in the accuracy of the estimation when the F ST is increased to 0.1, compared to the case F ST = 0.05. Additionally, we note that it does not seem necessary to go beyond the value F ST =0.1 to get a very high discrimination power. This also means that the variability in our estimator D, when F ST =0.1, is not due to an uncertainty on the habitat of origin of the genotyped individuals, but rather to the sampling variability due to the limited number of these individuals.

Fig. 5 also shows the effect of the number of habitats, the number of loci and the number of alleles. Increasing the number of alleles per locus or the number of loci have a comparable effect, with in both cases an increase in the discrimination power. For F ST values close to 0, the discrimination power converges to 1/H, meaning that very small genetic differentiation between the subpopulations leads to equiprobable habitat of origin. In such case, the data contain no information on the origin of the trapped individuals, and the estimation of D is therefore not possible. For larger values of the F ST , increasing the number of sites always leads to lower discrimination power.

The results in Fig. 6 show, for F ST = 0.1, the effect of the number of traps and of the total number of genotyped individuals on the quality of the estimator. In all cases, the median of D is close to the true value, in each of the three regions. For a fixed total number of genotyped individuals, increasing the number of traps J from 10 to 20 does not increase significantly the quality of the estimator. It can even decrease it in some situations, e.g., the interquartile range is about twice larger with 20 traps than with 10 traps in the habitats. A possible explanation for this counterintuitive effect is that, when the number of traps is increased, the number G of individuals genotyped per trap decreases; in the regions where no traps are added (here, the habitats, see Fig. 1), this can lead to more uncertainty on the estimator.

The results of Fig. 6 confirm that increasing the number of genotyped individuals leads to far better estimations: in both cases J = 10 and J = 20, the interquartile ranges are divided by more than two in average when the total number of genotyped individuals is increased from J ×G = 500 to J ×G = 2000, and the distance of the outliers to the median (i.e., approximatively to the true value) is reduced.

Discussion

Using broadly-recognized population models based on a mechanistic description of individual movements [START_REF] Roques | Modeling animal movements using stochastic differential equations[END_REF][START_REF] Smouse | Stochastic modelling of animal movement[END_REF], we have developed an approach to estimate the local effect of the environment on individual mobility, based on genetic data. In an Comment citer ce document : Roques, L., Walker, E., Franck, P., Soubeyrand, S., Klein, E. (2016). Using genetic data to estimate diffusion rates in heterogeneous landscapes. Journal of Mathematical Biology, 73 (2), 397-422. DOI : 10.1007/s00285-015-0954-4

Using genetic data to estimate diffusion rates in heterogeneous landscapes 13 environment made of three types of regions, each one associated with a different level of mobility -or diffusion -we successfully estimated the diffusion parameters D(x) in each region.

The reaction-diffusion framework enabled a fast computation of expected population densities, making parameter estimation possible in a reasonable time. Genetic data had already proved their effectiveness in the estimation of dispersal kernels (Robledo-Arnuncio and Garcia 2007; [START_REF] Klein | Seed dispersal kernels estimated from genotypes of established seedlings: does density-dependent mortality matter?[END_REF]) in more empirical models. Our results show that successful estimation of parameters of mechanistic population models is also possible using genotype measurements during dispersal and allele frequency data before dispersal, without the need of abundance or mark-recapture data. Genetic data lead to likelihoods of arriving from a given source, and are insensitive to the global population size (i.e., to the parameter α in our approach) and to the relative efficiency of the different traps (the parameters β τ ). This advantage of working on probabilities of originating from the different sources has already been shown for kernel estimations (Robledo-Arnuncio and Garcia 2007; [START_REF] Klein | Seed dispersal kernels estimated from genotypes of established seedlings: does density-dependent mortality matter?[END_REF]. A related approach, proposed by [START_REF] Ovaskainen | Bayesian methods for analyzing movements in heterogeneous landscapes from mark-recapture data[END_REF] allowed to estimate the parameters of a diffusion model, based on mark-recapture data, with a single type of marks. In our framework, the genotype information, given the allele frequencies in the different habitats constituting our study-site, can be seen as mark-recapture data, with several types of marks (one per habitat) and some uncertainty on the marks of the captured individuals. Mark-recapture experiments with several types of marks should lead to good estimation results, as they would combine the advantages of our method (insensitivity to several parameters) and of traditional mark-recapture experiments (no unknown external sources, perfect knowledge of the frequencies).

The genetic differentiation between the subpopulations corresponding to the different habitats of origin plays a key role in the quality of the estimation, as shown by the strong effect of the F ST , especially on the variability of the estimator. Low F ST values are associated with larger standard deviations and interquartile ranges, which can be explained by a lower posterior probability associated to the true habitat of origin of a genotyped individual. Defining the discrimination power of our genetic system as the average of this posterior probability, we could disentangle the effect of the F ST and that of the other sources of uncertainty in our estimator of the diffusion parameters D(x). With an intermediate level of genetic differentiation (F ST = 0.05), the discrimination power was high (0.9 with 20 traps and 2000 genotyped individuals), and the quality of our estimator of D(x) was comparable to the case F ST = 0.1. With an F ST of 0.1, the discrimination power was close to 1, which is almost exact; the remaining uncertainty in the determination of D(x) may therefore be sampling variance due to the finite number of genotyped individuals per trap. It could also be due to the lack of uniqueness in the inverse problem of determining D(x), even with infinite population sizes. From a theoretical Comment citer ce document : Roques, L., Walker, E., Franck, P., Soubeyrand, S., Klein, E. (2016). Using genetic data to estimate diffusion rates in heterogeneous landscapes. Journal of Mathematical Biology, 73 (2), 397-422. DOI : 10.1007/s00285-015-0954-4 viewpoint, the unique determination of diffusion and conductivity coefficients based on a finite set of measurements is a difficult problem, as illustrated by the Calderón problem [START_REF] Calderón | On an inverse boundary value problem[END_REF] of determining the electrical conductivity of a medium for which uniqueness is only proved with infinitely many observations [START_REF] Sylvester | A global uniqueness theorem for an inverse boundary value problem[END_REF][START_REF] Nachman | Global uniqueness for a two-dimensional inverse boundary value problem[END_REF].

As expected, the size of the post-dispersal sample has a clear effect on the uncertainty of the estimation. The quality of the estimation, and especially its variability, is clearly improved when the number of genotyped individuals is increased. If the total number of genotyped loci and the number of traps were fixed, the trade-off between increasing the number of loci and the number of genotyped individuals per trap would depend on the main source of uncertainty: increasing the number of loci per individual increases the discrimination power while increasing the number of genotyped individuals per trap reduces the sampling variance. The role of the number of traps is less obvious. Intuitively, increasing the number of traps leads to a better coverage of the study site, which should have a positive impact on the estimation. However, with a fixed total number of genotyped individuals, this leads to a decrease in the number of genotyped individuals per trap and can therefore produce more uncertainty in the regions where the number of traps has not been increased (the habitats in our simulations). It should be noted however that too few traps may lead to identifiability problems, as would in the extreme case of a unique trap placed at equidistance between two habitats in an otherwise homogeneous landscape.

Remarkably, the estimation of D(x) in the matrix remains accurate, with no outliers among the estimators of D(x) even with low F ST values. This may be the consequence of the larger area of the matrix compared to the other regions in the landscape leading to a stronger effect of the value of the diffusion parameter in this region. Conversely, larger standard deviations and interquartile ranges are observed in the habitats, for all F ST values. This cannot be fully explained by the smaller area of the habitat region, as the estimation on the barrier is more accurate, with a comparable area. Most of the individuals trapped in a given habitat come from the same habitat (about 90%). The remaining individuals being sparse, this leads to higher relative variance in the proportions of individuals trapped in the habitats than in the other regions, which can explain the lower accuracy of the estimation of the diffusion parameter in the habitats. Based on these observations, we suspect that the estimation of a single coefficient in a homogeneous environment would most likely be reliable, even with low F ST values, and that placing the traps far from the release sites should lead to a better estimation of the coefficient in such case.

In addition to the Hardy-Weinberg equilibrium and independence of loci, an important assumption in our approach was that the allele frequencies were exactly known in the habitats. In practice, the frequencies are determined from previously sampled populations. The sample size is known to have an important effect on successful assignment of genotyped individuals [START_REF] Cornuet | New methods employing multilocus genotypes to select or exclude populations as origins of individuals[END_REF]. Reducing this size should lead to some uncertainty in the allele frequencies with an effect comparable to that of decreasing the F ST , i.e, lowering the discrimination power of the genotype data. A problematic case noted in [START_REF] Paetkau | Microsatellite analysis of population structure in canadian polar bears[END_REF] and [START_REF] Cornuet | New methods employing multilocus genotypes to select or exclude populations as origins of individuals[END_REF] while studying assignment methods is when some individuals carry an allele which has not been detected in the sample corresponding to their population of origin, leading to a null posterior probability that these individuals come from their true habitat. In such case, [START_REF] Paetkau | Microsatellite analysis of population structure in canadian polar bears[END_REF] suggest to add the genotype of these individuals to the population samples defining the allele frequencies in all of the habitats.

In our study, the allele frequencies are determined before dispersal, in individuals from the same generation as the trapped individuals. At each generation, the allele frequencies are modified, due to drift and gene flow among habitat patches. Thus, using allele frequencies sampled from previous generations could lead to an inaccuracy in the frequencies which depends on the gene flow and on the number of generations before the capture of the genotyped dispersers. To estimate the effect of the population flows on the allele frequencies after one generation, we computed the quantity C h τ /C τ , for τ = h, corresponding to the proportion of individuals captured in Ω h , whose habitat of origin is Ω h . The values clearly depend on the proximity of another habitat, with a ratio of 0.98 in Ω 1 , which is the farthest from the other habitats, and Comment citer ce document : Roques, L., Walker, E., Franck, P., Soubeyrand, S., Klein, E. (2016). Using genetic data to estimate diffusion rates in heterogeneous landscapes. Journal of Mathematical Biology, 73 (2), 397-422. DOI : 10.1007/s00285-015-0954-4 of 0.82 in Ω 3 and Ω 5 which are close to each other (the other ratios are 0.87, 0.85 and 0.94 in Ω 2 , Ω 5 and Ω 6 , respectively). Thus, the allele frequencies may remain stable after several generations if the habitats are sufficiently far from each other or become rapidly inaccurate in the opposite situation. A difficulty in estimating whether two habitats are far from each other is that the population flows are not known a priori. Our approach was based on an unbiased Brownian motion description of the individual trajectories, which leads to a Fokker-Planck reaction-diffusion equation. It can easily be extended to include a bias in a direction (v 1 , v 2 )(t, x) modelling attractiveness or repulsiveness of some elements of the landscape, by adding a term -∂(u v 1 )/∂x 1 -∂(u v 2 )/∂x 2 in the reaction-diffusion equation (2.2). The method proposed here also readily applies to other linear dispersal terms, such as the Fickian diffusion ∇ • (D(x)∇u) which is nevertheless more adapted to describe electric and thermic conductivity [START_REF] Roques | Modelling the impact of an invasive insect via reaction-diffusion[END_REF][START_REF] Roques | Modèles de réaction-diffusion pour l'écologie spatiale[END_REF]). More general Lévy processes than Brownian diffusion, corresponding to movements with large jumps, could have been considered as well; in such cases, the Laplace diffusion operator ∆ would be replaced by a fractional Laplace operator [START_REF] Valdinoci | From the long jump random walk to the fractional laplacian[END_REF]. Integral kernel-based dispersal terms which can account for long distance dispersal events could also be considered in place of the diffusion approach [START_REF] Kot | Dispersal data and the spread of invading organisms[END_REF].

The purpose of our study was to propose a rigourous method for the estimation of dispersal parameters in stochastic differential equation-based models of individual movement. Here, the method performance analysis assumed the same type of model for simulations and for inference. In practice, the real dispersal process may differ from the assumptions of the model. Future work should focus on the robustness of the method when the model assumptions are violated.

Another possible extension of our study is to include the estimation of the relative predispersal densities in the different habitats of origin and/or the areas of these habitats. It is not straightforward however that these problems are identifiable.

As suggested by Remark 3, pre-dispersal subpopulations may not coincide with the habitats. The location of the subpopulations may also not be known a priori. A challenging extension of our approach would consist in clustering the subpopulations and inferring the allelic frequencies (following a Structure-like approach, see [START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF] together with the estimation of D(x) in a full Bayesian approach, based on pre-dispersal genotype data in the habitats and genotype data of trapped dispersers.

Comment citer ce document : Roques, L., Walker, E., Franck, P., Soubeyrand, S., Klein, E. (2016). Using genetic data to estimate diffusion rates in heterogeneous landscapes. Journal of Mathematical Biology, 73 (2), 397-422. DOI : 10.1007/s00285-015-0954-4 and w 0 (x) = 0. Similarly, the cumulated population density of individuals coming from Ω h is: (6.32) and w h 0 (x) = 0.

w h t (x) = t 0 u h (
Conflict of Interest: The authors declare that they have no conflict of interest.

Fig. 1

 1 Fig.1Study site Ω = (0, 1) × (0, 1) with six habitats Ω h of centers x h , for h = 1, . . . , 6. The grey region Q corresponds to a barrier to dispersal. Blue crosses correspond to the positions of the first 10 traps θτ , τ = 1, . . . , 10; red crosses correspond to the positions of a supplementary set of ten additional traps θτ , τ = 11, . . . , 20.

Remark 4 (Fig. 2

 42 Fig. 2 Probability w h ∞ (x)/w∞(x) that an individual trapped at the position x ∈ Ω comes from the habitat Ω h ; see formula (4.17). White regions indicate that the probability is smaller than 0.005. The value of w h ∞ (x)/w∞(x) has been computed at t = 20, with the true parameter values (see Appendix D).

Fig. 3

 3 Fig.3Five sections of the 3-dimensional log-likelihood function ln(L( d1 , d2 , d3 )). The color scale corresponds to the gap between the log-likelihood and the maximum log-likelihood, where the maximum is taken over the parameter region ( d1 , d2 , d3 ) ∈ (d 1 -5, d 1 + 5) × (d 2 -5, d 2 + 5) × (d 3 -5, d 3 + 5).

Fig. 4

 4 Fig. 4 Effect of the F ST on the quality of the estimator D(x) in each of the three regions: D1 in the matrix (in white), D2 in the habitats (in green) and D3 in the barrier (in grey). The blue dashed lines correspond to the true values of D in the three regions: D 1 = 0.01, D 2 = 0.005 and D 3 = 0.001.

Fig. 6

 6 Fig. 6 Effect of number of traps and of the total number of genotyped individuals on the quality of the estimator D(x) in each of the three regions: D1 in the matrix (in white), D2 in the habitats (in green) and D3 in the barrier (in grey). The blue dashed lines correspond to the true values of D in the three regions: D 1 = 0.01, D 2 = 0.005 and D 3 = 0.001.

Table 1

 1 Summary of the notations used in the main text.

	Notation	Explanation
	u(t, x)	density of dispersers
	u h (t, x)	density of dispersers coming from habitat Ω h
	u 0 (x)	pre-dispersal density
	u h 0 (x)	pre-dispersal density of individuals coming from habitat Ω h
	α	pre-dispersal density in the habitats
	w∞(x)	cumulated density of dispersers
	w h ∞ (x)	cumulated density of dispersers coming from habitat Ω h
	βτ	capture rate in trap θτ
	Cτ	expected number of individuals captured in trap θτ
	C h τ	expected number of individuals coming from habitat Ω h
		captured in trap θτ
	D(x) and D 1 , D 2 , D 3	diffusion parameters (mobility)
	ν	life expectancy of the dispersers
	Ω	study site
	Ω h , h = 1, . . . , H	habitats (subsets of Ω)
	θτ , τ = 1, . . . , J	traps (subsets of Ω)
	x h	position of the center of the habitat Ω h
	xτ	position of the center of the trap θτ
	G	number of individuals genotyped in each trap
	λ = 1, . . . , Λ	index for the loci
	a = 1, . . . , A λ	index for the alleles
	(a 1 , a 2 )	couple of alleles at a given locus
	p hλa	frequency of allele a of locus λ in habitat Ω h
	F hλ	allele frequencies at locus λ in habitat Ω h
	G iτ	genotype of the i th genotyped individual in trap θτ
	M	measurement set consisting of all the genotypes in all traps

  s, x) ds, for all t > 0, x ∈ Ω.

				(6.31)
	This function satisfies:
	∂w h t ∂t	= ∆(D(x) w h t ) -	w h t ν	+ u h 0 (x), t > 0, x ∈ Ω,

The research leading to these results has received funding from the French Agence Nationale pour la Recherche, within the ANR-12-AGRO-0006 PEERLESS, ANR-13-ADAP-0006 MECC and ANR-14-CE25-0013 NONLO-CAL projects and from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n.321186 -ReaDi -Reaction-Diffusion Equations, Propagation and Modelling.

Appendix A: gradual release of the pre-dispersal populations

The equation (2.2) describes a simultaneous release of all the individuals at t = 0. To account for a possible gradual release of the individuals, the equation (2.2) can be replaced by: (6.26) where the term u 0 (x) f (t) describes the release of the individuals; u 0 (x) still corresponds to the pre-dispersal density and the function f (t) is the release rate. It can be described by any nonnegative function or distribution with integral 1 and with support in [0, T ], T corresponding to the end of the release period. In this framework, the density of dispersers coming from habitat Ω h satisfies the equation:

where u h 0 is still given by (2.7).

Appendix B: precise shape of the diffusion terms

In our numerical computations, we took

Version postprint

Comment citer ce document : Roques, L., Walker, E., Franck, P., Soubeyrand, S., Klein, E. (2016) for the function µ defined by (see Fig. 7):

Appendix C: computation of the F ST

The index F ST is used as a measure of genetic differentiation among the subpopulations. It was computed as follows: we set

where Λ is the number of loci, A, the number of alleles per locus whose frequency is measured and H the number of subpopulations, and .28) This formula corresponds to Nei's G ST for a single locus [START_REF] Nei | Analysis of gene diversity in subdivided populations[END_REF], with numerator and denominator averaged over the Λ loci. In our computations, all the subpopulations had the same size; in other situations, the weight 1/H in the above formulas for H S and H T should be replaced by the relative sizes of the subpopulations.

Appendix D: numerical computation of the cumulated population densities

In order to compute the cumulated densities w ∞ (x) and w h ∞ (x), we used the time-dependent partial differential equation solver Comsol Multiphysics c applied to the evolution equations (6.30) and (6.32) below at large time (t = 20), with default parameter values (finite element method with second order basis elements) and a triangular mesh adapted to the geometry of our landscape and made of 5296 elements.

We defined the cumulated population density at intermediate times t and position x by: w t (x) = t 0 u(s, x) ds, for all t > 0, x ∈ Ω. (6.29) Integrating (2.2) between 0 and t > 0 we note that w t (x) satisfies the following equation:

∂w t ∂t = ∆(D(x) w t ) -w t ν + u 0 (x), t > 0, x ∈ Ω, (6.30)


Sébastien Rimbert
email: sebastien.rimbert@inria.fr

Oleksii Avilov

Laurent Bougrain

DISCRETE MOTOR IMAGERIES CAN BE USED TO ALLOW A FASTER DETECTION

Motor imagery (MI) modifies the neural activity within the primary sensorimotor areas of the cortex and can be measured through the analysis of electroencephalographic (EEG) recordings. It is particularly interesting for Brain-Computer Interface (BCI) applications. In most MI-based BCI experimental paradigms, subjects realize continuous motor imagery (CMI), i.e. a repetitive and prolonged intention of movement, for a few seconds. The system detects the movement based on the event-related desynchronization and the event-related synchronization features in electroencephalographic signal. Currently, improving efficiency such as detecting faster a motor imagery is an important issue in BCI to avoid fatigue and boredom. The purpose of this study is to show the difference, in term of classification, between a discrete motor imagery, i.e. a single short MI, and a CMI. The results of experiments involving 16 healthy subjects show that a BCI based on DMI is as effective as a BCI based on CMI and could be used to allow a faster detection.

INTRODUCTION

Motor imagery (MI) is the ability to imagine performing a movement without executing it [START_REF] Avanzino | Motor imagery influences the execution of repetitive finger opposition movements[END_REF]. According to Jeannerod [START_REF] Jeannerod | Mental imagery in the motor context[END_REF], MI represents the result of conscious access to the content of the intention of a movement, which is usually performed unconsciously during movement preparation [START_REF] Lotze | Motor imagery[END_REF]. MI has two different components, namely the visual-motor imagery and the kinesthetic motor imagery (KMI) [START_REF] Neuper | Imagery of motor actions:Differential effects of kinesthetic and visual motor mode of imagery in single-trial EEG[END_REF]. KMI generates an event-related desynchronization (ERD) and an event-related synchronization (ERS) in the contralateral sensorimotor area, which is similar to the one observed during the preparation of a real movement (RM) [START_REF] Pfurtscheller | Eventrelated eeg/meg synchronization and desynchronization: basic principles[END_REF]. More precisely, compared to a resting state taken before a motor imagery, several power modulations are observed in in the alpha [START_REF] Jonathan Wolpaw | Brain-Computer Interfaces: Principles and Practice[END_REF][START_REF] Butler | Mental practice with motor imagery: evidence for motor recovery and cortical reorganization after stroke[END_REF][START_REF] Blokland | Decoding motor responses from the eeg during altered states of consciousness induced by propofol[END_REF][START_REF] Ahn | Performance variation in motor imagery brain-computer interface: a brief review[END_REF][START_REF] Thomas | Investigating brief motor imagery for an erd/ers based bci[END_REF] and in the beta [START_REF] Pastotter | Oscillatory correlates of controlled speed-accuracy tradeoff in a response-conflict task[END_REF][START_REF] Parkes | Combining eeg and fmri to investigate the postmovement beta rebound[END_REF][START_REF] Schnitzler | Involvement of primary motor cortex in motor imagery: a neuromagnetic study[END_REF][START_REF] Erbil | Changes in the alpha and beta amplitudes of the central eeg during the onset, continuation, and offset of long-duration repetitive hand movements[END_REF][START_REF] Kilavik | The ups and downs of beta oscillations in sensorimotor cortex[END_REF][START_REF] Pfurtscheller | Contrasting behavior of beta event-related synchronization and somatosensory evoked potential after median nerve stimulation during finger manipulation in man[END_REF][START_REF] Salenius | Modulation of human cortical rolandic rhythms during natural sensorimotor tasks[END_REF][START_REF] Hohne | Motor imagery for severly motor-impaired patients: Evidence for brain-computer interfacing as superior control solution[END_REF] bands of the electroencephalographic signal measured over the sensorimotor area corresponding to the body part involved in the motor imagery. Firstly there is a gradual power decrease in the alpha and in the beta bands, called ERD. Secondly, a low power level is maintained during the movement. Finally, from 300 to 500 milliseconds after the end of the motor imagery, there is a power increase called ERS or post-movement beta rebound with a duration of about one second. Although several studies showed an activity uniquely in the contralateral area [START_REF] Pfurtscheller | Motor imagery and direct brain-computer communication[END_REF], other studies showed that ERD and ERS are also in the ipsilateral area [START_REF] Fok | An eeg-based brain computer interface for rehabilitation and restoration of hand control following stroke using ipsilateral cortical physiology[END_REF].

Emergence of ERD and ERS patterns during and after a MI has been intensively studied in the Brain-Computer Interface (BCI) domain [START_REF] Jonathan Wolpaw | Brain-Computer Interfaces: Principles and Practice[END_REF] in order to define detectable commands for the system. Hence, a better understanding of these processes could allow for the design of better interfaces between the brain and a computer system. Additionally, they could also play a major role where MI are involved such as rehabilitation for stroke patients [START_REF] Butler | Mental practice with motor imagery: evidence for motor recovery and cortical reorganization after stroke[END_REF], monitoring consciousness during general anesthesia [START_REF] Blokland | Decoding motor responses from the eeg during altered states of consciousness induced by propofol[END_REF] or the recovery of the motor capacity after neurological damage. For example MI training is a promising approach in facilitating paretic limb recovery.

Currently, most of the paradigms based on MIs require the subject to perform the imagined movement several times for a predefined duration of a few seconds. In this study, such a task is commonly referred to as a continuous motor imagery (CMI). However, first the duration of the experiment is long. Second a succession of flexions and extensions generate an overlapping of ERD and ERS patterns making the signal less detectable. In fact, one simple short MI, referred in this article as a discrete motor imagery (DMI), could be more useful for two reasons. Firstly, a DMI could be used to combat fatigue and boredom for BCI users improving ERD and ERS production [START_REF] Ahn | Performance variation in motor imagery brain-computer interface: a brief review[END_REF]. Secondly, the ERD and ERS generated by the DMI could be detectable at a higher quality and more rapidly compared to a CMI. This was found in a previous study that established a relationship between the duration of the MI and the quality of the ERS extracted [START_REF] Thomas | Investigating brief motor imagery for an erd/ers based bci[END_REF]. It also showed that a brief MI (i.e. a 2-seconds MI) could be more efficient then a sustained MI. Our main hypothe-sis is that a DMI generates robust ERD and ERS patterns which could be detectable by a BCI system. Results indicate that a DMI produces a robust ERS and is as detectable as a CMI.

MATERIALS AND METHODS

Participants: 16 right-handed healthy volunteer subjects took part in this experiment (9 men and 7 women, from 19 to 43 years old). They had no medical history which could have influenced the task such as diabetes, peripheral neuropathology, renal insufficiency, anti-depressant treatment or motor problem. All subjects gave their agreement and signed an information consent form approved by the ethical INRIA committee before participating. This experiment follows the statements of the WMA declaration of Helsinki on ethical principles for medical research involving human subjects [START_REF] Medical | World medical association declaration of helsinki: ethical principles for medical research involving human subjects[END_REF].

Real movement: The first task consisted of an isometric flexion of the right index finger on a computer mouse. A low frequency beep indicated when the subject had to execute the task. The right-click is recorded as a trigger and has allowed to know exactly when the participant executes the RM.

Discrete imagined movement: The second task was a DMI of the previous real movement. A low frequency beep indicated when the subject has to execute the task.

Continuous imagined movement: The third task was a CMI during four seconds of the real movement of the first task. More precisely, the subject imagined several (around four) flexions and extensions of the right index finger. This way, the DMI differed from the CMI by the repetition of the imagined movement. For this task, two beeps, respectively with low and high frequencies, separated by a four second delay, indicated the beginning and the end of the CMI.

Protocol: Each of the three tasks introduced in section corresponds to a session. The subjects completed three sessions during the same day. All sessions were split into several runs. Breaks of a few minutes were planned between sessions and between runs to avoid fatigue. Before each session, the task was described, and the subject practiced the tasks. At the beginning of each run, the subject was told to relax for 30 seconds. Condition 1 corresponded to RMs was split into 2 runs of 50 trials. Conditions 2 and 3 corresponded to discrete and continuous imagined movements, respectively, was split into 4 runs of 25 trials. Thus, 100 trials were performed by subjects for each task. Each experiment began with condition 1 as session 1. Conditions 2 and 3 were randomized to avoid possible bias cause by fatigue, gel drying or another confounding factor. For conditions 1 and 2, the timing scheme of a trial was the same: one low frequency beep indicated the start followed by a rest period of 12 seconds. For condition 3, a low frequency beep indicated the start of the MI to do during 4 seconds, followed by a rest period of 8 seconds. The end of the MI is announced by a high frequency beep (Fig. 1). Behavioral data: A custom-written scenario for OpenViBE [START_REF] Renard | Openvibe: An opensource software platform to design, test and use brain-computer interfaces in real and virtual environments[END_REF] was designed to automate the generation of beeps, and to record triggers and EEG signals. The triggers corresponding to the right-click allowed us to detect potential behavioral errors. All non realized movement were removed from the analysis. For all three tasks, we used a fixed preparatory period duration in which the subjects could anticipate the GO signal.

EEG data: EEG signals were recorded through the OpenViBE platform with a commercial REFA amplifier developed by TMS International. The EEG cap was fitted with 9 passive electrodes re-referenced with respect to the common average reference across all channels over the extended international 10-20 system positions to cover the primary sensorimotor cortex. The selected electrodes are FC3, C3, CP3, FCz, Fz, CPz, FC4, C4, CP4. No additional filtering was used during the recording. Skinelectrode impedances were kept below 5 kΩ. Incorrect trials were removed from the analyses.

ERD/ERS patterns: To evaluate more precisely the modulation produced by the tasks, we computed the ERD/ERS% using the "band power method" [START_REF] Pfurtscheller | Eventrelated eeg/meg synchronization and desynchronization: basic principles[END_REF] with a matlab code. First, the EEG signal is filtered between 8-30 Hz (Alpha+Beta) for all subjects using a 4thorder Butterworth band-pass filter. Then, the signal is squared for each trial and averaged over trials. Then it is smoothed using a 250-millisecond sliding window with a 100 ms shifting step. We have chosen a specific sliding window because the nature of the real and imagined movement, as well as the components ERD/ERS that underline them, require a short window. Finally, the averaged power computed for each window was subtracted and then divided by the averaged power of a baseline corresponding to 2 seconds before each trial. Finally, the averaged power computed for each window was subtracted and then divided by the averaged power of a baseline corresponding 2 seconds before each trial. This transformation was multiplied by 100 to obtain percentages. This process can be summarized by the following equation:

ERD/ERS% = x 2 -BL 2 BL 2 × 100 , (1)
where x 2 is the average of the squared signal over all trials and samples of the studied window, BL 2 is the mean of a baseline segment taken at the beginning of the corresponding trial, and ERD/ERS% is the percentage of the oscillatory power estimated for each step of the sliding window. It is done for all channels separately. ERD and ERS are difficult to observe from the EEG signal. Indeed, an EEG signal expresses the combination of activities from several neuronal sources. One of the most effective and accurate techniques used to extract events is the average technique [START_REF] Quiroga | Single-trial event-related potentials with wavelet denoising[END_REF]. We decided to use this technique to represent the modulation of power of the Al-pha+Beta rhythms for three tasks (Real Movement, Discrete Motor Imagery and Continuous Motor Imagery).

Common Spatial Pattern: We used the algorithm called Common Spatial Pattern (CSP) to extract motor imagery features from EEG signals; this generated a series of spatial filters that were applied to decompose multi-dimensional data into a set of uncorrelated components [START_REF] Blankertz | Optimizing spatial filters for robust EEG single-trial analysis[END_REF]. These filters aim to extract elements that simultaneously maximize the variance of one class, while minimizing the variance of the other one.

Feature extraction and linear discriminant analysis: We trained a linear discriminant classifier to distinguish the features of motor imageries from the ones of a resting state. We applied the common spatial pattern algorithm to obtain 3 pairs of linear combinations from the 8-30 Hz filtered EEG signals. Then for each linear combinations we computed the logarithm of the variance for a studied window. We considered a 2 seconds window taken 3 seconds before the GO signal for the resting state. The 1 second window before the GO signal is not taken into consideration because the beep can generate an audio ERP and the subject usually prepares the movement in advance. So it not really a resting state. The features of a DMI is computed from 0.2 to 1 second after the GO signal (Fig. 5, Tab. 1). The features of a CMI is computed from 0.2 to 3 seconds after the GO signal (Fig. 5, Tab. 1). Table 1: Grand average accuracies obtained by linear discriminant analyses for the 3 conditions (RM, CMI and DMI) and 3 frequency bands (alpha, beta and alpha + beta).

RESULTS

To prove the usability of a DMI in BCI-domain, firstly we computed ERD and ERS patterns to study the relative power (8-30 Hz) for the electrode C3. Secondly, we verified the detectability of a DMI in calculating classification rate.

ERD and ERS modulation: To verify if a DMI generates ERD and ERS patterns which could be detectable by a CMI, we studied the relative power (8-30 Hz) for the electrode C3. Electrode C3 is suitable for monitoring right hand motor activity. A grand average was calculated over the 16 subjects. We used a Friedman's test to analyze whether ERS were significantly and respectively different during the three conditions. Due to eyes-closed experiment, the alpha band is disturbed (confirmed by the time-frequency analysis) and not considered for this study. Consequently values corresponding to the desynchronization during the real and imagined movements will appear smaller. Moreover a visual inspection of time-frequency analysis shown modulations of alpha + beta power between 8-30 Hz for all the 16 subjects (Fig. 3).

ERD and ERS modulation during/after a real movement:

The ERD/ERS% averages (Fig. 4) indicate that one second after the cue, the power in the 8-30 Hz band increases by around 60%, reaches its maximum and returns to the baseline 7 seconds after. The evolution from ERD to ERS is rapid (less than one second) and should be linked to the type of movement realized by the subjects. Interestingly, each subject (except Subject 9 and Subject 16) has the same ERD/ERS% profile (i.e. a strong rebound) after the real movement (Fig. 3).

ERD and ERS modulation during/after a discrete motor imagery: The ERS post-DMI reaches 18% which is weaker compare to the other tasks (Fig. 4). Some subjects (S1, S2, S5, S6, S10) have a stronger robust ERS produced by DMI while others have no beta rebound. Some subjects (S9, S10, S15) have a strong ERD after the task (Fig. 3). This variability between subjects could explain the weakness of the ERS post-DMI. The presence of ERD and ERS during/after a DMI suggest that a DMI could be used in BCI-domain.

ERD and ERS modulation during/after a continuous motor imagery: During the CMI, the subjects imagined several movements in a time window of 4 seconds. Fig. 4 shows a global decrease of activity during the CMI and stronger modulation in 8-30 Hz after the CMI. The results of the grand average shows a low desynchronization during this time window. It is interesting to note that some subjects (S7, S9, S11) have no desynchronization during the CMI task and could have a negative effect on the classification phase. Other subjects (S2, S15) have a different profile which shows that a first ERS is reached one second after the beginning of the CMI, then the power increases and decreases again, being modulated during 3 seconds. Indeed, this global ERD can be considered as the concatenation of several ERDs and ERSs due to the realization of several MIs. The variability between subjects during this period could also have a negative effect on the classification rate.

Detection results: Discrete motor imageries generate robust ERD and ERS (see previous section). In this section, we will study if they are detectable enough to have a faster detection in BCI than using continuous motor imageries. For each subject, 4 runs of 25 trials were available. The process of the cross validation consisted in using trials of 3 runs for train classifier and 1 run for test it. Four permutations were possible and we averaged accuracies obtain by the 4 classifiers on their testing run for a better evaluation. This method of cross validation was chosen because of its proximity to online condition. Figure 5 presents the accuracy for each subject and the mean accuracy. The mean accuracy for RM, DMI and CMI are respectively 78,4%, 71,9% and 71,4%. The detection of real movement is easier than the one of motor imageries. The difference between the two motor imageries is not statistically significant at a level of 5%. The precision of the ME is 0,78 (22% of false positive) and the recall is 0.79. The precision of the DMI is 0.71 (29% of false positive) and the recall is 0.74. The precision of the CMI is 0,73 (27% of false positive) and the recall is 0,68. Interestingly, some subjects (S5, S6, S9, S16) have a better detection for the DMI task. In Tab. 1, we computed the grand average of accuracies for the three conditions (RM, DMI and CMI) for three frequency bands (alpha, beta and alpha + beta). It appears that the 8-30 Hz frequency band increases the classification rate. Furthermore, the comparison of the classification rate between a DMI and a CMI on the same period (0.2s-1s) shows an equivalence.

DISCUSSION

The subjects carried out voluntary movements, DMI and CMI of an isometric flexion of the right hand index finger. Results show that the power in the 8-30 Hz band is modulated during the three tasks. The comparison between ERSs suggests that subjects on average have a stronger ERS during a CMI than a DMI. However, this is not the case for all subjects. Furthermore, the detection rate for a DMI is as effective as for a CMI.

EEG system: It is well established that a large number of electrodes allows having a good estimation of the global average potential of the whole head [START_REF] Dien | Issues in the application of the average reference: review, critiques, and recommendations[END_REF]. Although we focused on specific electrodes, our results were similar by using method of the derivation, which corresponded with the literature. We chose to study C3 without derivation because we are interested in designing a minimal system to detect ERD and ERS during general anesthesia conditions.

ERD/ERS modulation during real movements: The results are coherent with previous studies describing ERD/ERS% modulations during motor actions. The weakness of the ERD can be linked to the instruction that was focused more on the precision than the speed of the movement [START_REF] Pastotter | Oscillatory correlates of controlled speed-accuracy tradeoff in a response-conflict task[END_REF]. However, although some subjects were making efforts to do a voluntary movement, we must consider that an isometric flexion movement on a mouse is a movement setting in the subject's memory. This can have an impact on the low ERD amplitude. We also showed that the rebound starts before the click. Since a mouse click, is a really fast movement, we expect that the beta rebound will appear fast as well [START_REF] Parkes | Combining eeg and fmri to investigate the postmovement beta rebound[END_REF].

ERS modulation during motor imageries: The results show that the ERS is lower after a DMI or a CMI than after a real movement, which has been already been demonstrated previously [START_REF] Schnitzler | Involvement of primary motor cortex in motor imagery: a neuromagnetic study[END_REF]. However, the novelty is the beta rebound is stronger on average after a CMI than DMI for a few subjects.

ERD modulation during continuous motor imagery: When the subjects performed the CMI, the ERD was highly variable during the first 4 seconds. For some subjects, our hypothesis is there are some intern-ERD and intern-ERS into this period. The difficulty is that the CMI involves several MI, that are not synchronized across trials, unlike the DMI which starts and ends at roughly the same time for each trial, due to the cue. Normally, for continuous real movement, the ERD was sustained during the execution of this movement [START_REF] Erbil | Changes in the alpha and beta amplitudes of the central eeg during the onset, continuation, and offset of long-duration repetitive hand movements[END_REF]. However, in our data it is possible to detect several ERDs during the 4 seconds of CMI where the subject performed 4 MIs. This assumes that the ERD and ERS components overlap in time when we perform a CMI. Several studies already illustrate the concept of overlap of various functional processes constituting the beta components during RMs [START_REF] Kilavik | The ups and downs of beta oscillations in sensorimotor cortex[END_REF]. Moreover, the beta rebound generated by a median nerve stimulation is reduced when the stimulation is made during different types of real or imagined hand movements [START_REF] Pfurtscheller | Contrasting behavior of beta event-related synchronization and somatosensory evoked potential after median nerve stimulation during finger manipulation in man[END_REF], [START_REF] Salenius | Modulation of human cortical rolandic rhythms during natural sensorimotor tasks[END_REF]. However, even if the components cancel each other out in the signal, it does not mean that the operation of the underlying processes are similarly affected. This interpretation assumes implicitly that the components are combining each other, which means that the temporal superposition of an ERD and an ERS would result in an intermediate amplitude signal. This could explain why the ERD during a CMI could be less detectable and more varied than the ERD during a DMI. To validate this hypothesis, we plan to design a new study to explore how two fast-successive movements (or MIs) can affect the signal in the 8-30 Hz frequency band.

Detection rate: We showed that the detection of a real movement was easier than discrete and continuous motor imageries but we could discuss about the weakness of the classification rate for the three task. Usually, a real movement has often a high classification rate and it is not the case in this study. However, it is important to remind that the subject performed real movement, DMI and CMI of an isometric flexion of the right hand index finger. The rapidity and the precision of the three tasks could be linked with the low classification rate. One limitation of this study is that all trials of a type had the same length and was not randomized within a block.

Establishing a link between the ERD/ERS users profiles and the detection rate: Our study shows results which could allow to understand more differences, in term of ERD and ERS, between subjects. Indeed, we showed that for a same task (RM, DMI and CMI), for some subjects a strong ERS appeared whereas for some others, no ERS appeared. We observed the same phenomenon for ERD. It could be interesting to establish a link between the particular ERD/ERS users profiles and the detection rate. The importance of BCI users profiles, especially for patients with severe motor impairments has already been established by other studies [START_REF] Hohne | Motor imagery for severly motor-impaired patients: Evidence for brain-computer interfacing as superior control solution[END_REF]. This is why, we expect designing an adaptive BCI based on the specific motor activity of the motor cortex. More subjects are necessary to precise this BCI user profile.

CONCLUSION

This article examined the modulation of power in EEG during a real movement, a discrete motor imagery (DMI) and a continuous motor imagery (CMI). We showed that during a real voluntary movement corresponding to an isometric flexion of the right hand index finger a low ERD appeared, and was followed by a rapid and powerful ERS. Subsequently, we showed that the ERD and ERS components were still modulated by both a DMI and a CMI. The ERS is present in both cases and shows that a DMI could be used in BCI domain. The classification results show no any difference between a CMI and a DMI and confirm that a DMI could have a future impact in BCI-domain to save time and avoid fatigue.

Figure 1 :

 1 Figure 1: Timing schemes of a trial for each task: Real Movement (RM, top); Discrete Motor Imagery (DMI, middle); Continuous Motor Imagery (CMI, bottom). The DMI and CMI sessions are randomized.

Figure 2 :

 2 Figure 2: Schematic representation of the experiment. A low frequency beep indicates the start of the (real or imagined) movement. A high frequency beep indicates the end of the continuous imagined movement. Depending on the task, the subject presses or imagines pressing the button of the mouse.

Figure 3 :

 3 Figure 3: Grand average ERD/ERS% curves (in black, Average) estimated for the RM, the DMI and the CMI within the alpha + beta band (8-30 Hz) for electrode C 3 . The average for each subject is also presented. A first beep indicated the start of the (real or imagined) movement. A second beep indicated the end of the continuous imagined movement.

Figure 4 :

 4 Figure 4: Grand average ERD/ERS% curves estimated for the RM (blue), the DMI (green) and the CMI (orange)

Figure 5 :

 5 Figure 5: Accuracies obtained by linear discriminant analyses for the 3 conditions (RM, CMI and DMI). The features of a RM and DMI are computed from 0.2 to 1 second after the GO signal. The features of a CMI are computed from 0.2 to 3 seconds after the GO signal.