Ali Boudani
email: aboudani@irisa.fr

Bernard Cousin
email: bcousin@irisa.fr

Jean-Marie Bonnin
email: jm.bonnin@enst-bretagne.fr

MPLS Multicast Traffic Engineering

In this paper, we present multicast traffic engineering and we compare it to unicast traffic engineering. We study the advantages given by the integration of multicasting and MPLS. We present current proposals for multicast traffic engineering and using MPLS network. We describe our approach, the MPLS multicast tree (MMT) protocol. In order to reduce forwarding states and enhance scalability, MMT utilizes MPLS LSPs between branching routers of the multicast tree. We present a simulator for MMT and finally we discuss some simulation results. 1

I. INTRODUCTION

Best-effort service model existing in the current Internet is inadequate in meeting the growing demands of the next generation applications. These applications require QoS guarantees and effective data delivery. In consequence, the network is required to provide various qualities of service (QoS) for applications sensitive to delay, jitter and packet loss.

In order to provide QoS to users across the Internet, either we could increase the bandwidth available in the network such that the extra capacity allows all users to meet their appropriate QoS, or we should suppose that bandwidth is limited and therefore network resource should be appropriately allocated among users. Thus, some form of resource provisioning is necessary to provide QoS across the Internet. One of the most promising tools for providing QoS across the Internet is traffic engineering (TE).

TE improves the management of data traffic within a network and in consequence provides better utilization of network resources. In the best effort service model, available network resources are not being used efficiently, resulting in higher delay and lower bandwidth while TE could provide better quality of service by reducing delay and packet losses and increasing throughput experienced by end users using the same network infrastructure. This results in a minimization of the vulnerability of the network to service outages arising from congestions or failures occurring within the infrastructure. In addition to the quality of service, the factors driving the need for better TE tools include interdependent tunable parameters, network growth, traffic variability and multicasting.

An important parameter in network performance is the link loads. Many applications, like video/audio on-demand or teleconferencing, can consume a large amount of network bandwidth because of, first the volume of the transmitted data and, 1 This work is partially funded by the LAFMI.

second the larger number of application members. Multicasting is a useful service to support such applications. When a multicast service is used to send a packet to several destinations, a single transmission is required on any link of the multicast tree associated to the application; while several independent transmissions would be required using unicast service.

MPLS [START_REF] Rosen | Multiprotocol label switching architecture[END_REF] as a traffic engineering tool has emerged as an elegant solution to meet the requirements of the backbone networks of the Internet. Multicast and MPLS are two complementary technologies: multicast trees could be supported by MPLS networks. MPLS will enhance the network performance and present an efficient solution for multicast scalability and control overhead problems. Multicast attempts to conserve network bandwidth, while traffic engineering attempts to provision the bandwidth in an appropriate fashion to users. II. TRAFFIC ENGINEERING TE is the process of controlling how traffic flows through a network in order to facilitate efficient and reliable network operations while simultaneously optimizing network resource utilization and traffic performance [START_REF] Xiao | Traffic engineering with MPLS in the Internet[END_REF].

TE is needed in the Internet mainly because current interior gateway protocols (IGPs) always use the shortest paths to forward traffic. While the shortest paths approach is very simple to scale to very large networks and conserves network resources, it does not always make good use of these resources and may also cause the following problems:

• The shortest paths from different sources overlap at some links, causing congestion on those links.

• The traffic from a source router to a destination router exceeds the capacity of the shortest path, while a longer path between these two routers is underutilized. The first problem can be solved by expansion of link capacity, or by application of classical congestion control techniques, or both. Classical congestion control techniques attempt to regulate the demand so that the traffic fits onto available resources. Classical techniques for congestion control include: rate limiting, window flow control, router queue management, schedule-based control, and others [START_REF] Yang | A taxonomy for congestion control algorithms in packet switching networks[END_REF]. The second problem, namely congestion resulting from inefficient resource allocation, can usually be addressed through TE.

A constraint-based routing (CBR) and an enhancement of existing IGPs may be needed to permit unicast forwarding through explicit routes.

To build the multicast tree, multicast TE uses some estimation of network resource utilization, constraint-based routing algorithm and explicit routes, when usual multicasting only uses the knowledge of the network topology and the shortest paths. In consequence, multicast TE uses the same as unicast TE to achieve efficient network resource utilization.

Multicast traffic has some specific characteristics due to the multicast routing protocols nature [START_REF] Ooms | Framework for IP multicast in MPLS[END_REF]. Some of the multicast routing protocols are based on reverse path forwarding (RPF) to setup forwarding states on intermediate routers between the source and the destinations. But RPF is based on the idea that paths are symmetric in the network. When routing constraints are introduced, there is no guarantee that the link utilization is symmetric. Hence, RPF will cause forwarding on a sub-optimal path (in QoS routing) or might even prevent receivers from receiving traffic from some (or all) sources (in policy routing). This check must turned off or the multicast routing protocol must be able to obtain the constraint RPF via a constraint based routing (CBR) API. Multicast trees should be constructed taking into consideration the dynamism in the receiver set, and the receiver's heterogeneity, that is, receivers with different service requirements in terms of delay or jitter. This will introduce important modifications to CBR and conventional multicast routing protocols. Indeed, fast recovery for paths failure is very important in multicast TE since this failure may influence all the tree and not only the link in failure. Otherwise, load balancing should be carefully used since a packet should not pass by the same link more than one time. Finally, multicast forwarding is done based on the multicast IP address and that's why it is very difficult to aggregate multicast traffic since receivers can be located anywhere in the Internet.

The multicast traffic engineering trees can be built by expanding the existing protocols. There are two categories of protocols depending on the tree setup:

• Sender initiated tree setup: this kind of tree can have limited number of receivers with very rare join and prune action. Multicast trees are computed by the first-hop router from the source (root), based on sender traffic advertisements.

• Receiver initiated tree setup: this kind of tree can have a large number of receivers and they join and prune quite frequently. Multicast trees are computed from receivers to the root. Each receiver-side router independently computes a QoS-accommodating path from the source, based on the receiver reservation. This path can be computed based on unicast routing information only, or with additional multicast flow-specific state information. In any case, multicast path computation is broken up into multiple, concurrent unicast path computations. Finally, MPLS label switching can be used to forward unicast traffic through explicit routes and multicast traffic down the explicit tree to avoid RPF checking.

MPLS shows several advantages over conventional network layer forwarding [START_REF] Xiao | Traffic engineering with MPLS in the Internet[END_REF], [START_REF] Rosen | Multiprotocol label switching architecture[END_REF], [START_REF] Boudani | A simulator for multicast routing over an mpls network[END_REF]. Focusing on the advantages of the layer two switching protocol over , Multicasting over MPLS networks can benefit from the multicast reduce of traffic on one hand, and MPLS flexibility, speed and quality of service on the other hand.

III. MPLS PROPOSALS FOR MULTICAST TE

IP multicast protocols have different characteristics (scalability, computational complexity, latency, control message overhead, tree type, etc...). A framework for IP multicast deployment in an MPLS environment is proposed in [START_REF] Ooms | Framework for IP multicast in MPLS[END_REF]. Issues arising when MPLS techniques are applied to IP multicast are overviewed. Following characteristics are considered: aggregation, flood and prune, co-existence of source and shared trees, uni/bidirectional shared trees, encapsulated multicast data and loop free ness, and RPF check. The pros and cons of existing IP multicast routing protocols in the context of MPLS are described and the relation to the different trigger methods and label distribution modes are discussed. The framework did not lead to the selection of one superior multicast routing protocol but it concluded that different IP multicast routing protocols could be deployed in the Internet.

Using PIM-SM [START_REF] Estrin | Protocol Independent Multicast-Sparse Mode (PIM-SM): Protocol Specification[END_REF] join messages to distribute MPLS labels for multicast routes is proposed in [START_REF] Farinacci | Using PIM to distribute MPLS labels for multicast routes[END_REF] (called hereinafter PIM-MPLS). A piggy-backing methodology is suggested to assign and distribute labels for multicast traffic for sparse-mode trees. The PIM-SM join message is expanded to carry an MPLS label allocated by the downstream LSR. MPLS is not used with all its efficiency as a TE tool since the multicast tree still constructed using the RPF tree checking without constraints. In [START_REF] Boudani | A simulator for multicast routing over an mpls network[END_REF], we proposed a simulator for this methodology by using the MNS [START_REF] Chun | Overview of MPLS network simulator: Design and implementation[END_REF] (MPLS network simulator). We think that the join message in PIM-SM should be expanded to carry the explicited routed path towards the RP. A PIM-SM router always sends join/prune towards the upstream router listed in the explicited routed path. It can also carry other constraints, such as color or bandwidth. A new message join-nak can be sent from upstream to downstream if the upstream can not satisfy the constraints listed in the join message.

In [START_REF] Acharya | Ip multicast support in mpls networks[END_REF], authors consider the problem of supporting Ip multicast efficiently within MPLS environment for both PIM dense mode and sparse mode. They suggest a data-driven, per source assignment of labels to traffic on the shared tree and they present a common scheme for implicitly distributing and binding labels to multicast FECs. Authors suppose also like the previous proposal that multicast trees will be constructed using the RPF tree checking without constraints.

In [START_REF] Lee | Engineering paths for multicast traffic[END_REF], authors propose to engineer paths for IP multicast traffic in a network by directing the control messages to setup multicast trees on engineered paths. This proposal partitions the multicast traffic engineering problem such that multicast routing protocols do not have to be modified to allocate resources for multicast traffic nor do resource allocation protocols such as RSVP or CR-LDP have to be able to setup forwarding states (in this case labels) like multicast routing protocols. Resources are allocated on the same trip that paths are selected and setup. This prevent the problem of data being forwarded on branches of the tree where resources have not being allocated yet. An important aspect of this proposal is that it enables multicast paths to be engineered in an aggregatable manner, allowing this solution to scale in the backbone. But while this proposal uses MPLS (label and explicit route object) to cause engineered paths to be selected, it forwards data using multicast routing.

Another interesting proposal is aggregated multicast [START_REF] Fei | Aggregated multicast: An approach to reduce multicast state[END_REF]. The key idea of aggregated mulicast is that, instead of constructing a tree for each individual multicast session in the core network, one can have multiple multicast sessions share a single aggregated tree to reduce multicast state and, correspondingly, tree maintenance overhead at network core. In this proposal we address two requirements: (1) original group addresses of data packets must be preserved somewhere and can be recovered by egress nodes to determine how to further forward these packets; (2) some kind of identification for the aggregated tree used by the group must be carried by the packets and transit nodes must forward packets based on this identification. In group to aggregated tree matching, complication arises when there is no perfect match or no existing tree covers a group (leaky matching). The disadvantage in leaky matching is that certain bandwidth is wasted to deliver data to nodes that are not involved for the group. Bandwidth can be a crucial factor for provisioning QoS in multicast networks and even for best effort Internet. To handle aggregated tree management and matching between multicast groups and aggregated trees, a centralized management entity called tree manager is introduced.

In [START_REF] Ooms | Mpls multicast traffic engineering[END_REF], extensions to CR-LDP are proposed to construct multicast trees immediately on L2. Thus the mapping of L3 trees onto L2, as described in [START_REF] Estrin | Protocol Independent Multicast-Sparse Mode (PIM-SM): Protocol Specification[END_REF] and [START_REF] Acharya | Ip multicast support in mpls networks[END_REF] is not needed. All of the descriptive parameters of the tree must be carried in the initial label request. Given this and given that it is highly undesirable to fragment such requests, the tree building process is primarily applicable to trees with a small number of receivers.

In [START_REF] Chung | Ldp extensions for mpls multicasting services[END_REF], [START_REF] Chung | Rsvp-te extensions for mpls multicasting services[END_REF], extensions to LDP and RSVP for MPLS multicasting services are proposed. In these two proposals multicasting functions of LDP and RSVP are independent of traditional IP-based multicast routing protocols (such as DVMRP, MOSPF, PIM, etc.) and multicast trees are calculated by a special entity.

To enable MPLS based multicasting, the tree building with join, leave, destroy and RPF messages should be directly implemented in LDP and RSVP. New messages and extension of existing messages are studied for insertion into LDP or RSVP protocols. Multicasting message (for join, leave and destroy operations) is created. Extensions to hello, notification, path messages, the label request, the label mapping and the multicast forwarding table are introduced. These two proposals require MPLS and multicast routing protocols to be merged, an exercise which tend to increase the complexity of multicast traffic engineering while not providing any means of aggregation of multicast traffic.

The complete tree information should be stored in all LSR-RP (branching nodes in the tree). Multicast hello messages are used to inform the LSRs of the multicasting source and group IP address of the multicasting tree. When the number of group grows the number of hello messages grows also. And since we will send for every source and group a notification message, when the number of group grows the number of notification messages grows also. It should be noted that a new table for multicast should be created independent of the existing unicast table. It is not very clear how the source will choose the LSR-RP.

Using MPLS with multicast has many benefits not only for reducing multicast forwarding states but also for traffic engineering and QoS issues. In this paper, we only focus on the scalability problem. We propose a novel approach that uses MPLS LSPs between multicast tree branching node routers in order to reduce forwarding states and enhance scalability.

In [START_REF] Boudani | A new approach to construct multicast trees in mpls networks[END_REF], we proposed a new approach to construct multicast trees in MPLS networks. Each domain contains a network information manager system (NIMS) for each group, charged to collect join and leave messages from all group members in that domain. The NIMS is elected through a mechanism similar to the one used to elect the Rendez-vous Point router in PIM-SM [START_REF] Estrin | Protocol Independent Multicast-Sparse Mode (PIM-SM): Protocol Specification[END_REF]. After having collected all join messages, the NIMS computes the multicast tree for that group in the domain (using the short path tree algorithm). Branching nodes, for any group, could be computed from the network graph. The NIMS sends then branch messages to all branching node routers to inform them about their next hop branching node routers. On receiving this message, a branching node router creates a multicast forwarding state for the multicast session. Packets will be sent from a branching node router to another until arriving at their destinations.

Already established MPLS LSPs are used between multicast tree branching node routers in order to reduce forwarding states and enhance scalability. When a multicast packet arrives to the ingress router of an MPLS domain, the packet is analyzed according to its multicast IP header. The router determines who are the next hop branching node routers for that packet. Based on this information, multiple copies of the packets are generated and an MPLS label is pushed into the multicast packet according to next hop branching node router. When arriving to a next hop branching node router, the label is pulled up and again the same process is repeated. This process should be repeated until the packet arrives to its destination (see Fig. 1). Only those routers acting as branching nodes for a group need to keep forwarding state for that group. Other routers between two branching nodes do not need to store multicast states. Unicast LSP is used between two branching node routers. This way the total number of multicast forwarding states may be significantly reduced. In our approach the same LSP could be used for multicast and for unicast traffic. Other approaches use different labels for multicast and unicast traffic in consequence they require specific encoding techniques and additional overheads. When arriving to a LAN, the packet un-labeled can be delivered by conventional multicast protocols according to IGMP [START_REF] Cain | Internet group management protocol, version 3[END_REF] informations.

IV. MPLS EXTENSION TO THE NETWORK SIMULATOR

MPLS extension to the network simulator NS does not work with multicast routing, particularly because (1) there is no label setup mechanism for multicast groups, (2) there is no multicast replicator to cooperate with MPLS classifier, and (3) MPLS header contains pointers, which do not work with multicast replicator. In this section, we describe the modifications needed to allow multicast packet transmission in MPLS networks without implementing a new protocol. Two main points are to be considered: information tables of MPLS nodes, and multicast packet transmission. Our major objectif was implementing MMT without major modifications of the unicast MPLS code already existing in NS.

A. Information tables of MPLS nodes

As mentioned in [START_REF] Farinacci | Using PIM to distribute MPLS labels for multicast routes[END_REF], an MPLS node contains three information tables: LIB (Label Information Base), PFT (Partial Forwarding Table ), and ERB (Explicit Routing information Base). To apply the MMT proposal, a mapping of the (S, G) session to more than one FEC on one hand, and a mapping of each FEC to one <incoming label, incoming interface> and thus to one <outgoing label, outgoing interface>, on the other hand, are needed. The information base at the MPLS nodes must be modified as shown in Fig. 2. For the first mapping, the SGF (Source Group FEC) table is defined. This table exists only in branching node routers and includes three fields: Source, Group, and FEC. This table is filled at each node after receiving the branch message sent by the NIMS. In branching node routers, more than one FEC could be associated to one session ¡S,G¿ packet, because the router receiving the packet could have several next-hop branching routers. In consequence a complete searching of all the entries of the SGF table should be done. Each FEC is mapped to a couple <incoming label, incoming interface> in the LIB table. For the second mapping, each <incoming label, incoming interface> is mapped to exactly one <outgoing label, outgoing interface>. The LIB table remains unchanged with one <incoming label, incoming interface> for one <outgoing label, outgoing interface>.

B. Multicast packet transmission

Data is processed exactly as in unicast MPLS packets with only one difference at branching nodes. In branching node packets are processed as follows: When a labeled packet arrives, a swap/pop operation is executed and the LIB table is examined. If as a result of a pop operation, the packet remains without label, a global search is done in the SGF table to attribute the packet to FECs and thus to <incoming label, incoming interface> couples. For each FEC, a packet copy is created, and then the incoming label is swapped with the corresponding outgoing label, and then transmitted to the outgoing interface.

V. EVALUATION AND SIMULATION

MMT was evaluated in [START_REF] Boudani | A new approach to construct multicast trees in mpls networks[END_REF] in terms of scalability (state and control messages overhead) and efficiency (tree cost and data processing). The state information requirement can be measured using the average multicast forwarding table size. The control messages overhead can be measured in terms of average number of control messages sent per link or the total percentage of bandwidth spent on control traffic. MMT allows only the shortest path trees, which are the most efficient for data forwarding. Besides, since we are using label switching at routers, our approach may be considered more efficient in data movement than other schemes.

Multicast address aggregation is important since multicast groups may share some links in their multicast trees. In conventional multicast, it is not possible to aggregate multicast IP addresses. Receivers can be located anywhere in the Internet, there is no other alternative than having one entry by multicast IP address in the multicast routing table. Since in our approach, we are using MPLS, the aggregation problem of multicast IP addresses can be transformed to a simple aggregation of labels.

A. Simulation Analysis

We simulate MMT in NS (Network Simulator) [START_REF] Fall | The NS Manual[END_REF] to validate the behavior of our approach and its effectiveness in state reduction and tree construction. The performance of MMT is compared to PIM-MPLS. PIM-MPLS in our simulations refers to the simulator described in [START_REF] Boudani | A simulator for multicast routing over an mpls network[END_REF]. In [START_REF] Boudani | A simulator for multicast routing over an mpls network[END_REF] we presented a simulator for multicast routing over an MPLS network where we choosed PIM-SM (source specific tree) as the multicast routing protocol. In this paper we present the MMT protocol simulation which will be compared to PIM-MPLS.

B. Simulation Scenario

We use in our simulation two network models generated by the GT-ITM generator [START_REF] Zegura | How to model an internetwork[END_REF]: each model has a graph of 100 nodes and all the links in the network are bidirectional links with 20Mbps bandwidth. The topology of the first model is generated by the first Waxman algorithm [19] and used as an example of a dense network (with 0.3 as the node degree distribution). The topology of the second model is generated by a pure random algorithm, it has 5 domains and is an example of a sparse network. Four domains contain receivers and sources only, while the fifth domain is considered as the core domain. The horizontal axis is the percentage of sources that are active in the network, and the vertical axis is the overall forwarding table size in the network. The poly-lines labeled PIM-x and MMT-x show the overall forwarding table size for PIM-MPLS and MMT protocols respectively when the number of receivers per group is x.

The forwarding table size grows with the number of active groups and the number of receivers. From Fig. 4 and Fig. 3 we can see that the relative state reduction of MMT is roughly 40% and 80% respectively compared to PIM-MPLS. We deduce also that our protocol is more suitable for sparse mode networks and for groups with few members.

VI. CONCLUSION AND FUTURE WORKS

In this paper, we have presented a framework for multicast traffic engineering. We have defined multicast traffic engineering and we have studied its particularity comparing to unicast traffic engineering. We have presented several MPLS proposals for multicast TE. We have described our approach, the MPLS multicast tree protocol which utilizes MPLS LSPs between multicast tree branching node routers in order to reduce forwarding states and enhance scalability. We have discussed the realization of a simulator for MMT and finally we have presented some simulation results.

Fig. 1 .

 1 Fig. 1. The MPLS multicast tree construction

Fig. 2 .

 2 Fig. 2. Structure of tables for MMT packet switching

[

  19] B. Waxman. Routing of multipoint connections. IEEE Journal on Selected Areas in Communications, 6(9):1617-1622, December 1988. Ali Boudani received his B. Eng degree from the Lebanese University, in 1994. He is currently a PhD student in Computer Engineering at the University of Rennes1 and he is working in the Labs of the IR-ISA/INRIA Rennes. His research interrests include multicast, QoS and traffic engineering. Bernard Cousin is professor at the University of Rennes 1 (France). He is member of a research group on networking at IRISA (CNRS-INRIA-University joint reseach laboratory in computing science located in Rennes). His main research topics include high speed networking, trafic engineering, multicast routing, QoS management and network security. Jean-Marie Bonnin received his PhD in 1998 at the university of Strasbourg (France). His PhD work was on large scale reliable multicast mechanisms. He is now associate professor at the ENST bretagne (France). He is member of the INRIA's ARMOR team, a research group specialized in computer science in Rennes. His main research topics include MPLS and trafic engineering, as well as the QoS management in 3G IP core networks.

TABLE I SUMMARY

 I Nr receivers are randomly deployed in the network graph. A receiver join randomly a group. TableIsummarizes the parameters used in the simulation.The forwarding table size in all routers in the network using the pure random sparse model is shown in Fig.3and using the waxman model is shown in Fig.4.

					OF SIMULATION PARAMETERS
	N	100			number of nodes in the network
	NT	10, 20, 30,	percentage of sources in the net-
		40, 50, 60	work (number of trees)
	Nr	3, 6, 9, 12,	number of receivers for each source
		15, 18
	500 1000 1500 2000 2500 3000 NT sources and 0 Forwarding table size PIM-3 MMT-3 PIM-6 MMT-6 PIM-9 MMT-9 PIM-12 MMT-12 PIM-15 MMT-15 PIM-18 MMT-18
		10	15	20	25	30	35	40	45	50	55	60
				Percentage of sources in the network

Fig. 3. Forwarding table size -pure random sparse mode model


V Lebastard

C Chevallereau

A Girin

F Boyer

P B Gossiaux

Localization of small objects with electric sense based on kalman filter

Electric fish feel the perturbations of a selfgenerated electric field through their electro-receptive skin. This sense allows them to navigate and reconstruct their environment in conditions where vision and sonar cannot work. In this article, we use a sensor bio-inspired from this active sense to address the problem of small objects reconstruction and electrolocation. Based on a Kalman filter, any small object in the surrounding of the motion controlled sensor can be encapsulated into an equivalent sphere whose location is well estimated by the filter. Experimental results illustrate the feasibility of the approach.

I. INTRODUCTION

Developed by several hundreds of fish who have coevolved on both African and South-American continents, electric sense have been discovered by Lissman in 1958 [START_REF] Lissmann | The mechanism of object location in gymnarchus niloticus and similar fish[END_REF]. In the African fish Gnathonemus Petersii for instance, the fish first polarizes its body with respect to an electric organ discharge (EOD) located at the basis of its tail (Fig. 1). This polarization which is applied during a short timepulse, generates around the fish a dipolar shaped electric field which is distorted by the objects present in its surroundings. Then, thanks to many electro-receptors distributed along its body, the fish "measures" the distortion of the electric fields and infers through a brain processing intensively studied by neurobiologists [START_REF] Caputi | The electric image in weakly electric fish. ii physical images of resistive objects in gnathonemus petersii[END_REF] an image of its surroundings. In order to study the environment modalities sensed by the fish, biologists have developed many neuro-ethological experiments in which fish are trained to recognize a given shape [START_REF] Der Emde | Electric fish measure distance in the dark[END_REF], [START_REF]Electric fish measure distance in the dark[END_REF]. More commonly, electric fish can easily navigate in the dark or/and turbid waters of confined unstructured environments as the roots of the trees of flooded tropical forests which are their natural habitat. Also, electric sense has a quite narrow but relevant niche since none other sense as vision or sonar can work in these conditions. Based on these potential interests, Mc. Iver and co authors have recently exploited an experimental set up of four point electrodes placed at the apexes of a rhombus in a rigid moving frame driven by a Cartesian robot [START_REF] Solberg | Robotic electrolocation: Active underwater target localization[END_REF]. In this device, two electrodes situated at the opposite apexes of the lozenge are polarized in voltage and play the role of the electric organ of discharge of the fish while the two opposite electrodes play the role of receivers V. Lebastard is the corresponding author: E-Mail: vincent.lebastard@emn.fr, Tel.: +00 33 2 51 85 83 08.

* F. Boyer, A. Girin and V. Lebastard are with IRCCyN-Ecole des Mines de Nantes, 4, rue Alfred Kastler B.P. 20722 -44307 Nantes Cedex 3 France.

⋆ C. Chevallereau is with IRCCyN-CNRS, Ecole Centrale de Nantes, 1 rue de la Noë, 44000 Nantes, France.

† P.B. Gossiaux is with SUBATECH-Ecole des Mines de Nantes, 4, rue Alfred Kastler B.P. 20722 -44307 Nantes Cedex 3. between which the tension is measured. Using the Cartesian motion control of this sensor in a tank, they successfully located small spheres through off-line particle filtering [START_REF]Active electrolocation for underwater target localization[END_REF], [START_REF] Solberg | Robotic electrolocation: Active underwater target localization[END_REF]. However, the body geometry of this device ideally reduces to four points electrodes so limiting the implementation of its principle on a real autonomous underwater vehicle. Thus, using an alternative sensor technology based on the measurements of the currents, in [START_REF] Baffet | Biomimetic localization using the electrolocation sense of the electric fish[END_REF], [START_REF] Lebastard | Underwater robot navigation around a sphere using electrolocation sense and kalman filter[END_REF], Kalman filtering technics have been applied to the reconstruction of the surrounding of an underwater vehicle of realistic geometry. In [START_REF] Baffet | Biomimetic localization using the electrolocation sense of the electric fish[END_REF] the surroundings is constituted of insulating walls and corners as those of a tank, while in [START_REF] Lebastard | Underwater robot navigation around a sphere using electrolocation sense and kalman filter[END_REF], the case of small (insulating and conducting) objects is tackled. However, all these results have been till today restricted to simulations. Also, the article here presented aims at experimenting these solutions on a electric test-bed. Going further, the article here presented mainly addresses, in experimental conditions, the problem of the reconstruction of small objects. By "reconstruction", we mean the estimate of their size and their location in the sensor frame. Because the equations of electricity cannot be integrated analytically for an object of arbitrary shape [START_REF] Sihvola | Polarizabilities of platonic solids[END_REF], we reduce our investigations to the case of spherical objects. However, we will see that other shaped objects as the cube, can be modeled by (and encapsulated in) equivalent spheres. The article is structured as follows. First we will briefly present the sensor technology (section II). Then, we address in section III the modeling of the sensor measurements generated by its surrounding. A solution to localize and reconstruct the small objects is detailed in section IV. It will be tested on an experimental test-bed in section V. Finally, the article ends with some concluding remarks and perspectives in section VI.

II. THE SENSOR

Here, we just remind the physical principle of the sensor and refer the reader to [START_REF] Servagent | Bioinspired sensor for electrolocation and navigation in conducting media[END_REF] for more details about its technology.

A. The sensor principle

From an external point of view, the sensor is constituted of a body with insulating wetted boundaries on which are fixed a set of m + 1 electrodes noted E 0 , E 1 ...E m (figure 2). Thanks to a tension generator, the potential on E α=1,2...m is imposed to be the common ground with respect to which E 0 is put under the controlled tension U . Once the sensor immersed in a fluid, this active device produces a field of current lines flowing through the sensor surroundings from E 0 to the other macro-electrodes E α=1,...,m . These current lines are closed through an electric circuit embarked inside the sensor. This circuit measures the vector of the currents I = (I 1 , I 2 , ..., I m ) flowing across the E α=1,2,...,m . When an object is placed around the sensor, it generates perturbations of the sensor electric lines portrait that are measured by this current measurement circuit. Finally, the whole device is equipped with an electronic board able to change (in real time) the electric connectivity of the electrodes in order to give to any of the electrodes, the role of E 0 in all what has been just said. This ability which we name "electric reconfigurability" allows one to enrich the measurement by changing the previous scalar control input U into a m × 1 vector U. The entire device can be embarked on a mobile body as typically an underwater motion controlled vehicle.

B. Sensor morphology

Based on this principle, we built a first generation of sensors named "slender probes" due to their high aspect ratio (length/thickness) morphology. On these probes, the macroelectrodes are rings or hemispheres as it is illustrated on the example of Fig. 3 (which is named the 4-electrode probe). Based on these design, the sensor has a range of about its body length. This range is that observed on fish of the mormyride family for instance. This limit defines the surroundings of the sensor beyond which the sensor cannot perceive any change in its environment.

III. ANALYTICAL MODEL OF MEASURE

The sensor is approximated by an array of m+ 1 spherical conducting electrodes aligned on a rigid line materializing its solid axis. In the following, the characteristics of the tested sensor are those of Figure 3 with R = 1cm for its radius, and L = 22cm for its total length. The conductivity of the water is γ = 0.04S/m (which corresponds to the ordinary tap water of our lab). By imposing a difference of potentials between each electrode and all the others, an electric field is produced around the sensor. This field produces a field of current lines flowing across the measurements electrodes. The corresponding measured currents are influenced by the presence of objects in the sensor surroundings (Fig. 45). In all cases, the effects of the electrically un-transparent objects onto the sensor measurements can be modeled by an expression of the general form (can be numerically done by boundary elements method in [START_REF] Porez | Multi-physic model of an electro fish-like robot : numerical aspect and application to the obstacle avoidance[END_REF]):

U = R(x)I, ( 1
)
where U is the m × 1 vector of the imposed potential differences, I is the m × 1 vector of the measured currents, R is a (m×m) resistance matrix of the scene, while x is the state vector gathering the (electric and geometric) parameters of the scene. For example, based on Fig. 4, the sphere location with respect to the sensor is parameterized by a distance d and an angle θ, while the parameters encoding the intrinsic properties of the sphere are its electric conductivity γ s and its radius a. As a result, for this scene the state vector of geometric parameters is: x = (d, θ, a) T . Now reconsidering (1), the measured currents I are given by

I = R -1 (x, γ, γ s )U (2)
which is highly non linear with respect to the state. Going further, in the following we shall consider the sensor moving in a cubic tank bounded by insulating walls in which are placed a small sphere. In these conditions, the resistance R offered to the sensor takes the most general form :

R = A + S + W x + W y + C, (3)
where:

• A models the sensor self resistance, which only depends on the geometry of the robot and is detailed in [START_REF] Lebastard | Underwater robot navigation around a sphere using electrolocation sense and kalman filter[END_REF], [START_REF] Jawad | Sensor model for the navigation of underwater vehicles by the electric sense[END_REF]. • S models the sphere influence which depends of d, θ, a, γ, γ s and is based on a model initially proposed in [START_REF] Rasnow | The effects of simple objects on the electric field of apteronotus[END_REF]. A further study of this influence is developed in [START_REF] Lebastard | Underwater robot navigation around a sphere using electrolocation sense and kalman filter[END_REF]. • W x (W y ) models the wall resistance. It depends on d x and θ c (d y and θ c ). A study of this influence based on the so-called "method of images", is developed in [START_REF] Jawad | Sensor model for the navigation of underwater vehicles by the electric sense[END_REF]. • C models the influence of the corner which depends on d x , d y and θ c as detailed in [START_REF] Jawad | Sensor model for the navigation of underwater vehicles by the electric sense[END_REF].

Finally, the model of measurements take the general form:

y = I = h(x, γ, γ s ), (4)
where x is the geometric state vector of the sensor surroundings. For the sake of simplicity, in all the following, we assume that: 1 • ) the objects are motion less in the tank frame, 2 • ) the sensor moves in its equatorial plane which is also a symmetry plane (top-down) for the objects; 3 • ) γ and γ s are a priori known. Thus, the reconstruction only concerns the geometric part of the state x, or "geometric state".

IV. KALMAN FILTER

The sensor is put into a tank of our test bench (section V-A) and its (planar) motion is imposed in two kinds of scenarios. In the first set of scenarios, the scene around the sensor is composed of a small object (sphere or cube) in an unbounded conductive fluid (Fig. 4) in which the sensor (or "robot") goes forward in straight line and close to the sphere. In the second type of scenarios, the scene is composed of four walls disposed in a square and of a small object (Fig. 5). In this environment, the robot follows the wall. Finally, in these two situations, the objective of the observer consists in estimating the geometric parameters of the environment gathered into the (geometric) state vector x.

A. Control and Robot Evolution model

In our study, the motion of the sensor is imposed. In order to mimic the fish locomotion with a simple model (the fish cannot swim along the lateral dimensions of its body), the sensor motion is modeled as that of a non-holonomic unicycle whose the motions parallel to the y axis are forced to zero. Then the velocities V (axial velocity) and ω (yawing velocity) define the vector of control inputs u = (V , ω) T . In the scenario of Fig. 5, the sphere and the corner (defined by the walls W x , W y and their intersection C) being fixed in space, the model of evolution of the scene with respect to the sensor can be written under the form:

                 d k+1 x = d k x -V k cos(θ k c )∆t d k+1 y = d k y -V k sin(θ k c )∆t θ k+1 c = θ k c + ω k ∆t d k+1 = d k -V k cos(θ k )∆t θ k+1 = θ k + (ω k + V k sin(θ k ) d k )∆t a k+1 = a k (5)
where ∆t is the sampling period, k denotes the time step at which the the state vector,the state vector of geometric variables is then

x k = (d k x , d k y , θ k c , d k , θ k , a k ) T .
In the first type of scenario (Fig. 4), the robot goes forward with an axial velocity V = V d > 0, while its yawing velocity is forced to zero (ω = 0). In any case, the evolution model can be written in the general form :

x k+1 = f (x k , u k ), (6)
where x k is the geometric state vector of the scene and the vector of control inputs is

u k = (V k , ω k ) T .

B. Kalman Filter

The objective of the observer is to estimate the state vector x of each scene. The chosen observer is an unscented Kalman filter [START_REF] Simon | A new extension of the kalman filter to nonlinear systems[END_REF]. This choice enjoys the advantages of the classical Kalman filter [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF], [START_REF] Nijmeijer | Nonlinear Dynamical Control Systems[END_REF] while avoiding the complex analytical linearization of the evolution equation ( 6) and mesured model 4. As any Kalman filter, our filter requires to define some matrices which tune its properties. The first one is the covariance matrix P k which is attached to the state variable x k and is a trust indicator associated to the estimated state, a small (large) covariance meaning that the state is precisely (imprecisely) known. The evolution of P k is defined by the unscented Kalman filter algorithm. Two other matrices are defined: Q m , which defines whether the model ( 6) is well known or not, and Q s which defines the level of noise on the measurements (4). The principle of the unscented Kalman filter is reminded on Fig. 6. Based on its covariance P k , the filter computes some sampling points (called "sigma points" χ k+1 ) around the estimated state x k using a deterministic sampling technique known as unscented transfrom. Then, it propagates the sigma points χ according to: 1 • ) the evolution model

χ k+1 = f (χ k , u k ), 2 •
) the model of measurements ỹ = h(χ k , γ, γ s ), which improves the convergence of the estimated state xk+1 and the covariance P k+1 .

V. EXPERIMENTATIONS

In order to illustrate the performances of the algorithms of object reconstruction and localization, we now report three experiments which have been carried out in a test bench presented in section V-A. Two of these experiments are carried out in the conditions of Fig. 4, first with an insulating sphere, and second, with an insulating cube. In these two experiments, the small object is modeled by an equivalent sphere whose size and location is estimated. The third (and last) experiment is carried out in the scenario of Fig. 5, where an insulating sphere is placed in a corner while the robot follows the walls of the tank. All, these experiments are carried out with the four-electrode sensor of Fig. 3. The third experiment focuses on the ability of the algorithm to extract the parameters of a small object, once immersed in a complex scene. For each of the experiments, we record the velocity (V and ω) and the measured currents I = (I 1 , I 2 , I 3 ) T .

A. Tank and cartesian robot

In order to test our electrolocation sensors and algorithms in controlled and repeatable conditions, an automated test bench consisting in a tank of one cubic meter side and a three-axis cartesian robot has been built (see Fig. 7). The robot fixed on top of the aquarium allows probes positioning in translation along X and Y with a precision of 1/10mm and the maximum speed available is 300mm/s (≃ 1km/h) for both translations. The orientation in the (X, Y ) plane is adjusted in 0.023 • using an absolute yaw-rotation stage and the velocity is 80 • /s (13.5tr/min). The three axes robot is motion-controlled using simulink software with the dSpace system. All probes tested are positioned in the aquarium at adjustable height (Z) using a rigid glass epoxy fibre tube. This vertical insulating tube allows the passage of electrical cables dedicated to the signals coming from the electrodes of the probe. The electrolocation signals first processed by analogue chain (amplification and filtering) then converted by a 16 bits ADC (DS2004 card) with a resolution of 0.3mV/bit (range +/ -10V ).

B. Reconstruction of a small object

In the following, any small object in the surroundings of the sensor is modeled as a sphere whose the size and location (related to the sensor) are reconstructed by the filter. Before testing the sphere localization, two preliminary tests have been carried out in order to assess whether the modeling assumptions (as replacing any object by a sphere for instance) are justified. All x (m)

I1 (A) (b)
The probe approaches a insulating cube.

Fig. 8. Measured (red) and modeled currents (black) with the real state vector x when the probe approaches a small object (sphere model).

comparing the actual measurements with those given by the model ( 4), once it is informed by the actual geometric (and electric) parameters of the scene. In the first test, a sphere is placed in the middle of the tank and the measurements are recorded in the conditions of Fig. 4. Then, the measurements are computed from the model ( 4) with the real extended state (x, γ, γ s ), and for a scene approximated by the sphere immersed in an unbounded fluid domain. In the second test, the measurements are carried out with a cube located in the middle of the tank, and compared with the same model as in the previous test (a sphere alone in an unbounded world).

Here, the sensor goes in straight line with V = 0.04m/s and ω = 0. The figures 8(a)-8(b) show the 3 measured currents of the sensor (red) as well as the calculated currents (black) (deduced from the model), when the probe approaches the object. In the conditions of the results of Fig. 8(a), the sphere was insulating and with a radius a = 0.0305m. The results show that the measured and calculated currents are enough close to each other to consider the model of the sphere enough good for its reconstruction. The plots or figure 8(b) have been obtained with a real insulating cube (of side-length a = 0.040m) located in the middle of the tank, and with the model of a sphere of radius a = 0.037. The measured and calculated currents are close to each other, so allowing us to conclude that the approximation of a cube by an equivalent sphere is justified. 1) Reconstruction of a sphere in a simple scene (Fig. 4): By "simple scene", we mean that the object (here a sphere) is alone in the sensor surroundings. In this first test devoted to the state reconstruction, the robot moves forward in the tank where an insulating sphere with a radius a = 0.0305m is immersed. The initial real and estimated states of the scene are fixed to : (d = 0.257 m, θ = 0.294 rad, a = 0.0305 m) and: d = 0.146 m, θ = 0.494 rad, â = 0.0105 m). Figure 9(a) shows time evolution of the real (red) and estimated (blue) states of the sphere. Figure 9(b) represents the real (cyan) and the estimated (blue) scenes at different times of the motion. At time t = 2s the sphere is quite good localized but its estimated size is bigger than the actual one. Between 2s < t < 8s, the estimated sphere is close to the real one. When the object is close to the emitter (t = 8s), the real sphere is well reconstructed and the estimated error is small. 2) Reconstruction of a cube in a simple scene (Fig. 4): This second test is carried out in the same conditions as the previous one, except that the object is now an insulating cube. Based on a sphere model, the filter so estimates the parameters of an equivalent sphere. The initial real and estimated states of the scene are fixed to: (d = 0.225 m, θ = 0.298 rad, a = 0.037 m (equivalent radius)) and: ( d = 0.145 m, θ = 0.498 rad, â = 0.001 m). Figure 10(a) shows the time evolutions of the real (red) and estimated (blue) state of the cube. Figure 10(b) represents the actual (cyan) and the estimated (blue) scenes at different times of the sensor motion. Like the previous test, at time t = 2s the cube is quite good localized but its estimated size is bigger than the side-length of the cube. When the object is close to the emitter (t = 8s), the equivalent sphere is well reconstructed and the estimated error is small. This test illustrates a common feature of all the tests carried out: the estimated equivalent sphere always encapsulates the cube, while its location is well estimated.

C. Localization and reconstruction of a sphere in a complex scene.

This test is carried out in the conditions of Fig. 5, with an insulating sphere localized in a corner while the robot follows the wall. This test focuses on the ability of the filter to extract the parameters of a small object in a complex scene. The sphere has a radius a = 0.0305m. The initial real state of the scene is: (d x = 0.3 m, d y = 0.3 m, θ w = 0 rad, d = 0.4 m, θ = 1.57 rad, a = 0.034 m) and the initial For the time between 0 < t < 8s, the sensor doesn't detect the sphere because it is too far. At t > 8s the sensor perceives the perturbation due to the sphere, then the estimated parameters of the sphere converge to the real ones. For the time 8 < t < 23s, the sensor is close to the sphere and for t > 20s the sphere is quite good localized. This test highlights that it is possible to extract (localize and reconstruct) a small object plunged in a complex scene, at least when the sensor is close to this object.

VI. CONCLUDING REMARKS AND PERSPECTIVES

Using a sensor bio-inspired from electric fish, we shown through experiments: first, that small objects can be reconstructed and localized as equivalent spheres; second, that small objects can be extracted from a complex scene. In spite of these encouraging results, many further works should have to be pursued in future. For instance, let us remark that with the Kalman filter, the current differences between a sphere and a cube are too small to allow the filter to distinguish between them. As a consequence, other strategies more oriented towards the haptic modality of the electric sense could be implemented in order to disambiguate the perception in this critical case. Another way of improvement of the presented result concerns the richness of the sensor measurements. Indeed, with the current sensor design the reconstruction of a sphere (or cube) suffers from intrinsic ambiguity due to the symmetry of the sensor. For instance, two identical objects symmetrically placed on both sides of the sensor generate the same measurements. This scenario is illustrated on Fig. 12, where the estimated state of the scene is initialized to a symmetric position so leading to a reconstructed sphere, symmetric of the real sphere with respect to the sensor axis. In order to disambiguate this kind of situation, we developed a sensor in which electrodes are separated into two lateral (left and right) measurement sub-electrodes. For such a more sophisticated sensor design, the model of the measurements of small objects has been recently proposed in [START_REF] Boyer | Model for a sensor bio-inspired from electric fish[END_REF]. Among its advantages, this design also allows to identify the electric nature -insulating or conducting -of the objects material, so allowing to release the a priori knowledge of γ and γ s . The implementation of these models in the unscented Kalman filters is today a work in progress. Finally, in the article here presented only the reconstruction problem is addressed. Also, in a second step we should have to include this filter in a global control feedback loop allowing a real underwater robot to navigate in a complex scene. However, we recently proposed another solution to the navigation problem. Based on the embodied philosophy, this alternative solution requires no model of the scene and considers the navigation problem as a low level reflex control problem. Thus, taking benefit of this recent result would allow one to directly couple the reconstructor here presented with this low-level reflex navigator.

VII. ACKNOWLEDGMENT

The ANGELS project [START_REF]The angels' project[END_REF] is funded by the European Commission, Information Society and Media, Future and

Fig. 1 .

 1 Fig.1. The electric fish, here the Gnathonemus Petersii. The electric field is distorted by the presence of an object (for instance, an insulating cylinder blocks the field lines whereas a conductive cube funnels them).

Fig. 2 .

 2 Fig. 2. General principle of the U -I sensor from the external (left) and internal (right) point of view.

Fig. 3 .

 3 Fig. 3. Picture (left) and schematic view (right) of a 4-electrode sensor.

Fig. 4 .

 4 Fig. 4. Simple scene: The sensor moves in an unbounded environment with a sphere, along the red dotted path.

Fig. 5 .

 5 Fig. 5. Complex scenes: The sensor follows the walls of a tank in which a sphere is placed, along the red dotted path.

Fig. 6 .

 6 Fig. 6. Principle of the observer.

Fig. 7 .

 7 Fig. 7. Electrolocation test bench.

Fig. 9 .

 9 Fig.9. Reconstruction of a sphere in a simple scene (V = 0.04m/s).

Fig. 10 .

 10 Fig.10. Reconstruction of a cube in a simple scene (V = 0.04m/s).

  Real (red) and estimated (blue) state of the sphere.

Fig. 11 .

 11 Fig.11. Reconstruction of a sphere in a complex scene (V = 0.04m/s).

Fig. 12 .

 12 Fig. 12. Real (cyan) and the estimated (blue) scene..

  these tests are carried out by

		6.05 x 10 -3								6.05 x 10 -3
	I3 (A)	5.95 6							I3 (A)	5.95 6
		-0.15	-0.1	-0.05	0	0.05	0.1	0.15	0.2	-0.15	-0.1	-0.05	0	0.05	0.1	0.15	0.2
					x (m)								x (m)
		6.18 x 10 -3								6.18 x 10 -3
		6.16								6.16
	I2 (A)	6.1 6.12 6.14							I2 (A)	6.1 6.12 6.14
		6.08								6.08
		-0.15	-0.1	-0.05	0	0.05	0.1	0.15	0.2	-0.15	-0.1	-0.05	0	0.05	0.1	0.15	0.2
		7.5 x 10 -3			x (m)					7.5 x 10 -3			x (m)
	I1 (A)	7.45								7.45
		7.4								7.4
		-0.15	-0.1	-0.05	0	0.05	0.1	0.15	0.2	-0.15	-0.1	-0.05	0	0.05	0.1	0.15	0.2
					x (m)
	(a) The probe approaches an insulat-
	ing sphere.


Stefano Mona
email: mona@mnhn.fr

On the role played by the carrying capacity and the ancestral population size during a range expansion

Keywords: Metapopulation, population structure, stepping stone, approximate Bayesian computation

Most species are structured and various population genetics models have been proposed to investigate their history. For mathematical tractability, most of these models make the simplifying assumption of equilibrium. Here we focus on the properties of a non-equilibrium spatial explicit model, range expansions (RE). Despite their abundance, many details of their genetic consequences need yet to be fully investigated. The model we studied is characterized by four main parameters: the effective population size of each deme (N), the migration rate per generation per deme (m), the time of the expansion (Texp) and the effective size of the deme from which the expansion started (Nanc). By means of extensive coalescent simulations, we focused on two aspects of range expansions for fixed Nm: i) the separate influence of N and m; ii) the role of Nanc. We compared our results to an equilibrium stepping stone model and found two main features typical of REs: an excess of rare variants for larger N and a complex interaction between N, Texp and Nanc in shaping the degree of population differentiation (which depends only on Nm in the stepping stone model). Finally, we developed an approximate Bayesian computation approach to jointly estimate N and m and to infer Nanc. When applied to pseudo-observed data sets, we could correctly recover both N and m (but not Nanc), provided a large number of demes were sampled. These findings highlight how it will be possible to estimate the dispersal rate in non-equilibrium metapopulations by using population genetics approaches.

Introduction

In recent years, advances in both theoretical and empirical population genetics have highlighted the importance of taking population structure into account when interpreting the variability of a species and inferring demographic parameters. Most (if not all) populations are indeed not isolated, but linked in a network of demes exchanging migrants to some extent (i.e., they belong to a metapopulation). Many recent works have shown that failing to take the structure of the metapopulation into account may lead to spurious signatures of effective population size changes through time [START_REF] Chikhi | The confounding effects of population structure, genetic diversity and the sampling scheme on the detection and quantification of population size changes[END_REF][START_REF] Peter | Distinguishing between population bottleneck and population subdivision by a Bayesian model choice procedure[END_REF][START_REF] Stadler | The impact of sampling schemes on the site frequency spectrum in nonequilibrium subdivided populations[END_REF] and to biased estimates of the mutation rate [START_REF] Navascues | Elevated substitution rate estimates from ancient DNA: model violation and bias of Bayesian methods[END_REF]. Moreover, models of population structure can also leave a genetic footprint similar to episodes of ancient hybridization [START_REF] Eriksson | The doubly conditioned frequency spectrum does not distinguish between ancient population structure and hybridization[END_REF].

Of particular importance in this context are the Range Expansions (RE). Episodes of RE have been postulated for many species [START_REF] Boessenkool | Relict or colonizer? Extinction and range expansion of penguins in southern New Zealand[END_REF][START_REF] Brown | Range expansion and hybridization in Round Island petrels (Pterodroma spp.): evidence from microsatellite genotypes[END_REF][START_REF] Cavalli-Sforza | The application of molecular genetic approaches to the study of human evolution[END_REF][START_REF] Klopfstein | The fate of mutations surfing on the wave of a range expansion[END_REF][START_REF] Rendine | Simulation and separation by principal components of multiple demic expansions in Europe[END_REF][START_REF] Ibrahim | Spatial patterns of genetic variation by different forms of dispersal during range expansion[END_REF] but they have been quantitatively tested in a spatially explicit framework only in few studies [START_REF] Barbujani | Indo-European origins: a computer-simulation test of five hypotheses[END_REF][START_REF] Francois | Demographic history of european populations of Arabidopsis thaliana[END_REF][START_REF] Gaggiotti | Disentangling the effects of evolutionary, demographic, and environmental factors influencing genetic structure of natural populations: Atlantic herring as a case study[END_REF][START_REF] Hamilton | Bayesian estimation of recent migration rates after a spatial expansion[END_REF][START_REF] Neuenschwander | Colonization history of the Swiss Rhine basin by the bullhead (Cottus gobio): inference under a Bayesian spatially explicit framework[END_REF], Ray et al., 2005[START_REF] Schneider | Signals of recent spatial expansions in the grey mouse lemur (Microcebus murinus)[END_REF][START_REF] Mona | Investigating sex-specific dynamics using uniparental markers: West New Guinea as a case study[END_REF]. REs leave distinct signatures in the genetic diversity of populations [START_REF] Barbujani | Indo-European origins: a computer-simulation test of five hypotheses[END_REF][START_REF] Edmonds | Mutations arising in the wave front of an expanding population[END_REF][START_REF] Ray | Intra-deme molecular diversity in spatially expanding populations[END_REF][START_REF] Sousa | Impact of range expansions on current human genomic diversity[END_REF][START_REF] Travis | Deleterious mutations can surf to high densities on the wave front of an expanding population[END_REF], proceeding as a series of founder events followed by continuous migrations among neighboring demes [START_REF] Austerlitz | Evolution of coalescence times, genetic diversity and structure during colonization[END_REF][START_REF] Ray | Intra-deme molecular diversity in spatially expanding populations[END_REF]. In a seminal work on the genetics of metapopulation, [START_REF] Wakeley | Nonequilibrium migration in human history[END_REF][START_REF] Wakeley | Gene genealogies in a metapopulation[END_REF] has elegantly separated the coalescent process of genes sampled in a structured population into two separate phases, namely the scattering and the collecting phase. The scattering phase, proceeding backward in time, starts at the present generation, ending when all sampled genes have either coalesced or migrated to different demes. At this point, the collecting phase begins: this period is characterized by a large number of migration events and a few coalescent events that can happen only when two lineages end up in the same deme. The coalescent process during the collecting phase in a structured population was shown to have the same properties as in an unstructured population, with time scaled by the number of demes, deme size, and migration rates [START_REF] Wakeley | Nonequilibrium migration in human history[END_REF]. A similar separation of time scales was also shown to hold in more complex and realistic models taking into account geographic features [START_REF] Wakeley | Gene genealogies in a metapopulation[END_REF], in extinction/recolonization models [START_REF] Wakeley | Metapopulation models for historical inference[END_REF][START_REF] Wakeley | Gene genealogies in a metapopulation[END_REF], in continuous populations [START_REF] Wilkins | A separation-of-timescales approach to the coalescent in a continuous population[END_REF][START_REF] Wilkins | The coalescent in a continuous, finite, linear population[END_REF] as well as in range expansions [START_REF] Mona | Genetic consequences of habitat fragmentation during a range expansion[END_REF][START_REF] Nullmeier | The coalescent in boundary-limited range expansions[END_REF][START_REF] Ray | Intra-deme molecular diversity in spatially expanding populations[END_REF][START_REF] Wegmann | Molecular diversity after a range expansion in heterogeneous environments[END_REF].

Despite the growing interest in the genetics properties of REs, many of their features are not yet completely understood. Most of the attention has focused so far on RE in homogeneous environments [START_REF] Ray | Intra-deme molecular diversity in spatially expanding populations[END_REF], where the driving parameter of the genetic diversity is Nm (effective population size times the migration rate). This result mirrors Wright's finite island model (FIM), where Nm determines not only the diversity of the meta-population [START_REF] Wakeley | Segregating sites in Wright's island model[END_REF], but also its genetic structure [START_REF] Wright | Evolution in Mendelian Populations[END_REF]. However, REs are spatial explicit nonequilibrium models where the colonization process starts from an ancestral deme. So far, the role played by the effective population size of the ancestral deme (hereafter N anc ) on shaping the genetic variation of the meta-population has not been thoroughly investigated. Moreover, both [START_REF] Ray | Intra-deme molecular diversity in spatially expanding populations[END_REF] and [START_REF] Wegmann | Molecular diversity after a range expansion in heterogeneous environments[END_REF] noted that, for low Nm, demes with large carrying capacity have a higher number of haplotypes compared to demes with low carrying capacity.

They explained this finding suggesting that gene lineages may co-migrate in the same deme and coalesce quickly with low N. This increases the probability of identity and decreases the number of haplotypes compared to demes where N is larger. Nevertheless, no further analyses were conducted to test this hypothesis. To better characterize the effect of different N values given a fixed Nm on the pattern of genetic diversity as well as to define the influence of N anc , we simulated REs using SPLATCHE2 [START_REF] Ray | SPLATCHE2: a spatially explicit simulation framework for complex demography, genetic admixture and recombination[END_REF]. We found that for a fixed Nm: i) intrademe genetic diversity depends on N (i.e., for increasing N values we found a monotonic increase in the number of haplotypes and in the number of segregating sites); ii) increasing N leads to a skew in the site frequency spectrum towards an excess of rare variants, which leads to more significantly negative values of Tajima's D. We compared these results to an equilibrium stepping stone model: fixing Nm we observed an increase in the genetic diversity for higher N as in the REs but we found no influence on the frequency of rare variants. We also identified the role played by the ancestral population size of the founding deme particularly on the degree of population structure of the whole metapopulation. We further demonstrated how it is indeed possible to correctly jointly estimate N and m in a RE using an approximate Bayesian computation framework (ABC), provided we have sampled many demes (at least more than one) and we have typed at least few independent loci.

Materials and Methods

We simulated a range expansion of haploid individuals over a two-dimensional array of demes arranged on a lattice using the software SPLATCHE2 [START_REF] Ray | SPLATCHE2: a spatially explicit simulation framework for complex demography, genetic admixture and recombination[END_REF]. Each deme exchanges migrants with its four nearest neighbors as in a 2D stepping stone model [START_REF] Kimura | The Stepping Stone Model of Population Structure and the Decrease of Genetic Correlation with Distance[END_REF] with bouncing edges. Simulations are done in two consecutive steps, namely the forward (demographic) and the backward (coalescent) steps. The forward simulation starts from a founding (or ancestral) deme, which sends migrants to its neighboring demes. The founding deme is characterized by an effective population size N anc , a parameter which does not play a role in the forward process but it will determine the rate of coalescence in the backward phase (see below).

Migrations to empty demes represent new colonization events. Each deme has an intrinsic growth rate r (the same value for all demes) and its density is logistically regulated by its carrying capacity [START_REF] Ray | Intra-deme molecular diversity in spatially expanding populations[END_REF]. We note that for this reason the effective population size in each deme varies slightly each generation, but the effect of such oscillations on the genetic variability are negligible. For this reason and for the sake of comparison to previous work and to equilibrium models we called hereafter the carrying capacity N. After the regulation step within each deme, migrants are sent to the four neighboring demes at rate m. The process is repeated generation by generation for each non-empty deme, resulting in a wave of advance of the whole population.

The second phase of the algorithm then starts at the present generation, proceeding backward in time. The actual number of individuals present in a deme is used to compute the probability of a coalescent event, and the migration rates determine the probability of each sampled genes to emigrate, backward in time, to the surrounding demes. The coalescent process stops after all genes have coalesced. The lineages which did not coalesce in the scattering or in the collecting phase will coalesce in the ancestral deme, with a rate determined by N anc .

The size of the array of stepping stones was fixed to 50×50 demes and the origin of the expansion was set to position <1;1> of the square lattice (North-West corner). We tested six carrying capacities (N=20; 50; 100; 200; 500; 1000), with migration rates varying accordingly to get a fixed Nm. To cover observed population differentiation in real data, we tested four Nm values (Nm=1; 10; 50; 100). This value would range from ~0.005 to ~0.33 in a finite island model.

Expansion times (Texp) were set to 50,000; 15,000; 3,000; 1,000 and 500 generations ago and the ancestral population size to Nanc =100; 1,000; 10,000. These parameters combinations were chosen to represent a post-glacial expansion of species with different generation time and to explore a broad range of plausible biological demographic scenarios. We simulated DNA sequences of 1,000 bp with a mutation rate of 0.001 per generation for the whole locus and performed 1,000 coalescent simulations for each demographic scenario. We sampled 9 fixed demes in the array, as displayed in Figure 1. From each deme we sampled 20 lineages. All these simulations were repeated under an instantaneous expansion model as in [START_REF] Hamilton | Bayesian estimation of recent migration rates after a spatial expansion[END_REF], where all the demes are colonized at the same time (hereafter, IE model). Finally, we simulated an equilibrium stepping stone models for the four Nm and the same combinations of N and m as above using fastsimcoal2 v2.5.1 [START_REF] Excoffier | Robust demographic inference from genomic and SNP data[END_REF]. Simulations of the equilibrium stepping stone were performed with a mutation rate of one order of magnitude lower than the one used in REs, in order to obtain comparable levels of variability. We computed several statistics to characterize the genetic diversity at the intra-deme level: the number of haplotypes H (i.e., the number of different sequences present in the sample), the number of segregating sites S (i.e., the number of polymorphic sites in the data), the mean pairwise difference π (i.e., the mean number of differences between all pairs of sequences in the sample). To detect departure from the coalescent equilibrium model (i.e., an isolated population of constant effective size) due to the range expansion as well as to the presence of metapopulation structure we computed two summary statistics and then tested their departure from mutation-drift equilibrium: Tajima's D, which is based on the comparison of two estimator of the scaled mutation rate θ (namely θ S and θ π ); Fu' Fs which evaluates the probability of having a number of haplotypes similar or smaller than the observed value given the observed number of the mean pairwise difference. To determine the extent of genetic differentiation between demes we computed the ɸst (based on molecular distance metrics [START_REF] Excoffier | Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data[END_REF]) and Fst (based on allele frequency). All summary statistics were calculated using Arlequin 3.5 [START_REF] Excoffier | Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows[END_REF]. Intrapopulation statistics were computed on the central deme of the array as well as on the nine pooled demes. ɸst and Fst were computed among the central and the bottom left deme (Figure 1), as well as among all demes.

To test the ability to estimate the four parameters of our RE model: N, m, Nanc and Texp, as well as the Nm product, we developed an ABC framework [START_REF] Beaumont | Approximate Bayesian computation in population genetics[END_REF][START_REF] Bertorelle | ABC as a flexible framework to estimate demography over space and time: some cons, many pros[END_REF][START_REF] Beaumont | Approximate Bayesian Computation in Evolution and Ecology[END_REF] in four steps.

Step 1: we created our reference tables for the later estimation process by performing 100,000 simulations sampling from uniform prior distributions for N (U:5-2,000), m (U:0.004-1), Texp (U:500-20,000) and Nanc (U:50-15,000). The resulting compound prior for Nm is very close to an exponential distribution with rate 0.002 bounded within 0-2,000.

We sampled 9 demes and 20 lineages per deme (Figure 1) and simulated 1,000 independent loci, each of 1,000 bp length with a mutation rate of 0.001 per generation for the whole locus as above.

Step 2: we computed two sets of summary statistics for the ABC analysis: i) SET1: the intrapopulation statistics (H, S, π, D and Fs) in the 9 demes pooled and the ɸst (for a total of six summary statistics); ii) SET2: the same five intra-population statistic computed in the central deme only. Each statistics was averaged over the number of loci.

Step 3: we generated pseudo observed datasets (pods) by randomly extracting 1,000 observations from the reference table, recording the corresponding vector of prior values.

Step 4: we analyzed pods using the ABC algorithm, with a weighted local linear regression [START_REF] Beaumont | Approximate Bayesian computation in population genetics[END_REF] performed on the 5% of the closest simulations. For each analyzed pods we retained the mode and the median estimate of the posterior distribution and plotted them against the "real" value recorded in step 3. We also computed the scaled root mean square error (SRMSE) [START_REF] Walther | The concepts of bias, precision and accuracy, and their use in testing the performance of species richness estimators, with a literature review of estimator performance[END_REF] both for the median and the mode as well as the 95% coverage. A power analysis was carried out by subsampling from the total number of loci: 1, 10 and 100 loci for both sets of summary statistics and sampling scheme. To determine which statistics is more informative and the importance of pooling demes we define three other sets of summary statistics: SET3 (summary statistics computed on all demes separately plus the ɸst, for a total of 46 statistics), SET4 (the same as SET1 without ɸst), SET5 (the same as SET1 without the neutrality tests).

Results

Unless specified, we present and discuss results in REs relative to: i) the nine pooled demes concerning the statistics computed at the intra-population level (Figure 1); ii) the central deme concerning the neutrality tests; iii) the central and the bottom left deme, concerning Fst and ɸst (Figure 1). The others sampling schemes displayed a similar pattern. We note that the grid of demes was not fully colonized for Nm=1 and Texp≤1,000, therefore we did not report results of these simulations.

For fixed Nm we observe a monotonic increase in genetic diversity (as measured by H, S and π) for growing N (Figure 2). Both H and S showed a dependency on N for all the Texp tested as well as for the equilibrium stepping stone model (Figure 2, Figure S1 andS2). Conversely, we observed a different pattern for π which is insensitive to N for low Texp but depends on it for higher expansion time (Texp =50,000) and in the equilibrium stepping stone model (Figure S3).

The excess of private haplotypes for larger N leads to an increase of H and S which is reflected in the more negative D and Fs (Figure 3, Table S1, S2, S3, S4). The influence of N on the neutrality tests decreases for high Texp and it is lost when the process arrives at equilibrium (Figure 3 and4). This result can be explained as follow. For fixed Nm more mutations will occur for higher N.

The length of the internal branches of the genealogy is constrained by Texp which implies that this excess of mutations will be concentrated in the terminal branches. These new mutations are therefore at low frequencies having almost no influence on π. Conversely, for large Texp and in the stepping stone the internal branches will be longer so that mutations will have time to occur.

In this case, these mutations will be at intermediate frequencies, having an impact on π (and, as a consequence, on the neutrality tests).The percentage of significant outcomes (as calculated by Arlequin 3.5), both positive and negative, is also related to N for all the four Nm tested (Figure 4, Table S1, S2, S3, S4). Low Nm can produce false signature of bottleneck in equilibrium model at the deme level [START_REF] Chikhi | The confounding effects of population structure, genetic diversity and the sampling scheme on the detection and quantification of population size changes[END_REF]. We found a similar result in REs model, where for Nm=1 many datasets display a significant signature of reduction of the effective population size (Table S3). However, when N grows from 20 to 1,000, the number of significant positive Tajima's D drop considerably (Table S3). Similarly, the power to detect a RE depends on the Texp [START_REF] Ray | Intra-deme molecular diversity in spatially expanding populations[END_REF]): here we show that even for large Texp the probability of detecting an expansion increases substantially for larger N (Figure 4), and the probability of erroneously inferring a bottleneck drops (Table S1, S2, S3, S4).

We present the ɸst and Fst distribution in REs for Nm=10 and Nanc=100 in Figure 5 and compare the results with the equilibrium stepping stone model. We observed a negative correlation between both measures of population differentiation and N when Texp≥ 3,000 (Figure 5). This result is consistent with the monotonic increase of diversity for larger N: higher heterozigosity tends to lower the degree of population differentiation, in particular when it is measured from allele frequencies [START_REF] Hedrick | A standardized genetic differentiation measure[END_REF]. However, ɸst does not depend on N in the equilibrium stepping stone, while the behavior of Fst is consistent to that of non-equilibrium models. We found a similar trend for all the Nm values we tested, ranging from 1 to 100 (Figure S4, S5, S6).

However, the decay of population differentiation for large N it is particularly strong for Nm=1 where the median value of Fst is ~0.8 for N=20, falling to ~0.35 for N=1,000. Such difference is weaker for ɸst (decreasing from ~0.8 to ~0.6), consistent with the suggestion that measures based on molecular difference can better correct for the heterozigosity excess [START_REF] Hedrick | A standardized genetic differentiation measure[END_REF]. For Texp≤1,000 the pattern is more complex and it is not consistent between the two differentiation measures. First of all, we observed a slight increase in the mean and in the variance of the distribution of ɸst for higher N (Figure 5, 6 and 7, Figure S5 andS6), while no dependency on N was found in the Fst. Secondly, the monotonic increase in ɸst for larger N is more pronounced when Nanc=10,000. This pattern is more pronounced in IEs than in REs, where we found similar results also for Texp=3,000 when Nanc=10,000 (Figure 7, Figure S6). Independently of Nm, no differences related to N were observed for the Fst. The genome wide distribution of measures of population differentiation is therefore dependent on the interaction between N, Texp and Nanc in RE model. Detecting outliers loci in genomic scan of selection will require a good knowledge of all these parameters when investigating structured species.

We performed a cross-validation ABC procedure to investigate the power to jointly estimate N and m in RE model. We presented the modes of the estimated N in Figure 8 and Figure S7 for all the combination of sets of summary statistics (SET1 and SET2) and number of loci (1,10,100 and 1,000). The black line represents a perfect match: the closer the points are to the line, the closer is the estimated mode to the "real" parameter. As expected, increasing the number of loci reduces the variance of the modes around the line. When we compared SET1 and SET2, we found that pooling more demes increases the performance of the ABC procedure, particularly when a small number of loci are used. A consistent pattern was found for all other parameters (Figure S8, S9, S10), with the exception of N anc , poorly estimated in all cases (Figure S11). We reported the Empirical population genetics is now experiencing a revolution. Thanks to the spread of next generation sequencing technology, an impressive amount of genetic data is now easy to collect in both model and non-model organisms. This huge amount of data carries precious information on the demographic and selective history of populations and species. However, having a lot of data has a potential drawback: the interpretation of genetic variation can be significantly wrong when a wrong demographic model is applied. Most (if not all) species are structured, i.e. they are composed by many demes exchanging migrants to various extent. The two structured models which have been explored the most in theoretical population genetics are the island and the stepping stone model, both of which are at equilibrium. It is therefore urgent to explore in greater details non-equilibrium structured model and to integrate spatial components. It is for this reason that here we explored RE. Previous studies have shown that in homogeneous environment the main driver of genetic diversity after a RE is the compound Nm product [START_REF] Ray | Intra-deme molecular diversity in spatially expanding populations[END_REF]. Here we focus on the influence of N for fixed Nm values and on the role played by Nanc in shaping the genetic diversity. Finally, we tested how well we can recover the demographic parameters after a RE using an ABC approach and how many demes it is important to sample to obtain accurate estimates. All those issues have not been investigated before and they are likely to be informative for empirical studies.

To better understand the specificity of REs we compared our results with an equilibrium stepping stone model, fixing the Nm product and exploring various combinations of N and m. We considered both the variability within deme and the population differentiation computed among demes. At the intra-deme level, we found a monotonic increase in genetic diversity in both RE and stepping stone for increasing N. However, we found two important characteristics of RE: i) an excess of of low frequency haplotypes for increasing N; ii) a complex interaction between N, N anc and Texp in shaping the degree of population differentiation. It is for these reasons that we could estimate N and m using an ABC approach. Below we provide some coalescent arguments to explain our findings.

Signatures of RE

Coalescent in meta-population can be described in terms of scattering and collecting phases [START_REF] Wakeley | Nonequilibrium migration in human history[END_REF]. In RE models the number of emigrants sent out each generation is Nm while the within deme probability of coalescent is determined by N only. Therefore, for a fixed Nm the scattering phase will be longer for higher N and fewer within deme coalescences will happen. This has three consequences, as for higher N: i) the number of lineages entering the collecting phase will be higher; ii) there will be more H and more private S since in RE models mutations can also occur during the scattering phase; ii) the ratio of external to internal branch lengths increases. Indeed, we observed a monotonic increase both in the number of haplotypes and in the number of segregating sites for larger N when we fixed Nm (Figure 2, Figure S1). This is true also for the equilibrium stepping stone and it is therefore related to the increase of the scaled mutation rate (Figure S2). However, the distribution of mean pairwise differences follows a different pattern: we observed a monotonic increase for larger N only when Texp ≥ 50,000 and at equilibrium. Both θ π and θ s are estimators of the scaled mutation rate θ and Tajima's D is based on their difference, which is centered around zero under neutrality [START_REF] Tajima | Statistical method for testing the neutral mutation hypothesis by DNA polymorphism[END_REF]. For a fixed Nm if S increases monotonically with N but π remains constant the difference between the two estimators will become negative, a signature typical of population expansion. For this reason, the power to detect an RE for fixed Nm depends heavily on N (Figure 4). This has practical implications: the signature of a RE will be lost more quickly in species characterized by smallsized demes. This is clearer when we compare the results of Nm=10 between RE and the equilibrium stepping stone, where N does not a play a role in determining the shape of the genealogy (Figure 3 and4). Moreover, in equilibrium model low Nm give rise to false signature of population bottleneck [START_REF] Chikhi | The confounding effects of population structure, genetic diversity and the sampling scheme on the detection and quantification of population size changes[END_REF]. This is true also in RE [START_REF] Ray | Intra-deme molecular diversity in spatially expanding populations[END_REF] but the number of positive outcomes of Tajima's D drops when the deme size increases (Table S3). We argue that in empirical studies, particularly when concerning endangered species, N and m must be estimated separately to fully recover the underlying demography.

Population differentiation

The collecting phase dominates the coalescent process in equilibrium model [START_REF] Wakeley | Nonequilibrium migration in human history[END_REF], while its length in RE/IE will depend both on N and on Texp. In RE/IE, backward in time lineages will contract towards the founding/ancestral deme: for a fixed Nm there will be less coalescence with higher N and more lineages entering the ancestral phase of the coalescent process. The number of lineages surviving until the ancestral phase is also influenced by Texp. Indeed, when Texp is low, there is less opportunity for the sampled lineages to coalesce in the array of demes.

The surviving lineages will coalesce in the ancestral deme according to a Kingman coalescent dependent on N anc only, accumulating mutations proportionally to N anc and the mutation rate.

Therefore, for fixed Nm larger N will increase the heterozygosity and the number of lineages reaching the founding deme, the two phenomena having opposite effect on the degree of population differentiation. The increased within deme heterozigosity tends to lower both ɸst and Fst [START_REF] Hedrick | A standardized genetic differentiation measure[END_REF][START_REF] Meirmans | Using the AMOVA framework to estimate a standardized genetic differentiation measure[END_REF], causing a negative correlation between N and the genetic differentiation between demes (Figure 5, S4). This is clearer for Texp≥3,000 because the difference in heterozigosity among the tested N is of higher magnitude. We found a similar behavior in the equilibrium stepping stone (Figure 5), consistently with the observed pattern of intra-deme diversity. Conversely, more lineages coming from different demes will arrive in the ancestral phase for larger N and Texp<3,000. Here, the occurring mutations enhance the between deme component of molecular variance, increasing ɸst but not Fst for increasing N (Figure 5, 6, 7, S5, S6). When Texp≤1,000 generations, under our parameters combination, the influence of N on the ɸst is remarkable (Figure 4). The larger the Nanc, the higher will be the average ɸst and its variance (Figure 6, 7, S5, S6). This effect is of greater intensity for IE (Figure 7, S6), where the number of lineages arriving in the ancestral phase is larger because the contraction phase typical of RE, where many coalescent happen before arriving in the ancestral deme [START_REF] Excoffier | Genetic consequences of range expansions[END_REF], is absent. These findings have important implications in the analysis of real data. For example, invasive species are typical instances of young REs: the genome wide distribution of the ɸst will differ enormously between species characterized by small vs large sized demes.

Many tests of selection are based on the idea of finding outliers in the distribution of some summary statistics such as the ɸst. A good knowledge of the effective deme size is therefore needed to generate a correct null distribution to avoid false positive results in genomic scan of selection.

Estimating parameters of a RE

The influence of the carrying capacity on the genetic variability both at the intra and interpopulation level suggested that it should be possible to jointly estimate N and m using summary statistics. To this end, we developed an ABC framework and check by cross-validation the power to estimate the parameters of a RE. We varied the number of simulated loci, the number of demes sampled and the set of summary statistics chosen. To be consistent with all the previous simulations, we computed inter and intra-population summary statistics both in a single deme and in a pool of lineages coming from all 9 demes. We measured the accuracy of the estimation of all paremeters by computing the SRMSE averaged through the prior distributions. For any ɸst) may provide more inferential power on the estimation of this parameter, which will be the subject of further investigations.

Conclusions

We show that for fixed Nm the carrying capacity of demes play an important role during a RE, both at the intra and at the inter-population level. Species are unlikely to be at equilibrium. We therefore stress the importance of estimating N and m separately as well as investigating the ancestral phase of the coalescent process. Our results highlight that we can actually disentangle N and m provided many demes and a consistent set of summary statistics is used in the estimation process (for example using an ABC algorithm). This would be extremely useful in conservation studies and also when investigating invasive species. Finally, we note that species undergoing recent REs will retain the signature of the ancestral demography, which could be uncovered focusing on the genome-wide distribution of population differentiation.

Titles and legends to figures

Figure 1 .

 1 Figure 1. The simulated array of 50x50 demes. Black dots: sampled demes. Gray dot: founding deme (origin of the range expansion). Grey circles: demes sampled to compute ɸst and Fst presented in the text.

Figure 2 .

 2 Figure 2. Distribution of the number of haplotypes (H) for Texp =3,000 and Nanc =100 for four Nm values computed on the pooled demes. N values on the x-axis.

Figure 3 .

 3 Figure 3. Distribution of Tajima's D for Nm=10 in a RE for various values of Texp and the equilibrium stepping stone model in the nine pooled demes

Figure 4 .Figure 5 .

 45 Figure 4. Proportion of significant negative outcomes of Tajima's D for the four Nm values as a function of N in the nine pooled demes

Figure 6 .

 6 Figure 6. Distribution of ɸst (on the left of the dotted line) and Fst (on the right of the dotted line) for Nm=10 in a RE for two values of Texp.

Figure 7 .

 7 Figure 7. Distribution of ɸst (on the left of the dotted line) and Fst (on the right of the dotted line) for Nm=10 in an IE for two values of Texp.

Figure 8 .

 8 Figure 8. Cross-validation of N. Modes (y-axis) are estimated under the ABC framework using the SET1 and plotted against the real values (x-axis). The analyses are performed with 1, 10, 100 and 1,000 loci (see text for details).

Acknowledgments

We thank Laurent Excoffier for stimulating the present study and Oscar Lao for interesting comments. We thank the anonymous reviewers for greatly improving a previous version of this manuscript. We are grateful to the Genotoul bioinformatics platform Toulouse Midi-Pyrenees for providing computing resources.

DATA ARCHIVING

There were no data to deposit.

SRMSE averaged over the priors for both the mode and the median for all the four parameters of the RE model and the Nm product in Table S5. An SRMSE of one means that mode (or median) is on average one time higher or smaller than the real value. SRMSE decreases with the number of loci and it is always smaller in SET1 compared to SET2. We noted that the accuracy in the estimate of Nm is slightly higher than that of N, while m it is considerably worse than both of them. Adding more loci did not improve the estimate of Nanc, which was unexpectedly poor.

Conversely, few loci were enough to correctly estimate Texp with any sets of summary statistics considered. We further tested three more combination of summary statistics. First, we evaluated whether pooling more demes increases the performance in the ABC estimation process. To this end, we used the intra-deme summary statistics calculated separately in the nine demes: adding ɸst we defined SET3, characterized by 46 summary statistics. SRMSE were similar to SET1 (Table S5), suggesting that it is important to sample more demes but not necessarily to pool them.

We note however that pooling reduces the number of predictors in the regression step of ABC and it could be particularly helpful when studying large datasets (i.e., composed by many demes).

We then investigated the role played by ɸst by removing it from SET1 (thus obtaining the SET4).

We obtained results very similar to SET2 (Table S5): migration rate is poorly estimated and, as expected, Nm shows an SRMSE twice as big as the one found in SET1. Finally, we obtained SET5 by removing D and Fs from SET1. We observed a decrease in the performance of the estimation of Texp, N and m, all of which showed an SRMSE two times bigger than SET1.

Conversely, the compound parameter Nm was unaffected.

Discussion

combination of loci/demes/summary statistics we obtained a similar accuracy in the estimation of N and Nm (Table S5). This shows unambiguously how N affects the distribution of genetic variability, playing a role that can be easily disentangled from m. Several independent loci need to be typed to take into account the stochastic variability of the coalescent process: however, it seems that going from 100 to 1,000 loci does only slightly improve the accuracy of our estimate, suggesting that at least under our model a relatively small set of loci should provide consistent results. When using summary statistics computed in pooled lineages (SET1) or in the nine demes independently (SET3) the performance of the ABC procedure increases drastically compared to SET2, where statistics from a single deme are used (Figure 8, S7, S8, S9, S10, S11). The number of demes studied is therefore important to correctly characterize the parameter of REs. Summary statistics convey different information on each parameter. As expected, to estimate the migration rate m is crucial to add the ɸst in the ABC algorithm (Table S5). Indeed even when we pooled demes the SRMSE is higher than one for m if we do not include the ɸst while the drop in accuracy is less drastic for the compound parameter Nm (SET1 and SET4, Table S5). Similarly, neutrality test statistics are also useful to disentangle N and m (Table S5, SET5), which is consistent with the finding that N influences the power to detect a RE for fixed Nm. It should be possible to devise an optimal sampling strategy in terms of number of demes and number of loci, but this is beyond the scope of this paper. We only highlight here the importance of sampling more than just one deme to estimate the parameters of a RE and to compute a measure of population differentiation to estimate the migration rate. Despite the large number of loci studied, we could not correctly estimate the N anc . Contrary to N, which plays a direct role on the intrapopulation diversity and on the measure of differentiations, N anc interact with N and Texp to shape Fst and ɸst. We argue that for this reason it may be not identifiable under our setting. At the same time, we note that other summary statistics (for example combing the variance of both Fst and

CONFLICT OF INTEREST

The authors declare no conflict of interest. Supplementary Information accompanies this paper on Heredity website (...)

Literature cited