Joël Chadoeuf

Alexandre Millon
email: alexandre.millon@imbe.fr

Jean-Luc Bourrioux

Thierry Printemps

Benoît Van Hecke

Vincent Lecoustre

Vincent Bretagnolle

Modelling unbiased dispersal kernels over continuous space by accounting for spatial heterogeneity in marking and observation efforts

Keywords: Montagu's harrier, Circus pygargus, citizen science, imperfect detection, bias correction Estimating unbiased dispersal kernels

Methods in Ecology and Evolution

unbiased dispersal kernels, a key component for modelling population dynamics and species distribution in a context of environmental change. Currently, our method assumes homogeneity in both habitat and dispersal behaviour across individuals. We discuss however how to relax these hypotheses to further investigate the effect of e.g. local conspecific density 50 or habitat quality on dispersal propensity.

1. Although a key demographic trait determining the spatial dynamics of wild populations, dispersal is notoriously difficult to estimate in the field. Indeed, dispersal distances obtained from the monitoring of marked individuals typically lead to biased estimations of dispersal kernels as a consequence of i) restricted spatial scale of the study areas compared to species potential dispersal and ii) heterogeneity in marking and observation efforts and therefore in detection probability across space.

2. Here we propose a novel method to circumvent these issues that does not require data on observation effort per se, to correct for the variability in detection of marked individuals across space. Observed dispersal events were weighted by the distribution of departure points and an eroded spatial window approach was applied so as to deal with border effect. We conducted a set of simulations which indicated that our method was successful in correcting the effect of spatially heterogeneous detectability and produce unbiased dispersal kernels.

3. We applied this method to a real dataset on Montagu's harrier (>5000 chicks tagged), providing ca. 6000 resightings collected in entire France by a network of 1200 volunteers within a citizen-science program. The median dispersal distance observed was 32 km (range: 0.1-627 km). Once corrected for spatial heterogeneity in marking and observation efforts and border effect, the modelled dispersal kernel indicated a median dispersal distance of 78-123 km depending on the spatial scale considered (constrained within French borders or not, respectively).

Synthesis and applications:

The current rise of citizen-science programs is likely to stretch our estimate of the ecologically-relevant spatial scale at which dispersal takes place for many taxa. Our method is particularly suited for such large scale data that typically suffer from high spatial heterogeneity in marking and observation efforts and offers the possibility to derive

INTRODUCTION

Dispersal is defined either as the movement of an individual from its birthplace to the location where it will reproduce for the first time (natal dispersal), or the movement between successive reproduction events (adult dispersal). Dispersal is a pivotal process underpinning spatial population dynamics (review in [START_REF] Clobert | Dispersal -Causes, Consequences and Mechanisms of Dispersal at the Individual[END_REF]. Indeed, movements of individuals, both within and among populations, sustain gene flow across space and time, thus driving genetic diversity as well as adaptation to local conditions [START_REF] Cain | Long-distance seed dispersal in plant populations[END_REF][START_REF] Postma | Gene flow maintains a large genetic difference in clutch size at a small spatial scale[END_REF]. Dispersal is a spatial demographic trait that affects extinction and colonisation rates of habitat patches in a metapopulation context [START_REF] Hanski | Metapopulation Ecology[END_REF][START_REF] Sutherland | Multi-scale processes in metapopulations: contributions of stage structure, rescue effect, and correlated extinctions[END_REF], and also the size of more continuously spread populations [START_REF] Doncaster | Balanced dispersal between spatially varying local populations: An alternative to the source-sink model[END_REF][START_REF] Schaub | Strong contribution of immigration to local population regulation: Evidence from a migratory passerine[END_REF][START_REF] Lieury | Geographically isolated but demographically connected: Immigration supports efficient conservation actions in the recovery of a range-margin population of the Bonelli's eagle in France[END_REF].

Ignoring spatial demographic processes when investigating population dynamics may indeed lead to spurious interpretations regarding e.g. the conservation status of a population [START_REF] Schaub | Massive immigration balances high anthropogenic mortality in a stable eagle owl population: Lessons for conservation[END_REF], the efficiency of management actions [START_REF] Lieury | Compensatory immigration challenges predator control: An experimental evidence-based approach improves management[END_REF], the estimation of fitness [START_REF] Dhondt | Summer dispersal and survival of juveniles Great Tits in southern Sweden[END_REF][START_REF] Tinbergen | Biased estimates of fitness consequences of brood size manipulation through correlated effects on natal dispersal[END_REF][START_REF] Doligez | Estimating fitness consequences of dispersal: a road to 'knowwhere'? Non-random dispersal and the underestimation of dispersers' fitness[END_REF] and more generally the evolution of populations in response to environmental changes [START_REF] Kokko | From individual dispersal to species ranges: perspectives for a changing world[END_REF][START_REF] Ronce | How does it feel to be like a rolling stone? Ten questions about dispersal evolution[END_REF]). In addition, the way dispersal is accounted for greatly affects the results of species distribution models predicting range shifts in response to climate change [START_REF] Guisan | Predicting species distribution: offering more than simple habitat models[END_REF][START_REF] Bellard | Impacts of climate change on the future of biodiversity[END_REF].

However, estimating dispersal is notoriously difficult and has been the subject of intense methodological developments over the last two decades [START_REF] Van Noordwijk | On bias due to observer distribution in the analysis of data on natal dispersal in birds[END_REF][START_REF] Van Noordwijk | Measuring dispersal as distance-dependent recruitment rates: Testing the performance of DDRR on simulated data[END_REF][START_REF] Koenig | Detectability, philopatry, and the distribution of dispersal distances in vertebrates[END_REF][START_REF] Paradis | Modeling large-scale dispersal distances[END_REF][START_REF] Schaub | Estimating true instead of apparent survival using spatial Cormack-Jolly-Seber models[END_REF]. The classical way for investigating the dispersal process involves the longitudinal monitoring of marked individuals over space. A major issue of such studies, however, arises from the finite nature of the study area, as the closer a departure point is from the border, the higher the probability of the corresponding arrival point to fall outside the study area. This results in the so-called border effect, which produces an inevitable reduction in the detectability of longdistance dispersal [START_REF] Barrowclough | Sampling bias in dispersal studies based on finite area[END_REF]. Indeed, in most instances, the spatial scale of the study area is defined according to logistical constraints so as to ensure an exhaustive marking of offspring and maximising recapture/resighting probability of previously marked individuals. However, there is growing evidence for the dispersal ability of individuals to largely overcome the limits of 'standard' study areas [START_REF] Van Houtan | Dispersal of Amazonian birds in continuous and fragmented forest[END_REF]. In birds particularly, data collected from individuals fitted with VHF, Argos or GPS devices have revealed much larger dispersal distances than those previously estimated using passive marks [START_REF] Van Vuren | Survival of dispersing and philopatric yellowbellied marmots: what is the cost of dispersal?[END_REF][START_REF] Koenig | Detectability, philopatry, and the distribution of dispersal distances in vertebrates[END_REF][START_REF] Hénaux | Tracking large carnivore dispersal using isotopic clues in claws : an application to cougars across the Great Plains[END_REF].

Weight and economic costs of such electronic devices currently limit the number of species on which such equipment can be deployed and the number of equipped individuals, therefore limiting the scope for evaluating population-level dispersal kernels [START_REF] Murray | A critical review of the effects of marking on the biology of vertebrates[END_REF][START_REF] Wikelski | Going wild : what a global small-animal tracking system could do for experimental biologists[END_REF]. There is therefore a pressing need to improve our ability to derive unbiased dispersal kernels from 'classically marked' individuals and beyond the scale of a single study area.

The rise of citizen-based science offers a unique opportunity to take up this challenge by providing observations over a large spatial scale, typically nation-wide scale. In birds, large spatial scale monitoring programs dedicated to the study of dispersal have started to emerge, such as in seabirds, waders or raptors [START_REF] Etheridge | The Effects of Illegal Killing and Destruction of Nests by Humans on the Population Dynamics of the Hen Harrier Circus cyaneus in Scotland[END_REF][START_REF] Gill | The effects of disturbance on habitat use by black-tailed godwits Limosa limosa[END_REF][START_REF] Barlow | Estimating dispersal distributions at multiple scales: within-colony and among-colony dispersal rates, distances and directions in European Shags Phalacrocorax aristotelis[END_REF]). Yet, increasing the spatial scale of the study inevitably comes at the expense of high spatial heterogeneity in marking effort and in the probability of detecting returning individuals, representing the second major issue when studying dispersal. Indeed, over large spatial scales, the probability of detecting a marked animal is typically low and cannot be considered constant across space. Data on distribution of the observation effort, that could help estimating this bias, are however usually missing in a citizen-science context. Different statistical methods have been proposed for estimating the bias in the observation process when investigating dispersal. Multi-states models within the Capture-Mark-Recapture framework (CMR) provide a way to estimate dispersal probability among a set of discrete and finite number of study areas [START_REF] Spendelow | Estimating annual survival and movement rates of adults within a metapopulation of roseate terns[END_REF][START_REF] Schwarz | A general methodology for the analysis of capturerecapture experiments in open populations[END_REF]. This approach is appropriate for species with highly patchy and restricted distributions such as rare colonial seabirds. However the number of parameters exponentially increases with the number of sites included, and moreover, this framework provides movement probabilities, not dispersal kernel (e.g. [START_REF] Lagrange | Estimating dispersal among numerous sites using capture-recapture data[END_REF]. For more continuously distributed species, the method classically used consists in increasing the scale of the sampling area and weighting the number of observed dispersal events (i.e. recruitments of individuals previously marked) by the total number of individuals caught or recruited (i.e. individuals not previously marked) at different distance classes. Obtaining such information on unmarked individuals entails substantial additional costs and thus reduces the scope of this method otherwise efficient in providing unbiased dispersal kernels [START_REF] Baker | Correcting biased estimates of dispersal and survival due to limited study area : theory and an application using wrentits[END_REF][START_REF] Van Noordwijk | On bias due to observer distribution in the analysis of data on natal dispersal in birds[END_REF][START_REF] Winkler | The natal dispersal of tree swallows in a continuous mainland environment[END_REF][START_REF] Van Houtan | Dispersal of Amazonian birds in continuous and fragmented forest[END_REF][START_REF] Van Noordwijk | Measuring dispersal as distance-dependent recruitment rates: Testing the performance of DDRR on simulated data[END_REF] proposed the use of distance-dependent recruitment rates to test whether the dispersal rule exhibited by the species differs from a uniform law, when the observation process is incomplete across space.

However, this method does not provide an estimation of dispersal kernels (see Appendix A where we demonstrate that the form of the distance-dependent recruitment rates depends on the spatial distribution of the observation effort). A promising method recently proposed is the development of a spatial extension of the Cormack-Jolly-Seber CMR model that allows one to jointly estimate true instead of apparent survival, dispersal and observation processes [START_REF] Schaub | Estimating true instead of apparent survival using spatial Cormack-Jolly-Seber models[END_REF]. However, a critical assumption of this model (primarily dedicated to the estimation of true instead of apparent survival) is that recapture probability is homogeneous across the study area, a clearly unrealistic assumption in nation-wide citizenbased programs. Therefore to date, there is no satisfactory methods producing unbiased dispersal kernels for species widely distributed over continuous space.

Here we developed such a method to estimate unbiased dispersal kernels at any spatial scale (including large nation-wide scale) and for any species distribution type (patchy or continuous), relying on three successive steps: 1) we minimised the bias induced by the border effect by considering the set of resighted individuals whose departure point lies within an eroded spatial window; 2) we estimated dispersal kernel parameters using a contrast function (i.e. an estimation function whose expectation is minimum at the true parameter set; Dacunha-Castelle & Duflo 1983) which converges to a weighted likelihood, up to a multiplicative constant proportional to the average recapture rate. We accounted for spatial heterogeneity in observation effort by weighting each dispersal event by the inverse of the density of departure points (i.e. marking effort). By doing this, we avoided the direct modelling of observation effort, typically intractable when dealing with citizen-science data;

3) finally, we provided an unbiased dispersal kernel by estimating the residual bias in dispersal parameters using a simulation and extrapolation method [START_REF] Cook | Simulation-extrapolation estimation in parametric measurement error models[END_REF].

As a first approach, we assumed that habitat quality was homogeneous across space and that all individuals had a similar average dispersal behaviour (but see Discussion for the possibility to relax these hypotheses).

The objectives of this paper are twofold. After elaborating on the method, its assumptions and mathematical developments, we checked for its robustness (validity, properties) with a simulation study based on a heavy-tailed distribution of dispersal events (log-sech) and spatially heterogeneous marking and observation efforts. Thereafter, we applied this method to a dataset gathered by a large-scale citizen-based science program conducted in France and designed to estimate the dispersal kernel of the Montagu's harrier Circus pygargus, a widespread species breeding in cereal intensive farming landscapes of Western Europe and expected to disperse over several hundreds of kilometres.

METHODS

DISPERSAL KERNEL ESTIMATION

Modelling dispersal data accounting for observation effort

Let consider ܫ individuals located at their departure points (e.g. birth nests) ܺ = ݔ( ଵ , … , ݔ ூ ) in an observation area ܹ ∈ ܴ ଶ . We suppose that ܺ follows a spatial Poisson process described by its intensity .)ݔ(ߣ The individual ݅ located at ݔ undergoes a single dispersal event to the arrival point ݕ with respect to a dispersal kernel

݃ ఏ బ such that )ݔ|ݕ( = ݃ ఏ బ ݕ( -,)ݔ where ߠ
is the set of parameters of the dispersal model ݃. The dispersal kernel is thus supposed to be stationary over W. The individual ݅ is observed at its arrival point when ݕ ∈ ܹ and when an observer detects it (either by direct capture or observation at distance) with a probability ݕ(߱ ), or the detection probability linked to observation effort. We define

)݅( = 1 if the individual ݅ is 'captured' by an observer in ܹ, )݅( = 0 if not. Let J = ∑ )݅( ூ ଵ
the number of individuals observed at their arrival points. The probability to observe an individual at ݕ is therefore:

)ݕ( = )ݕ(߱  ݃)ݔ(ߣ ఏ బ ݕ( -ݔ݀)ݔ ௫∈ௐ  ݔ݀)ݔ(ߣ ௫∈ௐ
with the set of arrival points in ܹ following a spatial Poisson process of intensity:

)ݕ(ߤ = )ݕ(߱ න ݃)ݔ(ߣ ఏ బ ݕ( -ݔ݀)ݔ ௫∈ௐ

Maximum Likelihood method (ML) leads to a biased estimation of dispersal kernel due to spatial heterogeneity in observation effort

A classical way of estimating the dispersal kernel is through the maximisation of the log-

likelihood ܮܮ = ∑ log /()ୀଵ ቈ ഇ (௬ ି௫ )  ഇ (௬ି௫ )ௗ௬
∈ೈ [START_REF] Clark | Invasion by extrems:population spread with variation in dispersal and reproduction[END_REF][START_REF] Gerber | High rates of gene flow by pollen and seed in oak populations across Europe[END_REF]. The law of large numbers implies that, when J → ∞, tends towards:

  )ݕ(߱)ݔ(ߣ ݃ ఏ ݕ( -)ݔ  ݃)ݖ(߱ ఏ ݖ( -ݖ݀)ݔ ௭∊ௐ log ( ݃ ఏ ݕ( -)ݔ  ݃ ఏ ݖ( -ݖ݀)ݔ ௭∊ௐ ݔ݀ݕ݀) ௬∊ௐ ௫∊ௐ   ݃)ݕ(߱)ݔ(ߣ ఏ ݕ( -ݔ݀ݕ݀)ݔ ௬∊ௐ ௫∊ௐ
However, because of the weightings by ,)ݕ(߱)ݔ(ߣ the ratio does not converge to a Kullback-Leibler divergence (up to an additive constant) and asymptotic unbiasedness is therefore not ensured [START_REF] Dacunha-Castelle | Probabilités et Statistiques: 2. Problèmes À Temps Mobile[END_REF]. Note that if ߱ is constant, that is to say when the observation process is spatially homogeneous (i.e. all individuals are detected with equal probability in any location of ܹ), then tends towards a weighted Kullback-Leibler divergence (up to an additive constant) of the form:

ܮܮ J →  )ݔ(ߣ  ݃ ఏ ݕ( -)ݔ  ݃ ఏ ݖ( -ݖ݀)ݔ ௭∊ௐ log ( ݃ ఏ ݕ( -)ݔ  ݃ ఏ ݖ( -ݖ݀)ݔ ௭∊ௐ ݔ݀ݕ݀) ௬∊ௐ ௫∊ௐ ܿ   ݃)ݔ(ߣ ఏ ݕ( -ݔ݀ݕ݀)ݔ ௬∊ௐ ௫∊ௐ
Thus, maximizing leads to asymptotically biased estimators of dispersal parameters unless observation process or departure point process is homogeneous.

Estimation of dispersal kernel using weighted Maximum Likelihood method (wML)

Below are the three successive steps we propose to account for the sources of bias identified above: the border effect and imperfect detection due to spatial heterogeneity in observation and marking efforts.

Step 1-Eroding the spatial window to minimise the bias induced by the border effect.

Let assume that the majority of dispersal distances are small relatively to the size of the spatial area ܹ considered, i.e. in a context of large-scale citizen science program. Let ܹ ∈ ܹ be the result of an erosion of ܹ by a distance ݀ so that the probability of an individual ݅ coming from ܹ and dispersing outside ܹ is small. We considered the subset of birds whose departure point lies in ܹ and arrival point in ܹ (we show in step 3 below how to correct for this potential bias).

Step 2-Estimating dispersal kernel parameters using weighted likelihood and accounting for local density of departure points.

The observation process ߱ cannot be considered as constant across space. To circumvent this issue, we weighted each dispersal event by the inverse of the density of departure points .)ݔ(ߣ By doing this, all points ݔ in ܹ have now a density probability proportional to ݃ ఏ ݕ( -)ݔ to be the departure point of the bird captured at

.ݕ Let ܼ ଵ (ߠ) = ∑ ୪୭ ( ഇ (௬ ି௫ )) ఒ(௫ ) /()ୀଵ
be the loglikelihood of observed dispersal distances weighted by the densities of departure points. Then

ܼ(ߠ) = ଵ ூ ܼ ଵ (ߠ) tends towards:  ௫∈ௐ బ  ݃)ݕ(߱ ఏ ݕ( -)ݔ log൫݃ ఏ ݕ( -)ݔ൯ ݔ݀ݕ݀ ௬∈ௐ  ݔ݀)ݔ(ߣ ௫∈ௐ ≈  ݕ݀)ݕ(߱ ௬∊ௐ  ݔ݀)ݔ(ߣ ௫∈ௐ න ݃ ఏ )ݑ( log൫݃ ఏ ()ݑ൯ ݑ݀ ௨∈ோ మ
which is a Kullback-Leibler divergence (up to a multiplicative constant

 ఠ(௬)ௗ௬ ∊ೈ  ఒ(௫)ௗ௫
ೣ∈ೈ and an additive one

݃ ఏ )ݑ( log ቀ݃ ఏ ()ݑቁ ݑ݀ ௨∈ோ మ
). The minimum of this divergence is obtained for ߠ = ߠ which ensures asymptotic unbiasedness, with ݑ a generic point anywhere across ܴ ଶ .

Here we want to estimate the parameter set ߠ by maximising ܼ(ߠ).

Random variables ܷ = )݅( ୪୭ ( ഇ (௬ ି௫ )) ఒ(௫ )
are independent and identically distributed. Let ߠ be the value maximizing .ܫ/)ߠ(ܼ Using a development similar to the ML method, one

obtains ߠ (ܹ) = ݊݅݉݃ݎܽ  )ݕ(߱   ߠ 0 ݕ( -)ݔlog( ߠ (y -ݕ݀ݔ݀))x ܹ∈ݔ 0 ܹ∊ݕ and √ܫ൫ߠ -ߠ (ܹ)൯
asymptotically converging towards a Gaussian distribution ܰ (0, ܦ ିଵ ′ܦܬܸ ିଵ ) where ܦ is the second derivative of ܷ(ܧ ) such as:

ܦ , = න න )ݕ(߱ ߲ ఏ ݕ( -)ݔ ߲ߠ ߲ ఏ ݕ( -)ݔ ߲݈  ఏ ݕ( -)ݔ ௬∈ௐ బ ݕ݀ ௫∈ௐ ݔ݀
with ܸ, the variance matrix of the ܷ gradient with respect to ߠ, such as:

ܸ , = න න )ݕ(߱ )ݔ(ߣ ߲ ఏ ݕ( -)ݔ ߲ߠ ߲ ఏ ݕ( -)ݔ ߲݈  ఏ ݕ( -)ݔ ௬∈ௐ బ ݕ݀ ௫∈ௐ

ݔ݀

Step 3-Correcting dispersal kernel parameters for the bias due to border effect Finally, ߠ (ܹ) → ߠ when ܹ → ܴ ଶ . In practice, the confidence interval around ߠ (ܹ) is computed by bootstrapping the dispersal data set, whereas the bias )ߠ(ܧ -ߠ (ܹ) is estimated as the difference between ߠ (ܹ) and the mean bootstrap parameter value [START_REF] Cook | Simulation-extrapolation estimation in parametric measurement error models[END_REF].

SIMULATING DISPERSAL DATA

We simulated a set of 5000 departure points by drawing coordinates from a non-stationary Poisson process including a gradient with densities declining eastward, within a 100×100 unit area. The detection probability was also spatially heterogeneous and set proportionally to the y-coordinates, therefore also producing a gradient with declining effort from north to south.

The combination of these two gradients (density of departure points, detection probability) created four distinct conditions spread over four quarters (starting from the north-west quarter, clockwise): 1. high nest density/high detection probability, 2. low nest-density/high detection probability, 3. low nest density/low detection probability, 4. high nest density/low detection probability (Fig. 1). These four conditions encompass the full range of spatial heterogeneity in data collection typically observed in large-scale monitoring schemes of animals. Dispersal distances were generated using a log-sech distribution (a = 0.75, b = 1; generating a median of 8.3 units). Dispersal events occurring outside the observation square were not considered. This simulation process was repeated 200 times to analyse the stability of parameter estimates.

A REAL CASE STUDY: DISPERSAL OF MONTAGU'S HARRIERS ACROSS FRANCE USING CITIZEN-SCIENCE

DATA The Montagu's harrier and its monitoring in France

The Montagu's harrier Circus pygargus is a medium-sized (ca. 300 g) ground-nesting raptor inhabiting open farming landscapes, within the Western Palearctic, and wintering in the Sahel region in Africa. About 3900-5100 breeding pairs are found in France [START_REF] Millon | Busard cendré. Rapaces Nicheurs de France. Effectifs, Distribution & Conservation[END_REF], essentially nesting in winter-sown cereal fields (wheat and barley). A substantial proportion of the broods (typically 2-5 chicks) does not fledge before crop harvesting takes place (20-50%; [START_REF] Millon | Comparative breeding biology of Hen Harrier and Montagu's Harrier: an 8-year study in north-eastern France[END_REF]. In response to this conservation issue, an extensive network of amateur ornithologists emerged from 1975, which undertakes the protection of the nests in close collaboration with farmers [START_REF] Santangeli | Identifying effective actions to guide volunteer-based and nationwide conservation efforts for a ground-nesting farmland bird[END_REF]. Between 2008 and 2015, an average of 1083 ± 191 nests were monitored annually (range: 681-1287).

Data from several small-scale ringing schemes suggested that if most birds tended to return in the vicinity of their natal nests, both natal and breeding dispersal can incur distances of several hundreds of kilometres. In 2007 was launched a wing-tagging program across the whole country, with the objective of unravelling the demography of the species at large spatial scale, including the estimation of dispersal kernels. Building upon the existing network of volunteers, a total of 5129 chicks have been wing-tagged from 1778 nests, between 2007 and 2010, using coloured plastic tags attached to the patagium of each wing (see www.busards.com for a complete description of the scheme & Appendix B).

Dispersal data

A dedicated website, including a Google© map interface with automatic recording of the geographical coordinates, has been launched so as to ensure an efficient collection of resightings across France and more widely across Europe and Africa (www.busards.com). We filtered the >6000 resightings recorded in the database between 2008 and 2013 so as to retain a single resighting per year (excluding the year of birth) for each bird observed during the breeding season (15 April-15 July, i.e. excluding movements during the migration period). In the case of multiple resightings of a same bird during the same breeding season, we selected the first occurrence, so that each bird can only produce one dispersal distance per year, but up to six dispersal distances over the study period. Birds from a same nest were considered as independent in regard to dispersal event. The actual reproductive status and sex of the birds when resighted were not considered. Our observation area was defined as France, thus observations made outside France were excluded from the dataset (5%, N = 68), and we ended up with a dataset gathering 1215 dispersal events.

We applied the aforementioned method to the harrier dispersal data. ܹ was defined as the minimum convex polygon containing all harrier nests. ܹ ∈ ܹ was an eroded area including nests located at ≥ 150 km from the border of ܹ. The threshold of 150 km was chosen from the dispersal distribution, so as to include ca. 80% of the observed dispersal events. We chose a heavy-tailed distribution, namely the log-hyperbolic secant, or log-sech distribution, to describe harrier dispersal distribution [START_REF] Van Houtan | Dispersal of Amazonian birds in continuous and fragmented forest[END_REF]. We assumed that dispersal distances ݀ ≥ 0 follow a log-sech distribution such that:

݃(݀) = 2 )ݎܾߨ( ⁄ )ܽ/ݎ( ଵ/ + )ܽ/ݎ( ିଵ/
where ܽ > 0 is the scale parameter (expressed in distance units), and ܾ > 0 sets the shape of the kernel (no dimension). We excluded the locations at sea from the area ܹ. Local nest density was calculated by a two-dimensional kernel density estimator [START_REF] Venables | Modern Applied Statistics with S, Fourth edi[END_REF]. Simulations and analyses were performed in R 3.1.3 (R Development Core Team 2015). Descriptive statistics (mean or median) are given together with their 95% confidence intervals. The code for conducting such analyses is fully available upon request.

RESULTS

ASSESSING THE PERFORMANCE OF THE METHOD USING SIMULATED DATA

As expected, spatial heterogeneity in marking and observation efforts (Fig. 1) led to a biased dispersal kernel, with relatively long dispersal events being under-represented (Fig. 2). The two parameters of the estimated log-sech distribution were biased (scale: a = 0.85, shape: b = 0.88) compared to the parameters used to simulate the data (a = 0.75, b = 1). Application of our method accounting for spatially heterogeneous marking and observation efforts did largely correct for this bias (mean value with 95% CI: a = 0.69 [0.49-0.93], b = 0.99 [0.94-1.05]). If the scale parameter a appeared as slightly underestimated, this had only little consequences in terms of departure from the simulated data (Fig. 2).

DISPERSAL OF MONTAGU'S HARRIERS ACROSS FRANCE

Considering the whole dataset of dispersal events occurring within France (N = 1215), Montagu's harrier dispersed over relatively long distance with median dispersal distance of 32 km and maximal dispersal distance reaching 627 km (Fig. 3). When spatially restricting the spatial scale of the dataset so as to limit the border effect (ܹ , N = 557), the distribution of dispersal events remained largely unchanged with a median distance slightly increasing up to 37 km (max. distance = 602 km; see Appendix C comparing histograms and quantile distributions between the two datasets). Parameters of the log-sech distribution fitted to the observed dispersal data (ܹ ) were a = 11.5×10 2 km and b = 0.66.

Implementing our correction method however led to a strong increase of the shape parameter b = 0.95 [0.90-1.03], whereas the scale parameter a = 11.5×10 2 km [10.0×10 2 -12.7×10 2 ] did not change. Corrected median dispersal distance was considerably increased and reached 123 km [105-160] (Fig. 4). The interquartile range lay between 34 and 473 km.

Note that the cumulative frequency only reached 0.8 at 700 km, i.e. at roughly the maximal dispersal distance observed within France, suggesting that another 20% of dispersal events are likely to be recorded further away (maximum fixed at 2000 km). In order to more directly compare the corrected kernel with the dispersal events observed (Fig. 3), we restricted the analysis up to a distance of 700 km. At this reduced spatial scale, the corrected median dispersal distance was 78 km [70-101] (interquartile range: 27-227 km, Fig. 4).

DISCUSSION

Here we proposed a method for deriving unbiased dispersal kernels over continuous space from data collected at large spatial scale. Such data typically suffer from high spatial heterogeneity in both marking and observation efforts, in addition to a border effect. Our method successfully corrected dispersal kernel parameters by accounting for heterogeneity in the distribution of departure points (i.e. local nest density) within an eroded spatial window.

Using simulations, we checked that this method was indeed able to provide unbiased and precise estimates of a probability density function fitted to dispersal distances. We further applied it to data collected on Montagu's harriers from a nation-wide citizen-science program and showed that corrected median dispersal distances were 2.4-3.8 larger than the median of observed dispersal distances.

ACCOUNTING FOR SPATIAL HETEROGENEITY IN MARKING AND OBSERVATION EFFORTS, AND BORDER EFFECT WHEN FITTING DISPERSAL KERNELS

From a technical viewpoint, the method we used is akin to a weighted likelihood method but present major improvements for dealing with the two classical issues encountered when modelling dispersal data: i) the finite nature of a study area (border effect) and ii) the spatial heterogeneity in detection probability of marked animals linked to marking and observation efforts. Indeed, all dispersal events within the study area cannot be observed. Therefore, in the contrast estimation function we used, the likelihood of each dispersal event was weighted by the inverse of the local nest density around the departure points, a critical advance since departure points (i.e. marking locations) are known while observation effort, the parameter necessary in alternative methods, is not.

The logic here was to account for the local nest density (i.e. potential departure points) when modelling the probability of birds to arrive at a given point. Indeed, by weighting the number of observed dispersal events by the number of departure points, and not the number of recruitments (see e.g. [START_REF] Baker | Correcting biased estimates of dispersal and survival due to limited study area : theory and an application using wrentits[END_REF], we avoid the need of collecting data on unmarked individuals and the direct modelling of observation effort, typically intractable using citizen-science data. Data on recruitment are particularly costly to collect, even for species nesting in a defined set of nestboxes (e.g. [START_REF] Van Noordwijk | On bias due to observer distribution in the analysis of data on natal dispersal in birds[END_REF][START_REF] Winkler | The natal dispersal of tree swallows in a continuous mainland environment[END_REF].

Our method does not require information on recruitment other than for marked individuals.

Alternative methods such as the distance-dependent recruitment rate allow to test for departure of the dispersal rule from a uniform law (i.e. birds chose their breeding location at random), and is robust to spatial heterogeneity in observation effort (van Noordwijk 2011).

However, this method cannot be used to derive dispersal kernels (Appendix A).

Second, the bias due to the finite nature of the observed area was corrected by a bootstrap procedure so as to get the unobserved tail of the distribution, assuming the form of the distribution was scale-independent. The advantage of this approach is that inference can be made on the frequency of dispersal occurring outside the study area. This however comes at the expense of relying on a particular statistical distribution. We used the log-sech distribution for its ability to fit heavy-tailed dispersal kernels [START_REF] Van Houtan | Dispersal of Amazonian birds in continuous and fragmented forest[END_REF]. Note however that our method can be fitted with alternative distributions (e.g. negative exponential, gamma, lognormal), and that a procedure can be performed to select the most appropriate distribution.

LIMITATIONS AND POSSIBLE EXTENSIONS OF THE METHOD

Our method relies on a non-stationary spatial Poisson distribution of departure points but can be extended to incorporate any point process distribution. However, the bootstrap step will be more tedious as it will need the simulation of the chosen distribution of departure points. We also assumed that the rules governing dispersal processes are stationary, i.e. are constant across space or across ecological conditions encountered by individuals, and that all individuals have a similar average dispersal behaviour. If dispersal behaviour varies among individuals, for example according to age or density of conspecifics in the natal area, the model can then be extended to investigate the effect of a covariate on the kernel parameters.

Let )ݔ(ݒ be this driving variable (continuous or categorical), the new kernel then becomes

݃ థ ᇱ ,ݔ( )ݕ = ݃ ఏ൫థ,௩(௫)൯ ,ݔ( )ݕ
where ߶ is the parameter set linking ߠ to .)ݔ(ݒ A comparison of dispersal kernels between birds born in low-density areas vs. birds born in high-density areas could then be made. Our method is also flexible enough to explore whether habitat selection, the final stage of the dispersal process [START_REF] Clobert | Dispersal -Causes, Consequences and Mechanisms of Dispersal at the Individual[END_REF], depends on local environmental suitability. In this case, let ℎ(߶, ))ݕ(ݑ be the probability to settle in site ݕ with environmental suitability .)ݕ(ݑ The model can then be extended by using:

݃ థ ᇱ ,ݔ( )ݕ = ݃ ఏ(௫,௬)(థ,௬)  ݃ ఏ(௫,௭)(థ,௭)ௗ௭ ௐ A LARGE-SCALE CASE STUDY USING CITIZEN-SCIENCE DATA: MONTAGU'S HARRIER DISPERSAL ACROSS FRANCE
Accounting for spatial heterogeneity in both marking and observation efforts increased the median dispersal distance from 32-37 km (using observed data) to 78-123 km, depending on the chosen spatial scale. Whether considering France or Western Europe, the interquartile range of dispersal approximately doubled, illustrating the necessity to clearly specify the population of interest when summarising dispersal propensity. Actually, the cumulated frequency of dispersal distances derived from the corrected log-sech distribution did not reach one at the maximal observed dispersal distance and the corrected log-sech distribution predicted that roughly 10-20% of dispersal events may occur beyond this distance (Fig. 4).

This may appear as substantial but it is important to note here that a number of dispersal events have indeed been recorded in neighbouring countries during the breeding season (Netherlands, Belgium, Germany, Switzerland, Scotland, Ireland, Czech Republic, Hungary, Spain) and were excluded from our analysis. Fourty-eight dispersal events from 34 different individuals have been actually recorded in neighbouring countries (4% of the data; Appendix D). Twelve individuals dispersed at >700 km from their birth nest (median = 561) with a maximal distance recorded at 1366 km (a bird resighted in Scotland).

ON THE USE OF PASSIVE MARKS TO STUDY DISPERSAL

Estimating animal dispersal has been identified as a major methodological challenge in ecology for decades [START_REF] Doligez | Estimating fitness consequences of dispersal: a road to 'knowwhere'? Non-random dispersal and the underestimation of dispersers' fitness[END_REF][START_REF] Schaub | Estimating true instead of apparent survival using spatial Cormack-Jolly-Seber models[END_REF]. Capture-recapture estimates based on marked animals have been criticised because of spatial sampling bias.

Alternative approaches have thus been proposed, but all have their limits. For instance, despite dispersal can be inferred from genetic data and provide insightful information about the spatial scale of gene flow, such methods cannot help to infer the shape of dispersal kernels and to assess the demographic contribution of dispersal, partly because they integrate dispersal events over several generations. Satellite tracking (Argos and GPS systems) has been presented as the panacea for the study of dispersal [START_REF] Koenig | Detectability, philopatry, and the distribution of dispersal distances in vertebrates[END_REF]. However, economic cost limits the ability to investigate juvenile dispersal (since young individuals suffer higher mortality rates compared to adults), and device weights further restrict the number of species that can afford to carry such material [START_REF] Wikelski | Going wild : what a global small-animal tracking system could do for experimental biologists[END_REF]). Therefore we argue here that passive marks remain a major tool for studying dispersal, especially combining citizen-science program allowing studies to be conducted at very large spatial together with the statistical method we propose here which correct for sample bias including heterogeneity in marking and observation efforts as well as for border effect.

Figure 1. Map of the simulated dispersal events (dark blue dots: departure points, light blue dots: arrival points). The simulations were generated so as to create spatial heterogeneity in the combination of nest density and detection probability, as typically observed in field data (see Methods). Note that many short-distance dispersal events are actually hidden by the dots.

  550

Figure 2 .

 2 Figure 2. Cumulated frequency of dispersal distances (in arbitrary unit) i) predicted from the

Figure 3 .

 3 Figure 3. Observed dispersal events of Montagu's harrier between natal nests (all nest

Figure 4 .

 4 Figure 4. Observed (black line) and modelled (red line, grey envelope for 95% CI) cumulated

Methods in Ecology and Evolution

ACKNOWLEDGMENTS

We are immensely thankful to the network of volunteers working on Montagu's harriers in France and sorry not to name the hundreds of them that contribute to the conservation of the species. This work has greatly benefited from the critical review of E. Gabriel and R.

Senoussi. We also thank A. van Noordwijk for his constructive comments.

DATA ACCESSIBILITY

Harrier data and R code are available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.g4n5q.

AUTHORS CONTRIBUTIONS

STATEMENT JC, AM, VL & VB designed the research. JC did the mathematical development and ran analyses and simulations. AM wrote the paper. JC & VB contributed to revisions of the manuscript. JLB, TP, BVH, AM & VB ensured the coordination of the volunteers' network and provided the data.


Faten Ajmi

Sarah Ben Othman

Hayfa Zgaya Biau
email: hayfa.zgayabiau@univ-lille2.fr

Slim Hammadi
email: slim.hammadi@centralelille.fr

Patient Pathway Workflow Model Identifying Overcrowding Indicators in Emergency Department

Keywords: Emergency department, Workflow, key performance indicators, overcrowding situation, Patient pathway

The overcrowding situation in emergency department is serious problem worldwide. In fact, the majority of developing countries are interested in developing their hospital system in order to anticipate the overcrowding situations and ameliorate the quality of patient care. In this paper, we focus on patient pathway in an Adult Emergency Departments (AED). This process presents the main burden for overall system operation due to its growing dynamics and the distributed organizational structure of hospitals. In order to determine the key performance indicators for the internal functioning of AED, we are building a workflow model based on the real patient pathway to direct the patients into the dedicated service while respecting the quality of care according to the expectations and requirements of the patient. Our approach is tested on a set of real database from AED of CHRU in Lille (France). Simulation results show that the proposed model can determine 3 keys performances indicators to optimize in order to improve the patient pathway in the AED.

INTRODUCTION

The ultimate goal of health care systems is to provide accurate, effective and rapid treatment of patients [START_REF] Kırış | A knowledge-based scheduling system for Emergency Departments[END_REF], particularly in emergency departments (ED). ED was first recognized in the 20th century [START_REF] Md | Anyone, Anything, Anytime: A History of Emergency Medicine[END_REF]. In France, and in the world, ED provides primary health care to patients whose diseases or symptoms require immediate attention and in some cases can put the patient life in danger.

In this paper we will focus only on the Adult Emergency Department (AED) of CHRU of Lille (city of north of France). From an outside view, an emergency service should be organized as an anthill. That is, the reception and treatment of emergencies are thought and coordinated. Consequently, the functioning in the emergency services must be perfectly organized and the team spirit particularly present. However, from a closer view the internal functioning of healthcare establishments in general and of ED in particular is much more complicated to master. Among the major problems of this system we can find, on the one hand, the lack of coordination between the various services or even in the same department and on the other hand the problem of the complexity of the patient pathway management, since an ineffective patient treatment procedure can unfortunately lead to pernicious results.

Part of the complexity of this care system is that the arrival time of the patients and their pathologies are impossible to predict. In addition, we have an incoming flow that is totally random. So in order to improve the efficiency of the service by identifying several performance indicators, we have to develop a reliable model of patient pathway in order to identify the key factors of dysfunction of the AED. The fig 1 summarizes the patient pathway from the moment of his arrival in AED until his first medical consultation "taken care by doctor". So after a significant waiting period depending on the crowding status of the AED and the urgency of the patient case, the patient is finally taken in charge by a doctor.

STATE OF THE ART

In order to improve the process of patient care in AED, as well as to optimize medical staff care activities, it's necessary to model the real patient pathway within the department studied. The idea is to orient patients to the dedicated service, allowing us to trace the flow of information and to kwon the material occupied in the service. Obviously maintaining a regular flow of information is the key to take a good decision for system evolution, hence the interest of setting up a connected Workflow system.

In the literature, the most adequate definition of Workflow for our work is given by (Thomas, n.d. 1995):"Workflow management is done by a proactive system for the management of a series of tasks which are transmitted to the appropriate participants in the right order and which are completed under allocated time". The application of this definition in the hospital system has been successful, which rendered the health system a very active area of research (RJ Salway et al. 2017), (Andrew S et al. 2018). The majority of this research has addressed about the subject of providing an efficient health care system for the patient and ensuring its continuity for the next generation [START_REF] Chen | A Survey of Recent Developments in Home M2M Networks[END_REF] [START_REF] Wan | Cloud-enabled wireless body area networks for pervasive healthcare[END_REF] [START_REF] Hossain | Cloud-Supported Cyber #x2013;Physical Localization Framework for Patients Monitoring[END_REF]. Thus, this system has undergone several changes over time, due to the advancement of technology up to the application of the smart home [START_REF] Chen | A Survey of Recent Developments in Home M2M Networks[END_REF]) [START_REF] Deen | Information and Communications Technologies for Elderly Ubiquitous Healthcare in a Smart Home[END_REF] [START_REF] Saleem | Human-oriented design of secure Machine-to-Machine communication system for e-Healthcare society[END_REF]. [START_REF] Bricon-Souf | Dynamic workflow model for complex activity in intensive care unit[END_REF]) used the Workflow system to model productive work in the intensive care unit, explaining that the flexibility of this system is crucial. [START_REF] Van Der Aalst | Workflow mining: A survey of issues and approaches[END_REF] proposed a Workflow mining, in order to facilitate the conception of the Workflows. [START_REF] Esdar | Understanding latent structures of clinical information logistics: A bottom-up approach for model building and validating the workflow composite score[END_REF] used the workflow to model the logistics of clinical information by aiming to describe and explain various phenomena of information provision in order to conduct clinical processes. [START_REF] Lee | Design and development of logistics workflow systems for demand management with RFID[END_REF] Examined how Artificial Intelligence technologies and RFID technology can improve workflow in logistics supply chains. [START_REF] Narock | An agent-based approach for capturing and linking provenance in geoscience workflows[END_REF] used the multi agent system with multiple agents to identify an efficient and evolutionary system in the provenance field. Besides, new approaches try to model and optimize simultaneously these patient care pathways using scheduling and workflow techniques joined into the same system (C. L. [START_REF] Meli | Load-sensitive dynamic workflow re-orchestration and optimisation for faster patient healthcare[END_REF].

Major problems in the management of health care production systems concentrate primarily up to: the control of hospital flows, anticipations of overcrowding and the control of process flows (Physical, information, equipment, financial...) [START_REF] Von Thiele Schwarz | Efficiency in the emergency department -A complex relationship between throughput rates and staff perceptions[END_REF]. Including restructuring which consists of the pooling of resources and the creation of a technical platform. However health professionals do not have adequate training to solve such problems, given their lack of methodologies and tools for decision-making and management adapted to the requirements implied by their external environment [START_REF] Ozkaynak | Characterizing workflow for pediatric asthma patients in emergency departments using electronic health records[END_REF].

Thanks to the collaboration between clinicians, healthcare, academic and industrial researchers, it is possible to model the current situation in AED in order to improve the quality patient care. Indeed, the model created will help us to identify situations of overcrowding in the service (thanks to the identification of corresponding indicators of overcrowding). Then, according to this modelling, we can establish a decision support system to provide strategies for avoiding overcrowding situations. Hence, it is essential to evaluate the relevance of the predictive models proposed below, then to develop recommendations for strategies for avoiding overcrowding situations. In this article, we focused on research based on the use of Workflow as a suitable model to represent the patient pathway and identify the key factors from the overcrowding situations.

The implementation of a structured research model is made thanks to the following points:

 Several visits were realized in the AED of CHRU of Lille with researchers and physicians allowing us to build a workflow model of the patient pathway and AED treatment processes as the support for a patient mapping form to collect database.  Each research assistant spent a continuous 15 days mapping exercise during which they independently observed and manually recorded patient pathway at the AED for a sample of presentations.  In order to validate an all-day observation, an AED staff meets the research assistant for controlling their observations of patient flow.  The analysis done by the researchers and AED staff medical conducted us to construct a process models. This model will be constructed with the BPMN graphical language which is an international standard.  The last point deal with the validation and simulation of the model. A comparison real patient pathway scenarios is necessary to identify the overcrowding indicators.

The Methodologies Included

The use of the Workflow approach and in particular the BPMN "Business Process Model and Notation" graphic tool, allows to represent the specificity of the business processes and to verify and validate this process directed by the staff medical.

The workflow presents a graphical model that translates the flow of the patient into the service. In addition, the use of the internationally standardized language (BPMN) guarantees common understanding by different users. There are 4 categories of BPMN elements, which represent the basic elements of graphical modelling, see table 1.

We have already compared the workflow model with the UML and Petri net model. We note that the wokflow model with the "Bonita Soft" tool is better thanks to an integrated engine, the automatic generation of the interfaces, the simulation of the modeling put in place and the possibility of connecting thanks to the connectors the Workflow system to other external systems.

The idea of modelling the patient pathway by using a computerized system "Bonita Soft" tool aimed as a first step to determine the bottleneck that presents the cause of global dysfunction in the service, as a second steps to:

 Ensure the control of all transactions planned at the moment, thanks to the insertion form and the possibility of interaction between the medical staff regardless of their location.  Ensure a fast registration of patients and their data in the system.

 Quickly identify the availability of medical staff wanted as needed.  Trace his pathway, therefore possibility to calculate the time planned for the remaining of his treatment.

Table 1. The BPMN graphical language

Description of CHRU of Lille

The CHRU of Lille is the third biggest hospital in France It includes a pediatric and adult emergency department. Hence, the complexity of control its management in different sectors, the most studied are the financial sector and the sector of the take care of patients. As we mentioned at the beginning, in this paper we study the modelling of the patient pathway in this huge hospital, especially in the AED. AED dispose materiel radiology, MRI and ultrasound providing the patient with maximum in situ care, there are 5 operating rooms, one recovery room and one surgical resuscitation room.

Patient care process in the AED

From the patient arrival in the AED he must undergo this process: 1-The patient has registered in the services electronic database. (Short Circuit or Long Circuit). 3-The patient is waiting to be treated.

There are different patients care processes depending on the severity of patient pathology.

Thanks to the visits of AED of CHRU Jeanne of Flanders at Lille. We can present the modelling of patient pathway in the AED as faithfully as possible in order to identify the various dysfunctions in this process. Especially to better encompass the problems previously faced by the service.

WORKFLOW MODELS

The model presented in this paper is generic to represent any patient pathway at the AED. Each healthcare unit is independent, in other words visiting one unit does not influence the patient visiting another unit. However, each instance of the model represents a full part of pathway care for each treatment. Thus, the same patient passes from one stage to another in the care process.

Although creating a Workflow process is a complex step and takes a lot of time it is obvious to find gaps between the real workflow processes and the scientific Workflow process performed. In this paper, we present the process of patient care in AED most accurately as possible. d-Out-of-CHU hospitalization: send patients to another hospital or clinic due to lack of places. e-Inadequate hospitalization: in the case of total bottleneck and for reasons of complexity of patient assignment in a service the patient can be sent to an inadequate service. f-The patient returns home: directly after a full consultation or after a hospital stays lasting less than 24 hours 3-Two tasks for patient reception: one dedicated for the entry of SAMU the other for the standard reception (for valid patient : who no need stretcher) 4-The surgical treatment shock task is to treat patients who have a very serious health condition that should not wait.

If the patient is treated less severely he will go through the waiting task (x minutes) and then to the reception and orientation zone (ROZ) task to give him a ground of appeal.

Depending on the severity of his state of health, he will be sent either in Short Circuit SC or Long Circuit LC. The same step for valid patients or patients on a stretcher (SAMU patient).

At this stage we have two different types of pathway depending if the patient is on a stretcher or he is valid. However, it all depends on the state of the service at that moment, the patient can wait or not wait before being taken care. Then the patient will be oriented or reoriented to one of these units where some of them are represented by a subprocess in our model:

5-Long Circuit (LC) task; Sub-process. 6-Short Circuit (SC)task; Sub-process. 7-The Waiting task after the "ROZ": for only the patient who is placed on a stretcher, they wait that there is available box to call them in the LC. 8-Operating room task. 9-Very short-duration hospitalization zone VSHZ task. 10-Short-duration hospitalization unit SDHU task.

We worked in an environment of uncertain care that increases the risks during a patient pathway. This uncertainty confirms the complexity of the emergency care system. These decisions are expressed in the model by six synchronization decision points Dpi (1< i <6) (see Fig. 2) these graphical connections are used to control the behavior of sequence flows in a process. Decision points Dpi (1 < i <6) of our global model are detailed as follows:

Dp1: represent the type of patient arrival. Valid patient or patient coming by SAMU "urgent" if the patient is urgent it goes directly to the SAMU reception where he is going to be directed either towards the medical shock or the surgical shock. If he is valid he goes directly to the standard reception.

Dp2: control the orientation of the patient according to his / her state of health either towards medical shock or the surgical shock or the waiting depending of its health state.

Dp3: This point is important because it determines the number of patients in the service who will subsequently oriented in the LC or in the SC. It is essential because it represents the input bottleneck point of the global model.

Dp5, Dp4: These points are central because they represent the second input bottleneck point of the global model. They guide patients in appropriate care according to their pathologies. In case of blocking of these points we will have "a zone E, LC and SC" complete / overpopulated and we cannot accept new entries in the LC/SC.

Dp6: this point is essential because it represents the output bottleneck point of the global model. It redirects the patients to one of these points:

 Post-emergency service: it is managing by interns 'physicians. It only takes emergency patients, includes 20 beds for a total of 50 patients in the AED per day.  Adapted service: concerned patients who are hospitalized in the VSHZ or SDHU and then sent in order to continue their treatment in a service adapted according to their pathologies  Home: the patient returns home with good quality of care. Also allows the exit of a patient and the recovery of another patient.

Sub-process SC

If the patient has a small medical problem, he is oriented to the SC which supports valid patients under the age of 65 who do not need for stretchers or extensive biological analysis. This is represented by a BPMN sub-process (call task).

Sub-process LC

If the patient is oriented to the LC, his is then oriented to a special care zone X{A,B,C,D} From the special care zone X{A,B,C,D}, the patient is oriented to the Urgent Box Unit (UBU) for urgent cases and to other box units, otherwise.

A set of decision points DSi (1≤i≤12):  DS1: Allows to orient the patient to the zone where it was sent from the global model  {DS2, DS3, DS4, DS5}: Separating the pathway of valid patients from the patient pathway most urgent one.

These DSi {DS1, DS2, DS3, DS4, DS5} allow to control the common process of the patient pathway in the LC from this each patient will follow its own pathway.

 {DS6, DS11}: allows the control and the differentiation of the paths of the patients for the remaining of the valid patient treatments  {DS7, DS8, DS9, DS10, DS12}: allows the control and the differentiation of the paths of the patients for the remaining of the urgent patient treatment.

The BPMN has some limitations to model a complex processes, especially the ability to represent queues and data-driven decision points. A tow interesting work in the literature has been developed to address these limitations:

 a conceptual design for an extension to BPMN (BPMN4SIM) using model-driven architecture.

SIMULATION & RESUTS

Based on the real data in the AED at CHRU of Lille, We note that it is rather difficult to collect exact and precise data notably because there are a many data censored. However, we selected two types of situations: Normal and Overcrowding Situations (NS, OS). In order to model the patient pathway in the AED, we choose firstly the workflow as an adequate model to represent these patient care activities in order to better understand and identify the dysfunction zones in the AED. Secondly, thanks to this model we use real data for analyse and identify several key performances indicators (KPIs) such us: Length of Stay indicator (LOS), Care Load Rate (CLR), and Left without Seen Rate (LSR). These KPIs are needed to indicate any possible issue with the current system. The evaluation of these indicators is crucial to determine in the AED the overcrowding situation. In addition, the analysis of the database and the modelling of the patient pathway presented in this paper are complementary to ensure the robustness of the solutions obtained in practice.

LOS: Represents the period that the patient spent in the service from its entry to its exit from the emergency system. CLR is the occupancy rate by the medical staff during a care treatment. The observation of the CLR curve suggests that the nature of its evolution is associated with the evolution of patient flow according to time. We note that a patient arrival flow of 5% following by another 10% in less than 20 minutes provoked a 10% increase for the CLR. This increase will last more than 35 minutes, imply that among these some of them have severity pathologies The LSR is the percentage of patients who leave the AED after the ROZ task and before the first medical consultation. Fig. 7 shows a comparison of a cumulative distribution of the LSR between NS & OS situations. We note that the LSR in OS begins to increase from the first 40th minutes of waiting; imply that for psychological raisons, the patient observes that as more and more new patients arrive in the AED, which will become more and more overcrowded. Although in the NS, the LSR begins to increase from the 70th minutes of waiting, imply that in NS the waiting time may increase excessively because of number of urgent patients to be treated as a priority which will increase the CLR indicator and makes the medical staff unavailable.

DISCUSSION

Improving the performance of any health care system requires first to know how to measure its performance properly. The different made visits to AED, the observations and the discussions with clinicians and the pathway modelization prove that the chosen KPIs (LOS, CLR, LSR) are suitable to study the AED overcrowding situation. Thanks to our model and to simulations results we conclude that the patient waiting time in AED may be excessive in the two situations OS and NS, because of the massive arrival of patient in AED or/and patient arrival in a vital emergency makes medical personnel unavailable. Then we note that the chosen KPIs are strongly correlated to congestion and patient satisfaction. The selection of the appropriate KPIs has always been a controversial subject, for the whys and wherefores remain unclear. Neither the scientific community nor practitioners are able to decide about the most appropriate KPIs, as each indicators presents at the same time benefits and drawbacks. For example, the LOS gives an overview on the entire system performance but does not allow to figure out local strengths and weaknesses. Another example the measure LSR depend on external factors on which the AED has no control. We note also, that the CLR is a crucial KPIs for the critical levels of the physician workload, but it does not give any information about the remaining care operations to be done on the patients during other important steps of the process (beyond the first consultation). We will therefore review studies that combine different KPIs. We will show relevant combination of KPIs and highlight potential interdependency between them in order to focus some universal quantitative measures of AED overcrowding. To minimize patient waiting times and optimally organize the medical staff care activity we need in addition to develop scheduling and orchestration approaches in order to minimize the waiting time of the patients during the workflow execution.

CONCLUSION

The use of the Workflow approach to model the patient pathway in the AED at CHRU of Lille was effective in facilitating understanding of the set of flow processes in the AED. Thanks to the global model and the two sub-processes we model the complete patient pathway, so with the simulation we identify the 3 key performance indicators (LOS, CLR, LSR) that their real time control is essential especially during an overcrowding situation at the service. This paper shows the complexity of modelling the patient pathway in an emergency department as large as the CHRU in Lille as well as the importance of determining the key performance indicators. So we conclude that overcrowding situation can be represented as an imbalance between the flow of care load and the ability to take care of patients over a period of time, it causes not only an excessive waiting time for the patients but also results harmful consequences for the proper functioning of the AED.

Fig 1 .

 1 Fig 1. Patient care pathway

Fig 2 .

 2 Fig 2. Global workflow for patient pathway model in the AED Fig.2 presents the simplest overall process of the patient pathway with 2 large sub-processes characterizing two famous circuits: Long Circuit see Fig.4 and Short Circuit see Fig.3. Our workflow is composed of several BPMN components including start and end events, several kinds of tasks and

Fig 3 .

 3 Fig 3. Short Circuit Workflow for patient pathway in the AED

  ( B. S. S. Onggo et al. 2017)  a conceptual design for a multi-agent system architecture in which each agent has a real time information about the patient such as the evolution of patient care state in order to orchestrate the workflow (S. Ben Othman et al. 2016)

  Fig 6. The CLR evolution according to the patient flow and time

Figure 5 .

 5 Figure 5. Cumulative distribution of the LOS between NS & OS

  Half an hour later, the nurse after sorting the different patients in order of gravity, orients each of them to the most appropriate circuit.

	2-
	Workflow	Activities
		Events,
		starts events,
		end events
		Flows
		Gateways
	Organisation Swim lanes
		or lanes
	Readability	Text
		Annotation
	Specifics	Int. Events
	Behaviours	Timers
		Mistakes
		loops

ACKNOWLEDGEMENTS

This work was supported and funded by the Federative Research Structures Technologies for Health and Drugs (SFR-TSM): http://sfr-tsm.eclille.fr/