Julien Clement
email: jclement@liafa.jussieu.fr

Carole Delporte-Gallet

Hugues Fauconnier

Mihaela Sighireanu
email: sighirea@liafa.jussieu.fr

Guidelines for Verification of Population Protocols

Keywords: population protocols, Petri Nets, vector addition systems, model-checking, LTL, Spin, Prism

We address the problem of verification by modelchecking of basic population protocol (PP) model of Angluin et al. This problem has received special attention the last two years and new tools have been proposed to deal with. We show that the problem can be solved using the existing modelchecking tools, e.g., Spin and Prism. For this, we apply the counting abstraction to obtain an abstract model of a PP which can be efficiently verified by the existing model-checking tools. Moreover, this abstraction is preserves the correct stabilization property of PP. To deal with the fairness assumed in the PP model, we provide two recipes. The first one gives sufficient conditions under which the PP fairness can be replaced by the weak fairness implemented in Spin. We show that this recipe can be applied to several PP models. In the second recipe, we show how to use the probabilistic model-checking and the tool Prism to deal completely with the PP fairness. The correctness of this recipe is based on existing theorems on finite discrete Markov chains.

I. INTRODUCTION

Population protocols [START_REF] Angluin | Computation in networks of passively mobile finite-state sensors[END_REF] are a formal model for sensor networks. Since such networks should become a part of the every day life, it is important to ensure their correctness. This task is now possible in a fully automatic way due to the important progress of verification techniques. Proposed in the early eighties, the model-checking technique [START_REF] Emerson | Characterizing correctness properties of parallel programs using fixpoints[END_REF], [START_REF] Queille | Iterative methods for the analysis of petri nets[END_REF] allows to check complex correctness properties for models with finite number of states. The system to be verified is formalized in some high level formal language, e.g., CSP or Petri Nets, to obtain a model. The requirements of the system are specified using some logical language, e.g. LTL [START_REF] Pnueli | The temporal logic of programs[END_REF], to obtain a property. Model-checking techniques explore exhaustively the model in order to check that each state of this model satisfies the given property (a reference for model-checking techniques is [START_REF] Clarke | Model Checking[END_REF]). Very performant modelcheckers, e.g., Spin [START_REF] Holzmann | Spin model checker, the: primer and reference manual[END_REF], are now available. They are able to deal with systems of more than 10 20 states. Moreover, the model-checking techniques have been extended to deal with more complex finite models like the probabilistic ones, e.g., in the Prism [START_REF] Kwiatkowska | Probabilistic symbolic model checking with PRISM: A hybrid approach[END_REF] model-checker.

Unfortunately, some models of realistic systems have a number of states which exceeds the capacity of finite modelcheckers. Thus, the research have focalized on designing abstraction techniques to reduce the model-checking problem of huge systems to the model-checking of smaller systems. One of the first abstraction technique proposed is the counting abstraction [START_REF] German | Reasoning about systems with many processes[END_REF]. It reduces systems with many identical processes running in parallel to a system which keeps track only of the number of processes which are in some particular state. This abstraction has shown its efficiency for the model-checking of safety properties in cache coherence protocols.

In this work, we apply finite model-checking tools and the counting abstraction technique to verify population protocols.

The model of population protocols (PP) [START_REF] Angluin | Computation in networks of passively mobile finite-state sensors[END_REF] involves individual agents with a very simple behavior, which can be seen as a finite state machine. An important property is the uniformity of the protocol, i.e., the fact that the protocol description is independent on the number of agents (called also the population size) or their identity. When two agents come into range of each other ("meet"), they can exchange information. The agents are anonymous and move in an asynchronous way. At the beginning of the protocol, each agent receives a piece of input. The goal of the protocol is to stabilize each agent into a state in which it outputs the value to be computed by the protocol. Thus, the specification of a population protocol includes the "correct stabilization" (CS) property which requires that each agent stabilizes its computation to an output value which corresponds to a function on the input assigned to agents.

Our contribution is to highlight the use of the existing model-checking techniques to check the correct stabilization property for protocols whose sizes is a fixed constant. For this, we use the counting abstraction to reduce the PP model to a vector addition system model, which is a model computationally equivalent to classical place/transition Petri Nets. We show in Section V that this abstraction allows the verification of the correct stabilization property in a more efficient manner. We also give sufficient conditions to check the CS property under a fairness constraint weaker than the one required in [START_REF] Angluin | Computation in networks of passively mobile finite-state sensors[END_REF]. Based on these theoretical results, we show in Section VI how to perform verification using the Spin tool. Spin is able to verify our benchmark of PP for population size greater than 10 3 . The second recipe is presented in Section VII. We highlight that the fairness of the PP model can be exactly captured for the CS property in a probabilistic model where probabilities strictly greater than some ε > 0 are assigned to transitions. Then, we show that the Prism tool is very efficient in the model-checking this property on the abstract model of PP. In addition, we provide a recipe to compute, using Prism, the average number of interactions before the convergence for a PP model. Section IX concludes our work and gives some directions of future research.

Related work: The verification of systems with a huge or unbounded number of processes using the counting abstraction has been introduced in [START_REF] German | Reasoning about systems with many processes[END_REF]. The PP enter in the class of systems considered by this work. However, the fairness assumption they consider is not the same. We use their abstraction technique to obtain a model on which the finite model-checking tools can be applied. Although we fix the number of processes, we deal with a more complex fairness condition.

In several papers, Pang et al. [START_REF] Pang | On automatic verification of self-stabilizing population protocols[END_REF], [START_REF] Bakhshi | Leader election in anonymous rings: Franklin goes probabilistic[END_REF], [START_REF] Liu | Verification of population ring protocols in pat[END_REF] consider the verification by model-checking of leader election PP on complete graph or ring networks. For complete graphs topologies, they apply directly Spin on the PP model, and thus they obtain results for small population sizes. For the leader election protocol in ring networks, they show that a special fairness is necessary to obtain the correct stabilization. Thus, Spin tool returns false negative for this protocol because it implements a weaker fairness. In, [START_REF] Bakhshi | Leader election in anonymous rings: Franklin goes probabilistic[END_REF], the authors apply Prism to verify the probabilistic version of leader election on rings. These works lead to the development of a new model-checker, called PAT [START_REF] Sun | Pat: Towards flexible verification under fairness[END_REF]. PAT does model-checking under different fairness conditions including the fairness needed by the leader election PP on rings. However, for the weak fairness, PAT is less efficient than Spin, and for the other fairness conditions, it is able to verify leader election PP with size smaller than 8. In this work, we consider only PP on fully connected graphs, and we manage to verify them for larger populations. Also, we prove theorems on counting abstraction for PP that can be used by any model-checker for the verification of PP models. [START_REF] Chatzigiannakis | Algorithmic verification of population protocols[END_REF] shows that the algorithmic verification of finite instances of PP is NP-hard. It provides an incomplete, polynomial algorithm for checking the correctness of the computation for PP whose size is under some fixed constant. We consider the same problem and we prove some theorems in order to use the off the shelf tools for model-checking. With these tools, we obtain better experimental results than those reported by [START_REF] Chatzigiannakis | Algorithmic verification of population protocols[END_REF].

II. LABELED TRANSITION SYSTEMS

At the semantic level, the population protocol model and the verification models we consider are finite labeled transition systems.

Definition 1: A finite labeled transition system (LTS) is a tuple T = S, I, Σ, ∆ where S is a finite set of states, I ⊆ S is the set of initial states, Σ is a finite set of labels, and ∆ ⊆ S × Σ × S is a finite set of labeled transitions. We also use the notations s ℓ -→ s ′ for (s, ℓ, s ′ ) ∈ ∆. The notation s -→ s ′ stands for ∃ℓ ∈ Σ such that (s, ℓ, s ′ ) ∈ ∆. We say that a label ℓ is enabled in a state s iff there exists s ′ such that s ℓ -→ s ′ . This model can be extended to obtain a probabilistic model by associating probabilistic distribution to transitions, as shown in Section VII.

An execution E of an LTS is a sequence of alternating states and labels E = s 0 , ℓ 0 , s 1 , ℓ 1 , . . . , s i , ℓ i , . . . such that s 0 ∈ I and s i ℓ i -→ s i+1 for all i ≥ 0. A computation C of an LTS is a maximal execution, i.e., either a finite execution s 0 , ℓ 0 , s 1 , ℓ 1 , . . . , s n such that for all s ∈ S and ℓ ∈ Σ, s n ℓ -→ s ∈ ∆, or an infinite execution. A label ℓ is enabled (resp. fired) in an execution E iff there exists i s.t. ℓ is enabled in s i (resp. ℓ i = ℓ).

We are interested in computations of LTS which do not avoid systematically some states or transitions of the LTS, i.e., they are fair. A fairness constraint restricts the set of system computations to ones which are fair. We give below three fairness constraints used in the remainder of the paper. Let C = s 0 , ℓ 0 , s 1 , ℓ 1 . . . , s i , ℓ i , . . . be a computation of an LTS.

Definition 2: C is strong locally fair (LF) iff for every ℓ, if ℓ is enabled in infinitely many states s i in C, then ℓ = ℓ j for infinitely many j in C, i.e., ✷✸ℓ enabled ⇒ ✷✸ℓ fired.

LF constraint corresponds to the classic strong fair [START_REF] Francez | Fairness[END_REF] constraint. Intuitively, it asks that if a label is enabled infinitely often in a computation, it is also executed infinitely often by the computation.

A stronger version of fairness is obtained by looking at the states in the neighborhood of a computation, i.e., states that can be reached using any transition from the states of a computation.

Definition 3: C is strong globally fair (GF) iff for every s and s ′ such that s -→ s ′ , if s = s i for infinitely many i in C, then s j = s ′ for infinitely many j in C i.e., s -→ s ′ ∧ ✷✸s ⇒ ✷✸s ′ .

In a GF computation C, any neighbor s ′ of a state s appearing infinitely often in C appears also infinitely often in C. Note that the neighbor state s ′ may be a direct successor of s through transitions in ∆ labeled with different labels. Also, GF does not require to have s ′ a direct successor of s in C. The GF constraint corresponds to the fairness constraint required for the PP in [START_REF] Angluin | Computation in networks of passively mobile finite-state sensors[END_REF].

Definition 4: C is weakly fair (WF) iff for every ℓ, if it exists j and s ′ i with s i ℓ -→ s ′ i ∈ ∆ for any i ≥ j, then s k ℓ -→ s k+1 for infinitely many k > j in C, i.e., ✸✷ℓ is enabled ⇒ ✷✸ℓ is fired.

WF constraint states that if a label becomes enabled forever after some steps, then it must be fired infinitely often. The weak fairness have been well studied and the verification under weak fairness is supported by performant tools like Spin [START_REF] Holzmann | Spin model checker, the: primer and reference manual[END_REF].

Let F be one of the fairness above (F ∈ {LF, GF, WF}). We use the notation comp F (T ) for the set of F computations of the LTS T .

From the algorithmic point of view, testing that a computation is strong fair (i.e., LF or GF) requires to compute its strongly connected components. Testing that a computation is weak fair can be done in constant time [START_REF] Holzmann | Spin model checker, the: primer and reference manual[END_REF], i.e., very efficiently.

III. POPULATION PROTOCOL MODEL

In this section, we briefly introduce the population protocol model. More details are available in [START_REF] Angluin | Computation in networks of passively mobile finite-state sensors[END_REF].

For a population of k agents, Π = {π 1 , . . . , π k } denotes the set of agents. Each agent in the system is modeled as a finite state machine, representing the program of the agent. These programs are uniform: each agent executes the same finite state machine which form does not depend on the number of agents in the system. This assumption makes the model strongly anonymous: the agents can not store a unique identifier. When two agents meet each other, they may interact and perform a change in their states. The underlying network of communication is a complete graph, every pair of agents may meet. The initial state of each agent is set by the protocol. The goal of a PP is to compute by each agent the same value called the output of the protocol. The value to be computed is a predicate on the initial state of the protocol.

The following definition [START_REF] Angluin | Computation in networks of passively mobile finite-state sensors[END_REF] indicates how to specify population protocols. In this definition, the initial state of an agent is fixed (using ι) by an input value in X received by the agent. In each state q, an agent outputs a value in Y given by the mapping ω. The interaction between agents is described by the relation δ: if two agents in states q 1 and q 2 meet each other, they change into states q ′ 1 resp. q ′ 2 , where (q 1 , q 2 , q ′ 1 , q ′ 2 ) ∈ δ. We also use notation (q 1 , q 2 ) -→(q ′ 1 , q ′ 2 ) for the elements of δ. In this paper, we restrict ourselves to PP whose specifications satisfy the following conditions: H 1 : the protocol computes predicates, i.e., Y = {0, 1}; H 2 : the transition relation is non-deterministic but symmetric, i.e., δ ⊆ P ≤2 (Q) × P ≤2 (Q) where P ≤2 (Q) is the set of non empty subsets of Q with cardinality at most 2. Hypotheses H 1 -H 2 simplify our proofs, but do not restrict the power of the PP model. To emphasize the symmetry of transitions required by H 2 , we use a special operator instead of tuples to specify a transition in δ, i.e., q 1 q 2 -→ q 3 q 4 (q 1 and q 2 or q 3 and q 4 may be equal). Note that a transition exists for any couple of states q 1 and q 2 ; when is no explicit transition with left side q 1 q 2 in δ, the default transition q 1 q 2 -→ q 1 q 2 is applied. We denote by U(δ) the elements of δ which are not identities, i.e., {q 1 , q 2 } = {q 3 , q 4 }.

The semantics of a PP P is given by an LTS T P = S, I, Σ, ∆ defined by:

• a state s (called configuration in [START_REF] Angluin | Computation in networks of passively mobile finite-state sensors[END_REF]) is a mapping Π → Q specifying the state of each agent; then S = 2 Π→Q , • the initial set of states I, is given by the image of the input of the protocol by the function ι, • the set of labels Σ ⊆ P 2 (Π) × δ, and

• the set of labeled transitions ∆ = {s ℓ -→ s ′ | ∃π i , π j ∈ Π s.t. t = s(π i ) s(π j ) - → s ′ (π i ) s ′ (π j ) ∈ δ and ℓ = ({i, j},t)}.
Intuitively, an interaction between two agents in the PP model is identified at the semantic level by the set of identities of the agents interacting and by the transition in δ which is applied by this interaction. The presence of an interaction in some state s of the protocol is denoted by a transition s -→ s ′ at the semantic level. An important semantic notion for the PP is the stable output of a computation. To define it, we use the natural ex-

tension of ω to Q k → Y k . A computation C = s 0 , ℓ 0 , s 1 , ℓ 1 , . . . stably outputs some b ∈ Y if there exists some i such that ω(s j ) = b k
for all j ≥ i, i.e., every agent eventually outputs b and never changes its output thereafter, i.e., C satisfies ✸✷ ω(s) = b k .

One of the characteristics of the PP model is that the order in which pairs of agents interact is unpredictable. In order to model this aspect, one may imagine the presence of a scheduler which is scheduling the interactions. This scheduler may force two agents to never interact. In the presence of such scheduler, the PP has no chance to compute its goal. Thus, Angluin et al. [START_REF] Angluin | Computation in networks of passively mobile finite-state sensors[END_REF] require that the scheduler allows only GF computations.

Under the GF requirement for the scheduler, a protocol of size k stably computes a function f : X k → Y if, for every input x ∈ X k , every GF computation starting in ι( x) stably outputs f ( x), i.e., ι( x) ∧ GF ⇒ ✸✷ ω(s) = f ( x) . The important result of Angluin et al. [START_REF] Angluin | The computational power of population protocols[END_REF] is that a predicate is stably computable by the PP model iff it can be defined as a first-order logical formula in Presburger arithmetic.

We shortly present two examples of PP that will be our running examples. For more examples, see e.g., [START_REF] Aspnes | An introduction to population protocols[END_REF].

Example 1 (Majority): Let X = {L, F}. The MAJORITY protocol stably computes the predicate |L| > |F| where |L| (resp. |F|) is the number of agents with input L (resp. F). The states of the protocol are Q = {L, F, 0, 1}, ι is the identity function, and ω maps L and 1 to 1, F and 0 to 0. The transition set δ is defined as follows: [START_REF] Aspnes | An introduction to population protocols[END_REF] notices that the execution in Figure 1 is not GF because the state (0, L, 1) appears infinitely often but not the state

t 1 : L F - → 0 0 t 2 : L 0 - → L 1 t 3 : F 1 - → F 0 t 4 : 0 1 - → 0 0
(L, L, F) (0, L, 0) (1, L, 0) (0, L, 0) (0, L, 1) ({1, 3},t 1 ) ({1, 2},t 2 ) ({1, 3},t 4 ) ({2, 3},t 2 ) ({1, 3},t 4 ) Figure 1
. LF execution which is not GF in the MAJORITY protocol.

(1, L, 1) which is its neighbor (by transition ({1, 2},t 2 )).

However, this execution is LF because each infinitely often enabled label is fired infinitely often. From the Figure 1, it appears that this execution is also accepted by a fairness constraint which requires that any two agents meet infinitely often.

Example 2 (Threshold): Let N ≥ 2 be an integer value and X = {0, 1}. The THRESHOLD(N) protocol stably computes the function |ι(1)| ≥ N. The states of the protocol are Q = {0, . . . , N}. ι maps all agents with input 1 in state 1 and all other agents in state 0. The output function ω maps a state q to the boolean value q = N. δ is defined as follows:

t q,q : q q - → q q + 1 if 1 ≤ q < N t N,q : N q - → N N if 0 ≤ q < N
Notice that the size of Q (resp. δ) for this protocol is linear (resp. quadratic) on N.

IV. ABSTRACT MODEL

We introduce an abstract model for the PP that facilitates the verification task. This model has been studied in [START_REF] Karp | Parallel program schemata[END_REF].

Intuitively, we abstract a state s of the T P by a vector of integers c indexed by Q such that c[q] is the number of agents in state q in s. This abstraction is possible due to the uniformity of the protocol: the behavior of an agent depends only on its state and not on its index.

Another intuition of this abstraction may be obtained by the translation of the PP model into a classic place/transition Petri net (PN) model. To each state in Q is associated a place in the PN, and to each transition in δ is associated a transition in the PN. Agents are represented by tokens in the PN; the firing of a transition in the PN corresponds to an interaction between agents in the PP model. Figure 2 pictures the PN obtained for the MAJORITY protocol with k = 3. The counting abstraction keeps track of the number of tokens in each place of the PN by associating a counter with each place (state in Q); the transitions of the PN are encoded into conditions and assignments on these counters.

Formally, the model obtained by this abstraction is a vector addition system model; it is equivalent (from the point of view of computation power) with the PN model, but more convenient to encode in the verification tools we consider.

i guard(φ i ) trans(φ i ) φ 1 : c L ≥ 1 ∧ c F ≥ 1 ∧ c ′ L = c L -1 ∧ c ′ F = c F -1 ∧ c ′ 0 = c 0 + 2 φ 2 : c L ≥ 1 ∧ c 0 ≥ 1 ∧ c ′ 0 = c 0 -1 ∧ c ′ 1 = c 1 + 1 φ 3 : c F ≥ 1 ∧ c 1 ≥ 1 ∧ c ′ 1 = c 1 -1 ∧ c ′ 0 = c 0 + 1 φ 4 : c 0 ≥ 1 ∧ c 1 ≥ 1 ∧ c ′ 1 = c 1 -1 ∧ c ′ 0 = c 0 + 1
The semantics of a VAS is given by an LTS T A = S, I, Σ, ∆ where:

• a state s is a valuation of counters in Γ, i.e., s : Γ → N; thus S = 2 Γ→N , • the initial set of states I is a set of valuations for the counters in Γ, • the set of labels Σ is defined by D,

• the set of transitions ∆ = {s ℓ -→ s ′ | s, s ′ : Γ → N, ℓ = (µ, τ), and s ≥ µ ∧ s ′ = s + τ}.

V. ABSTRACTING POPULATION PROTOCOLS

We now formalize the abstraction of the PP model into a VAS model and present the main property of this abstraction which allows us to perform the verification of the CS property on the abstract model. • Γ = {c q | q ∈ Q}, and

• D = {(µ t , τ t ) | δ ∋ t = q 1 q 2 -
→ q 3 q 4 and µ t = {q → (q 1 ≡ q) + (q 2 ≡ q) }, τ t = {q → -(q 1 ≡ q) -(q 2 ≡ q) +(q 3 ≡ q) + (q 4 ≡ q) } } where the expression (q ≡ q ′ ) returns 1 if q is the same as q ′ and 0 otherwise. We denote by α(t) the couple (µ t , τ t ) defined above. The abstraction function α comes with a concretization function γ which maps back the VAS model into the PP model. Figure 3 illustrates the relation between the mapping α and γ.

From the definition of α, it results that the LTS of a VAS A abstracting a PP model P has a number of states in O(k |Q|-1 ) while the LTS of P has O(|Q| k ) states. But k varies for a fixed Q and, in general, k ≫ |Q|. Thus, the number of states T A is smaller that the one of T P when A = α(P ). Also, α translates PP transitions in U(δ) into transitions of VAS with a not null translation component.

A couple of mappings (α, γ) induces a couple (α LT S , γ LT S ) defined on the LTS representing the semantics of PP and VAS models. Let T P = S, I, Σ, ∆ be the LTS of some PP P , and T A = S # , I # , Σ # , ∆ # be the LTS of A = α(P ). The α LT S mapping is defined by two mappings α S : S → S # and α Σ : Σ → Σ # such that:

α S (s) = {q → |{i | s(i) = q}|} α Σ ({i, j},t) = α(t) = (µ t , τ t )
The properties of the induced mapping α LT S are described by the following proposition:

Proposition 1: The components of T P and T A with A = α(P ) satisfy the following identities:

S # = s∈S {α S (s)}, Σ # = ({i, j},t)∈Σ {α Σ (t)} I # = s∈I {α S (s)}, ∆ # = s {i, j},t ----→s ′ ∈∆ {α S (s) α Σ (t) ----→ α S (s ′ )}
In the following, we remove the subscripts of the α and γ mappings if they are obvious from the context. Because one abstract state s # represents all the concrete states s which have the same number of agents in the same state, we have that one abstract execution E # in T A represents several concrete executions in E 1 , E 2 , . . . in T P . Figure 3 illustrates this fact. The concrete states situated at the same distance form the initial state are permutations of each other (relation represented by a dotted arrow in Figure 3). Also, the labels of transitions in the concrete executions abstracted by E # correspond in their second component t.

Using the induced abstraction on LTS, we obtain that an abstract computation C # = s # 0 , ℓ # 0 , s # 1 , ℓ # 1 , . . . stably computes a predicate f : X k → Y if there exists some i such that, for all

j ≥ i Σ q,ω(q)= f ( x) s # j (c q ) = k (1)
or equivalently, We denote by α( f ) the state predicate in equation 1. Then, C # stably computes f if it satisfies ✸✷α( f ).

Σ q,ω(q) = f ( x) s # j (c q ) = 0 (2) P ∈PP A ∈VAS α γ s # 0 s # 1 s # 2 . . . s # i . . . α(t j0 ) α(t j1 ) α(t j2 ) α(t ji ) s 0 γ s 1 γ s 2 γ . . . s i γ . . . ℓ 0 ,t j0 ℓ 1 ,t j1 ℓ 2 ,t j2 ℓ i ,t ji s ′ 0 γ s ′ 1 γ s ′ 2 γ . . . s ′ i γ . . . ℓ ′ 0 ,t j0 ℓ ′ 1 ,t j1 ℓ ′ 2 ,t j2 ℓ ′ i ,t ji

B. Properties of the abstraction

We are now ready to show that the verification of the correct stabilization property on the PP can be transferred into a verification of the same property in the abstract model. This result is given by the following theorem: Theorem 1 (Property transfer): A PP P stably computes a function f iff its VAS abstraction α(P ) stably computes

α( f ), i.e., T P satisfies GF ⇒ ✸✷ ω(s) = f ( x) iff T α(P ) satisfies GF ⇒ ✸✷α( f )
To obtain the proof of this theorem, we need the following two lemmas which relate GF computations in the LTS of the concrete and abstract models.

Proposition 2: If C is a GF computation of P , then α(C) is a GF computation in α(P ). Proof: Let C # = α(C) and let s # be a state in C # such that s # appears infinitely often. Recall that s # represents a finite set of states in C, denoted by γ(s # ). Then, there is at least one state s in γ(s # ) which appears in C infinitely often. Let s # 1 be a state neighbor of s # , i.e., s # -→ s # 1 . Since the transitions in P are anonymous and symmetric, the state s has also a neighbor s 1 such that s 1 ∈ γ(s # 1 ). From the hypothesis that C is a GF computation, we obtain that s 1 appears infinitely often in C, and then s # 1 = α(s 1 ) appears infinitely often in C # . Proposition 3: If C # is a GF computation of α(P ), then there exists C ∈ γ(C # ) s.t. C is a GF computation of P . Proof: Notice than not all computations in γ(C # ) are GF computations even if C # is GF computation. Indeed, suppose that some state s # i appears infinitely often in C # , and let z # i be its neighbor by some transition labeled by t. If s i is the i th concrete state in some C ∈ γ(C # ), i.e., s i ∈ γ(s # i ), the neighbors of s i by t are states z i such that z i ∈ γ(z # i ), and s i {u,v},t ----→ z i for some couple of agents (u, v). The computation C may not execute infinitely often some of the neighbors in γ(z # i ), but still its abstraction C # is still a GF computation. Then, we have to prove that there exists some computation C in γ(C # ) that contains infinitely often any neighbor of any state s i which appears infinitely often. To build such a sequence, we observe that the maximal number of neighbors of s i by some transition t ∈ δ is a constant ξ fixed by the population size k. Thus, we build C from C # by choosing some initial state s 0 ∈ γ(s # 0 ) and by choosing agents to meet in the transitions of C # such that all the interactions are considered by all concretizations of s # i which appears infinitely often in C. This choice is possible due to the symmetry property of transitions in P .

Proof of Theorem 1: (⇒) Suppose that any GF computation of T P satisfies ✸✷ ω(s) = f ( x) , and consider one GF computation C # in T α(P ) . From Proposition 2, there is some C which is a GF computation of P in γ(C # ). By hypothesis, C satisfies ✸✷ ω(s) = f ( x) , i.e., there exists some i such that for all j ≥ i, s j satisfies ω(s j ) = f ( x). Since C # has the same transitions than C, we obtain that for all j ≥ i, s # j (c q ) has a strictly positive value only if ω(q) = f ( x). From the definition of α( f ) (Equation 1), it follows that s # j satisfies α( f ).

(⇐) Suppose that any GF computation of T α(P ) satisfies ✸✷α( f ), and consider some GF computation C in T P . From Proposition 1, C # = α(C) is a GF computation of T α(P ) . By hypothesis, C # satisfies ✸✷α( f ), i.e., there exists some i such that for all j ≥ i, s # j satisfies α( f ). From the definition of α( f ) , it results that for any s in γ(s # j ), there is no agent with state q such that ω(q) = f ( x), which also means that ω(s j ) = f ( x). Then, in computation C, all the states s j for j ≥ i satisfies ω(s) = f ( x).

C. Verification under weak fairness

As mentioned in Section II, the verification of a property under some strong fairness (LF or GF) constraint is more expensive than the verification under the WF constraint. Thus, we study the conditions on which the WF constraint can be used for the verification of PP. The problem of using the WF constraint is that we can obtain false negatives, since the set of WF computations includes the set of GF computations. For example, a false negative is signaled for the MAJORITY protocol because the loop shown at Example 1 is not a WF computation.

The results of the model-checking under GF and WF constraints are equal for computations which are "insensitive" to the fairness constraint. We provide two sufficient conditions on the transition system T P of some PP model P to obtain this equality. These conditions constrain the strongly connected components (SCC) of the transition system of the protocol.

Condition SCC 1 : T P has only sink strongly connected components.

Condition SCC 2 : All the SCC of T P have the same output value.

Theorem 2 (Correction of verification under WF):

Let P be a PP whose T P satisfies SCC 1 or SCC 2 . If P stably computes f under the WF constraint then P stably computes f under the GF constraint.

Proof : Recall that the correct stabilization property has the form ✸✷φ. Its model-checking algorithm consists in finding a point in each execution from which φ is satisfied in each configuration.

Condition SCC 1 means that each execution has the form of a finite sequence (without loops) followed by a SCC. For both GF and WF constraint, this kind of executions is fair. The point to be searched by the model-checking algorithm is in a state before or at the entry in the sink SCC. Thus, the result of the model-checking is given by the satisfaction of φ in all the states of the sink SCC.

Condition SCC 2 allows executions that loop infinitely in some SCC. The GF constraint forces the exit from any nonsink SCC, but the WF constraint does not if the SCC has at least two transitions. However, since all SCC have the same output, the result of the model-checking is the same in both cases and equal to the satisfaction of φ by the common output value.

Notice that the constraint SCC 1 is not satisfied by a PP model P which contains only identity transitions, i.e., ones not in U(δ), because identity transitions generates self loops in T P . This check is syntactic, and can be done in constant time. For PP models where U(δ) = δ, the condition SCC 1 is not very restrictive. For example, it is satisfied by all our benchmark protocols except MAJORITY. The condition SCC 1 can be tested using classical tools on LTS [?], as well as the HSF-Spin [START_REF] Edelkamp | Directed explicit-state model checking in the validation of communication protocols[END_REF], and extension of Spin which computes the SCC of any Spin model. Notice however that this test is as expensive as testing the strong fairness.

VI. VERIFICATION USING SPIN

The Spin model-checker [START_REF] Holzmann | Spin model checker, the: primer and reference manual[END_REF] is one of the most efficient model-checkers for finite systems. The language for model description in Spin is called Promela. It provides intuitive constructs for modeling finite state automata extended with finite domain variables (e.g., finite counters). The general scheme of translation from our VAS model to the Promela code is given in Figure 4. The do loop statement contains several branches (introduced by ::) which are chosen to be executed in a non-deterministic way; a branch is executed only if its guard (expression before →) is true.

The property to be checked on the model is specified using LTL. To check an LTL formula ϕ on a model M, Spin executes the following two steps:

(1) ltl2ba: the negation of ϕ is translated into a Büchi automaton which thus accepts all the traces which are not satisfying ϕ; the size of this automaton 2 O(|ϕ|) where |ϕ| is the number of symbols in ϕ, This step is called one time for each formula to be checked, so its result can be used for several models M.

(2) check: the product of the Büchi automaton with the input model M is done on-the-fly. ϖ is satisfied by M if there are no cycles including an an accepting state of the Büchi automaton. Otherwise, a trace of M not respecting the formula ϕ is provided. The time complexity of this step is then 2 O(|ϕ|) O(|M|), where |M| is the sum of transitions and states in M.

The exponential size of the Büchi automaton generated by the ltl2ba step prevents us to encode the GF constraint in the LTL formula in order to select only GF computations. Indeed, the size of an LTL formula encoding the GF constraint is linear in the size of the LTS of M. To give an idea of the explosion of this step, we provide in Table II the execution times of the ltl2ba step for formulas encoding the LF and the WF constraints. These experiments have been done on an Intel i686 double core with 4GHZ and memory limited to 512MB. All experiments taking more than 30 minutes have been stopped (entries "-" in the table).

Fortunately, Spin has an option to select WF computations during the check step. This selection is done very efficiently by Spin. From the Theorem ??, the result of model-checking with Spin under the WF constraint is valid for models satisfying SCC 1 or SCC 2 . Thus, we have also to do the check on one of these conditions.

A. Benchmark

The PP models in our benchmark are given on Table I together with the characteristics of their VAS model: the number of transitions (|D|) and the number of counters (|Γ|). The majority and threshold-n protocols are our running examples. The broadcast protocol computes the OR predicate on the input values. The mod3 protocol instantiate the one given in [START_REF] Angluin | The computational power of population protocols[END_REF] to compute the predicate testing that the quotient by 3 of the number of agents with input 1 is equal to the number of agents in state 1. The flock-n protocols are the "flock of birds" protocol given in [START_REF] Chatzigiannakis | Algorithmic verification of population protocols[END_REF] which computes the same predicate that threshold-n but in a different manner.

protocol computes |D| |Γ| majority |ι(1)| > |ι(0)| 4 4 mod3 [ |ι(1)| / 3 ] = |ω(1)| 3 4 broadcast [ |ι(1)| ≥ 1 ] 1 2 flock-2 |ι(1)| ≥ 2 3 3 flock-4 |ι (1)

B. Experimental results

Table III presents the experimental results obtained for the phase check on several instances of PP in our benchmark. The execution memory is limited to 512MB. The worse execution times observed for different inputs are given. The entries annotated with (+) (resp. (-)) denote a positive (resp. negative) answer of the checking phase. The column SCC 1 indicates whether the model satisfies condition SCC 1 .

As expected, the MAJORITY protocol does not satisfies the condition SCC 1 and Spin returns a negative result. For k = 3, the execution obtained as witness for the negative result is the one shown in Figure 1. Notice that all other protocols verify condition SCC 1 and the checking times are very performant.

VII. PROBABILISTIC VERIFICATION

A. Probabilistic model

An alternative way to introduce fairness constraints in the PP model is to consider probabilistic schedulers. Angluin et al. consider this alternative in [START_REF] Angluin | A simple population protocol for fast robust approximate majority[END_REF] and define the notion of probabilistic computation for a PP model. Definition 8 (Probabilistic computation): A population protocol P = Q, X,Y, ι, ω of population size N computes a predicate f : X N -→ {0, 1} in w steps with error probability ε if for all x ∈ f -1 ({0, 1}), the state s reached after w steps satisfies the following properties with probability 1 -ε:

1) all agents agree on the correct output, i.e., ∀i ∈ {1, . . . , N} f ( x) = ω(s(π i )), and 2) the previous property is also true for every state reachable from s.

Definition 9 (Probability of a computation): For a given scheduler, the probability that a PP P = Q, X,Y, ι, ω computes a predicate f on input x is the probability of all computations starting in state ι( x) that stabilize with output f ( x).

Let us introduce now two different probabilistic schedulers. The simplest and more natural one has been considered in [START_REF] Angluin | Computation in networks of passively mobile finite-state sensors[END_REF] and schedules ordered pairs of agents.

Definition 10: A random pairing (RP) scheduler choses, in each state of a computation, the ordered pair of agents which interact in a random independently and uniformly manner from all ordered pairs of agents.

The modeling of the RP scheduler needs the LTS of a PP because it makes its decision based on the state. It cannot be modeled on an abstract VAS model because the VAS model collapses transitions whose set of agents is the same. But working on the LTS of PP suffers from the state explosion problem when the population size increases. For example, the MAJORITY protocol can only be verified with existing probabilistic model-checkers for population sizes less than 30 agents. Therefore, we define a scheduler that can be used in conjunction with the counting abstraction because its decision is based on the transitions of the PP which are also transitions of the abstract VAS model.

Definition 11:

The random ruling (RR) scheduler choses, in each state of the computation, one enabled transition t ∈ δ of the PP in a random uniformly and independently manner.

Note that for each state s, the set of enabled transitions is a subset of δ. Therefore, the probability of each transition to be chosen by the scheduler is greater than a strictly positive constant (1/|δ|). This fact allows us to obtain the following property as a consequence of the ergodic theorem on finite discrete Markov chains: Proposition 4: A population protocol P computes a predicate f almost surely with a random pairing scheduler if and only if P computes f almost surely with a random ruling scheduler.

Notice that, for the non-probabilistic model, the LF constraint, which also considers the firing of transitions of P , is not sufficient to verify the correct stabilization property for all PP, as illustrated in Example 1. In the probabilistic case, considering states (random pairing) or transitions (random ruling) is equivalent for almost sure computations of predicates.

In fact, as in the non-probabilistic case, to verify that a predicate f is computed almost surely when probabilities are strictly positive consists in searching the closed strongly connected components in the transition system. Since the LTS in the probabilistic and non-probabilistic cases are the same, we obtain the following result: Proposition 5: A PP P stably computes a predicate f ⇔ P computes f almost surely with a RP scheduler ⇔ P computes f almost surely with a RR scheduler.

Hence, we can use probabilistic model-checkers to verify that P stably computes a predicate f even for a non randomized model of population protocols. Indeed, we have to verify that P computes f almost surely with some scheduler assigning strictly positive probabilities to transitions.

B. Experiments with PRISM

Prism [START_REF] Rutten | Mathematical Techniques for Analyzing Concurrent and Probabilistic Systems[END_REF], [START_REF] Kwiatkowska | Probabilistic symbolic model checking with PRISM: A hybrid approach[END_REF] is a tool for formal modeling and analysis by model-checking of systems which exhibit random or probabilistic behavior. Models are described using the PRISM language, a simple, state-based language. The property specification language incorporates the several temporal logics including LTL, as well as extensions for quantitative specifications and costs/rewards.

We encode in Prism the abstract VAS model used for Spin, without probabilities, but specified as a discrete Markov chain. Then, Prism associates to each transition enabled at some state a uniform distribution of probabilities. Based on the results of the previous section, we model-check the VAS model against the property ✸✷α( f ) with probability at least one, i.e., almost surely.

Table IV presents the experimental results obtained with Prism. The time reported at each entry is the maximal time observed on the different inputs for a given population size. Entries marked by '-' denote a timeout (fixed to 30 minutes) of the experiment. Entries marked by (*) correspond to a partial experiment where some inputs of the protocol produce a timeout.

VIII. EXTENSIONS

The previous sections demonstrate that standard verification tools may help to prove that PP stably computes some predicate. In fact, these tools may be used to verify other properties of population protocols. Indeed, it could be interesting to verify more general properties on states that can be expressed by LTL formulas. For example, such properties are (1) all agents are infinitely often in some state, or (2) for PP with privacy [START_REF] Delporte-Gallet | Secretive birds: Privacy in population protocols[END_REF], to verify that the privacy of computations is ensured.

A probabilistic model-checker like Prism allows to compute the expected average time of convergence, i.e., the number of interactions that occur before convergence. We illustrate this last extension of use to compute in a random pairing scheduler model the average time of convergence for two classical problems:

• the epidemic broadcast: some information broadcasted by some agents needs to reach all agents; • the leader election: starting from any numbers of leaders, the protocol converges to a state with at most one leader. Notice that for an appropriate choice of probabilities for transitions, the computation of the average time of convergence can be done on the abstract VAS model. Figure 5 (resp. Figure 6) gives the Prism code of the epidemic broadcast (resp. the leader election). In both cases, we specify probabilities for transitions corresponding to the random pairing scheduler. For example, for the epidemic broadcast (Figure 5), this probability is 2*IN*(N-IN))/(N*(N-1)), where IN is the number of infected agents (i.e agents already having the information) and N is the population size.

Angluin et al. [START_REF] Angluin | Fast computation by population protocols with a leader[END_REF] prove that, for a population size N, the expected average convergence time is:

• for the epidemic broadcast, Θ(N log(N)), and • for the leader election, Θ(N 2 ). Using Prism and its reward computation option, we obtain these results experimentally, as shown on Figures 7 and8

IX. CONCLUSION

We show in this paper that the counting abstraction is very useful to deal with the verification by model-checking of PP on complete graphs. Indeed, it allows to use the existing model-checking tools without facing the state explosion problem for population sizes of order of 10 3 . Nevertheless, the fairness constraint of PP for stably computes property is very difficult to verify at least when using model checkers like Spin. We provide sufficient conditions to replace this fairness constraint by the weak fairness and thus to obtain scalability of verification. Moreover, since the fairness condition can be replaced by randomization, probabilistic model checkers like Prism allow to directly verify the stably computation property with good time performances.

Therefore, we recommend to use the existing modelcheckers for the verification of finite instances of PP. We show that such an experience may bring to more interesting ones, e.g., the computation of the expected average time of convergence.

Finite model checkers allow only verification of PP with a fixed population size and the main extension of this work would be to deal with unbounded population protocols.

Definition 5 :

 5 A population protocol (PP) is specified as a six-tuple P = (Q, X,Y, ι, ω, δ) which contains a finite set Q of possible agent states, an input assignment ι : X → Q, an output assignment ω : Q → Y , and a transition relation δ ⊆ Q 4 .

Definition 6 :Figure 2 .

 62 Figure 2. The Majority PP for k = 3 abstracted into a Petri net.

A. Abstraction function Definition 7 :

 7 The counting abstraction α maps a PP model P = Q, X, ι, ω, δ into a VAS model A = Γ, D as follows:

Figure 3 .

 3 Figure 3. Abstraction and concretization mappings.

  do :: guard(φ 1 ) → trans(φ 1 ) ... ... ... :: guard(φ m ) → trans(φ m ) od

Figure 4 .

 4 Figure 4. Promela code scheme for VAS model of PP.

  const int N; // Population size global IN: [1..N] init 1; // number of infected nodes // -transition T1: IN, NI-> IN,IN module T1 [] (IN >= 1) & (N-IN >= 1) -> (2*IN*(N-IN))/(N*(N-1)): (IN' = IN + 1) + 1-(2*IN*(N-IN))/(N*(N-1)): true; endmodule

Figure 5 .

 5 Figure 5. Prism code for epidemic broadcast.

  const int N; //Population size const int L=N; // Number of initial leaders global CL: [1..N] init L; // -transition T1: L,L-> L,N module T1 [] (CL >= 2 ) -> (CL*(CL-1))/(N*(N-1)): (CL' = CL -1 ) + 1-(CL*(CL-1))/(N*(N-1)): true; endmodule

Figure 6 .

 6 Figure 6. Prism code for leader election.

  .

Figure 7 .

 7 Figure 7. Expected average convergence time for epidemic broadcast.

Figure 8 .

 8 Figure 8. Expected average convergence time for leader election.

Table I TIME

 I PERFORMANCES (IN SEC.) FOR THE STEP ltl2ba.

	protocol	|D|	WF	LF
	majority	4	0.057	0.021
	mod3	3	0.087	0.021
	broadcast	4	0.057	0.021
	flock-2	3	0.043	0.021
	flock-4	9	0.377	0.091
	flock-8	28	-	-
	threshold-2	3	0.033	0.030
	threshold-4	7	-210.515
	threshold-8 14	-	-

Table III PERFORMANCES

 III (IN SEC.) FOR CHECKING CORRECT STABILIZATION.

Table IV PERFORMANCES

 IV (IN SEC.) FOR MODEL-CHECKING WITH PRISM.

	protocol		population size k
		11	101	1001
	majority	0.140	4.824	(*) 27.202
	mod3	0.228	0.570	(*) 70.486
	broadcast	0.137	0.385	14.982
	flock-2	0.284	14.933	(*) 35.398
	flock-4	0.661 (*) 122.240 (*)13205.884
	flock-8	0.561	(*) 4.891 (*) 1030.536
	threshold-2 0.167	10.242	(*) 24.455
	threshold-4 0.531	(*) 34.441	-
	threshold-8 0.719	-	-


Said Msaidie

Axel Ducourneau

Gilles Boëtsch

Guy Longepied

Kassim Papa

Claude Allibert

Ali Ahmed Yahaya

Jacques Chiaroni

Ph.D Michael John Mitchell
email: michael.mitchell@univmed.fr

Gilles Boetsch

Genetic diversity on the Comoros Islands shows early seafaring as a major determinant of human biocultural evolution in the Western Indian Ocean

Keywords: Y chromosome, mitochondrion, Indian Ocean, East Africa, Comoros

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The Indian Ocean can be considered as a closed sea, an afro-Asiatic Mediterranean, [START_REF] Beaujard | The Indian Ocean in Eurasian and African World-Systems Before the Sixteenth Century[END_REF][START_REF] Chaudhuri | Trade and Civilisation in the Indian Ocean: An Economic History from the Rise of Islam to 1750[END_REF] around which populations have migrated and mixed. In contrast to the Atlantic Ocean, which was a formidable natural barrier to East-West migration, the Indian Ocean with its seasonal monsoon winds favoured such exchanges, and most of the early trade routes were maritime. The Comoros archipelago is situated in the western Indian Ocean, midway between the island of Madagascar and the coast of East Africa at the northern end of the Mozambique Channel. The archipelago is composed of four main islands Grand Comore (Ngazidja), Anjouan (Ndzuani), Mohéli (Mwali) and Mayotte (Maore).

The settlement of the four islands was an integral part of migration within the Indian Ocean, since they represent a potential maritime crossroads, and juncture, between Bantu African, Middle Eastern and Southeast Asian (SEA) spheres of influence. The modern Comorian population is the result of a long-term process of biocultural admixture, mainly related to ancient trade and colonisation in the Indian Ocean.

The Comoros and Madagascar share obvious signs of SEA influence including the cultivation of rice (phased out during 20 th century), bananas and coconuts, and the use of outrigger canoes. Evidence from plant translocation suggests a migration from SEA 1 500 Years Before Present (YBP). [START_REF] Beaujard | The Indian Ocean in Eurasian and African World-Systems Before the Sixteenth Century[END_REF][START_REF] Allibert | KOMR" dans l'océan indien[END_REF]4 Clear genetic evidence for the SEA influence has been found on neighbouring Madagascar. [START_REF] Hurles | The dual origin of the malagasy in island southeast Asia and East Africa: evidence from maternal and paternal lineages[END_REF][START_REF] Regueiro | Austronesian genetic signature in East African Madagascar and Polynesia[END_REF][START_REF] Soodyall | Polynesian' mtDNA in the Malagasy[END_REF][START_REF] Tofanelli | On the origins and admixture of Malagasy: new evidence from high resolution analyses of paternal and maternal lineages[END_REF] Based on Y chromosome and mitochondrial variation, ethnic groups with the strongest SEA biocultural features on Madagascar were estimated to have approximately 50% SEA ancestry. [START_REF] Hurles | The dual origin of the malagasy in island southeast Asia and East Africa: evidence from maternal and paternal lineages[END_REF][START_REF] Tofanelli | On the origins and admixture of Malagasy: new evidence from high resolution analyses of paternal and maternal lineages[END_REF] In contrast to Madagascar where the language, Malagasy, is an Austronesian language with origins in SEA, the languages spoken on the Comoros are of Bantu origin. They are distinct from, but have close affinity to, Swahili, both branching from the precursor Sabaki language, 1 000 -2 000 YBP. [START_REF] Nurse | Swahili and Sabaki: A Linguistic History[END_REF] The cultural contributions of Middle Eastern civilisation are equally evident on the Islands. By 2 000 YBP, a thriving commercial maritime network already existed, extending from the Middle East to India, and as far South as Tanzania on the East African coast. The name "Comoros" is from the Arabic Kmr, meaning "light in the sky". [START_REF] Allibert | KOMR" dans l'océan indien[END_REF] From 1 300 YBP, the Comoros archipelago served as a stepping-stone, for Middle Eastern traders operating along the East African coast, and for Southeast Asian traders travelling to Madagascar and the East African coast. [START_REF] Allibert | Austronesian Migration and the Establisment of the Malagasy Civilization: Contrasted Readings in Linguistics, Archaeology, Genetics and Cultural Anthropology[END_REF][START_REF] Allibert | The Early Pre-Islamic History of the Comores Islands: Links with Madagascar and Africa[END_REF] By 1 000 YBP, the Shirazi, traders with origins in the Persian city of Shiraz in present day Iran, had established themselves on the island of Kilwa. The Shirazi were responsible for the generalisation of Islam on the Swahili coast by 500 YBP. They had built mosques on Kilwa, Zanzibar and Anjouan by 800 YBP. [START_REF] Allibert | La chronique d'Anjouan par Said Ahmed Zaki[END_REF] Islam remains the religion of the Islands today.

An unambiguous genetic signal from the Middle East has not, however, been detected in East Africa further south than Ethiopia, [START_REF] Luis | The Levant versus the Horn of Africa: Evidence for Bidirectional Corridors of Human Migrations[END_REF][START_REF] Semino | Origin, diffusion, and differentiation of Ychromosome haplogroups E and J: inferences on the neolithization of Europe and later migratory events in the Mediterranean area[END_REF] or in the ethnic groups sampled on Madagascar. [START_REF] Hurles | The dual origin of the malagasy in island southeast Asia and East Africa: evidence from maternal and paternal lineages[END_REF][START_REF] Tofanelli | On the origins and admixture of Malagasy: new evidence from high resolution analyses of paternal and maternal lineages[END_REF] The Lemba people of South Africa, carrying a putative semitic Y chromosome, currently provide the only evidence for gene flow from the Middle East into southern Africa. [START_REF] Spurdle | The origins of the Lemba "Black Jews" of southern Africa: evidence from p12F2 and other Y-chromosome markers[END_REF] The alleles of some autosomal genes found in the ex-patriot Comorian population living in Marseilles indicate a genetic contribution from Western Eurasia, [START_REF] Badens | Molecular basis of haemoglobinopathies and G6PD deficiency in the Comorian population[END_REF][START_REF] Chiaroni | Adsorption of autoantibodies in the presence of LISS to detect alloantibodies underlying warm autoantibodies[END_REF][START_REF] Gibert | HLA-DRB1 frequencies of the Comorian population and their genetic affinities with Sub-Saharan African and Indian Oceanian populations[END_REF] but the populations living on the Comoros have until now not been studied.

In this context, the peopling of the Comoros is evidently integral to the movements of men and women across the entire Indian Ocean. To gain insights into this process, we therefore determined the Y chromosomal and mitochondrial genetic variation on the three Bantu-speaking islands of the Comoros Republic.

Material and methods

Sample group

In February and March 2006, we obtained blood samples from 577 unrelated Comorian men (n=381) and women (n=196). We sampled the populations of three of the four islands of the Comoros archipelago (Grand Comore -170 men, 67 women, Anjouan -104 men, 69 women and Moheli -107 men, 60 women). In 2006, this represented approximately 0.1% of the total Comoros population of 690 000 people. Blood was collected in EDTA vacutainer tubes and DNA extracted using the salting-out method. [START_REF] Miller | A simple salting out procedure for extracting DNA from human nucleated cells[END_REF] Samples were collected from multiple towns and villages on each island (supplemental Figure 1). Recruitment was achieved through contacts established by medical personnel who originated from each community sampled. Each donor included in the study had four grandparents who were born on the same island and were native speakers of the island's language (Shingazidja, Shindzuani or Shimwali). Informed consent was obtained from all participants.

Y chromosome haplogroups and haplotypes

We typed 68 binary polymorphisms mainly by PCR-RFLP (Figure 1 and

Mitochondrial haplogroups

We typed 31 coding region polymorphisms in the mitochondrial genome mainly by PCR RFLP (Figure 3, Supplementary Table 3). We also sequenced a 501 bp fragment including HVS-I from all M, N(xR) and R samples (Supplementary Table 4 and GenBank: HM565257-HM565275). Choice of markers and branch designations was based on published data [START_REF] Behar | The dawn of human matrilineal diversity[END_REF] and trees presented at http://www.ianlogan.co.uk.

Data analysis

The genetic structure (haplogroup number, haplogroup diversity, population differentiation, Fst and Rst) of the study population sample was analysed using the ARLEQUIN package v. 3.01, 21 and phylogenetic comparisons of Y STR haplotypes and mitochondrial SNP haplogroups were examined by Multidimensional Scaling (MDS) based, respectively, on Rst or Fst distances, [START_REF] Reynolds | Estimation of the Coancestry Coefficient: Basis for a Short-Term Genetic Distance[END_REF][START_REF] Slatkin | A measure of population subdivision based on microsatellite allele frequencies[END_REF] using SPSS 10.00 software. Admixture fractions were estimated using ADMIX 2.0. [START_REF] Dupanloup | Inferring admixture proportions from molecular data: extension to any number of parental populations[END_REF] Based on our haplogroup and MDS analyses, and historical and linguistic data, we chose Borneo, Iran and East Africa (Y: Kenya and Tanzania and mitochondrial: Mozambique) as the most likely parent populations from published data. [START_REF] Hurles | The dual origin of the malagasy in island southeast Asia and East Africa: evidence from maternal and paternal lineages[END_REF][START_REF] Luis | The Levant versus the Horn of Africa: Evidence for Bidirectional Corridors of Human Migrations[END_REF][START_REF] Hill | A mitochondrial stratigraphy for island southeast Asia[END_REF][START_REF] Metspalu | Most of the extant mtDNA boundaries in south and southwest Asia were likely shaped during the initial settlement of Eurasia by anatomically modern humans[END_REF][START_REF] Regueiro | Iran: tricontinental nexus for Y-chromosome driven migration[END_REF][START_REF] Salas | The making of the African mtDNA landscape[END_REF]

Results

Y chromosome diversity on the Comoros

We analysed 381 Y chromosomes from the Comoros and identified 28 distinct haplogroups belonging to 11 of the 20 major clades of the Y chromosome tree as shown below (Figure 1). [START_REF] Karafet | New binary polymorphisms reshape and increase resolution of the human Y chromosomal haplogroup tree[END_REF] These fall into four groups, based on the geographical distribution of haplogroups around the Indian Ocean: sub-Saharan African 59.6%; Western and Southern Asia 29.7%; Southeast Asia 6% and uncertain origin 4.7%. Four clades, E, J, O and R, have frequencies greater than 5% and represent 87.4% of the sample.

The paragroups C*(xC1-5), F*(xM282,M427), J* and K*(xLMNOPQRST), cannot be assigned an origin with certainty. Nevertheless the high frequencies of C*-M216 (Borneo -2.5-25%) and K* (2-30%), [START_REF] Hurles | The dual origin of the malagasy in island southeast Asia and East Africa: evidence from maternal and paternal lineages[END_REF][START_REF] Kayser | peopling of Indonesia: genetic contributions from pre-neolithic hunter-gatherers, Austronesian farmers, and Indian traders[END_REF] in SEA, make an SEA origin probable. J* has been found in Bali (1.5%) but also on the island of Soqotra (71%) situated in the Gulf of Oman between Somalia and Yemen. 31,32 F*(xM282,M427) has been found mainly in the Indian subcontinent. [START_REF] Sengupta | Polarity and temporality of highresolution Y-chromosome distributions in India identify both indigenous and exogenous expansions and reveal minor genetic influence of Central Asian pastoralists[END_REF] A West or Southwest Asian origin is therefore more likely for the F* and J* chromosomes.

Y STR analysis revealed a generally high variance (Table 1), which coupled with the large number of Y haplogroups, suggests that genetic drift has not drastically reduced genetic diversity on the Comoros Islands.

Sub-Saharan African Y chromosomes

The most common Comorian haplogroups, E1b1-M2 (41%) and E2-M90 (14%), are those that are frequent in sub-Saharan Africa. [START_REF] Luis | The Levant versus the Horn of Africa: Evidence for Bidirectional Corridors of Human Migrations[END_REF][START_REF] Cruciani | A back migration from Asia to sub-Saharan Africa is supported by high-resolution analysis of human Y-chromosome haplotypes[END_REF][START_REF] Passarino | Different genetic components in the Ethiopian population, identified by mtDNA and Y-chromosome polymorphisms[END_REF][START_REF] Scozzari | Combined use of biallelic and microsatellite Y-chromosome polymorphisms to infer affinities among African populations[END_REF][START_REF] Underhill | Y chromosome sequence variation and the history of human populations[END_REF] They are present, respectively, at 56% and 6.4%, on Madagascar. [START_REF] Tofanelli | On the origins and admixture of Malagasy: new evidence from high resolution analyses of paternal and maternal lineages[END_REF] Two haplogroups were identified under E1b1-M2, derived for markers M191 (22%) and U209 (9%). The haplogroup E1b1a-M191 has been found in east and west sub-Saharan Africa, 19% in Tanzania and 57% in Benin. [START_REF] Luis | The Levant versus the Horn of Africa: Evidence for Bidirectional Corridors of Human Migrations[END_REF] The marker U209 was identified in Afro-Americans, [START_REF] Sims | Sub-populations within the major European and African derived haplogroups R1b3 and E3a are differentiated by previously phylogenetically undefined Y-SNPs[END_REF] and has not, until now, been tested for in African populations.

The low incidence of E-M293 (0.8%) and A-M91 (0%) on the Comoros contrasts strongly with the frequency of these haplogroups in East African populations. E-M293 is found mainly in East Africa, Kenya and Tanzania (18%). [START_REF] Henn | Y-chromosomal evidence of a pastoralist migration through Tanzania to southern Africa[END_REF] Furthermore, on the African mainland M293 chromosomes carry either 10, or 13 and more repeats at the DYS389I STR locus, [START_REF] Henn | Y-chromosomal evidence of a pastoralist migration through Tanzania to southern Africa[END_REF] while, on the Comoros, they have 12 repeats. Haplogroup A has a frequency of 14% in Kenyan Bantu and 7% in Tanzania. [START_REF] Luis | The Levant versus the Horn of Africa: Evidence for Bidirectional Corridors of Human Migrations[END_REF] Other haplogroups of likely sub-Saharan African origin on the Comoros are E-SRY 4064 (xM2,M35,M75) (1.3%) and B2a (1.6%). B2a has a low frequency in southern Iran and Qatar, [START_REF] Regueiro | Iran: tricontinental nexus for Y-chromosome driven migration[END_REF]39 but this is thought to be a consequence of the Arab slave trade. We therefore treat B2a as an African chromosome in this study.

Y chromosomes from around the Arabian Sea.

The northern Y chromosomes on the Comoros, E-V22, E-M123, F*(xF2, GHIJK), G2a, I, J1, J2, L1, Q1a3, R1*, R1a*, R1a1, and R2 (29.7%), make up a diverse group. G2a, J1 and J2 (16.5%) are thought to have originated in the Middle East. [START_REF] Semino | Origin, diffusion, and differentiation of Ychromosome haplogroups E and J: inferences on the neolithization of Europe and later migratory events in the Mediterranean area[END_REF][START_REF] Cinnioglu | Excavating Y-chromosome haplotype strata in Anatolia[END_REF] J1-M267 has mainly spread south and west into the Arabic Peninsula, and into North and Northeast Africa, while J2-M172 lineages have expanded north into Europe and east into Asia. [START_REF] Luis | The Levant versus the Horn of Africa: Evidence for Bidirectional Corridors of Human Migrations[END_REF][START_REF] Semino | Origin, diffusion, and differentiation of Ychromosome haplogroups E and J: inferences on the neolithization of Europe and later migratory events in the Mediterranean area[END_REF]39,[START_REF] Arredi | A predominantly neolithic origin for Ychromosomal DNA variation in North Africa[END_REF][START_REF] Hassan | Y-chromosome variation among Sudanese: restricted gene flow, concordance with language, geography, and history[END_REF][START_REF] Tofanelli | J1-M267 Y lineage marks climate-driven prehistorical human displacements[END_REF] The M78 subclade, E-V22, and E-M123 are believed to have originated in Northeast Africa, with E-V22 spreading to the west of North Africa and to the Arabic peninsula by the Levantine corridor (United Arab Emirates 6.7%), 39,44 while M123 spread mainly to the East (Yemen 8%; Oman 12%; Turkey 5.5%; Iran 1%). [START_REF] Luis | The Levant versus the Horn of Africa: Evidence for Bidirectional Corridors of Human Migrations[END_REF][START_REF] Regueiro | Iran: tricontinental nexus for Y-chromosome driven migration[END_REF]39,[START_REF] Cinnioglu | Excavating Y-chromosome haplotype strata in Anatolia[END_REF] In contrast the haplogroups L1, Q1a3, R1, R1a, R1a1 and R2 (10.5%), are thought to be of Central or Southern Asian origin and describe clines of decreasing frequency from India and Pakistan towards the Middle East. [START_REF] Sengupta | Polarity and temporality of highresolution Y-chromosome distributions in India identify both indigenous and exogenous expansions and reveal minor genetic influence of Central Asian pastoralists[END_REF] A comparison of the relative incidences of E-M78(V22), E-M123, G, J, L, Q and R on the Comoros with populations around the Arabian Sea shows greatest similarities with Southern Iran and, to a lesser extent, Turkey (Supplementary Figure 2). [START_REF] Regueiro | Iran: tricontinental nexus for Y-chromosome driven migration[END_REF][START_REF] Cinnioglu | Excavating Y-chromosome haplotype strata in Anatolia[END_REF] The higher affinity to South Iran is also evident in the MDS analysis with the Comoros Y-STR data for the E-V22, E-M123, G, J, L Q and R haplogroups (Figure 2a). In the MDS, Comoros shows greatest affinity with UAE and South Iran. Southern Iran is the site of the first towns to develop in the Southern Middle East 2 000-3 000 years ago (Supplementary Figure 2).

A possible source of the Northern Y chromosomes is therefore the Shirazi traders from Southern Iran who established trading posts on the Comoros by 800 YBP. [START_REF] Allibert | La chronique d'Anjouan par Said Ahmed Zaki[END_REF] It has previously been estimated that, at 9 Y-STR loci, 0-1 mutation will most likely separate the descendants of a single Y chromosome haplotype after 40 generations (1 000 -1 200 years). [START_REF] Capelli | Moors and Saracens in Europe: estimating the medieval North African male legacy in southern Europe[END_REF]47 Compatible with a Shirazi origin, we found that, at 9 Y-STR loci (DYS19, 389AB, 389CD, 390-393, 438 and 439), 42% of the Comoros Northern chromosomes differ by 0-1 mutation from chromosomes in Southern Iran. [START_REF] Alshamali | Local population structure in Arabian Peninsula revealed by Y-STR diversity[END_REF][START_REF] Roewer | A Y-STR database of Iranian and Azerbaijanian minority populations[END_REF]

Southeast Asian Y chromosomes

We found the O1 lineage (6%) in the Comoros sample, providing genetic evidence for an SEA influence. Haplogroup O has been found at highest frequencies in East Asia and Island Southeast Asia. [START_REF] Capelli | A predominantly indigenous paternal heritage for the Austronesian-speaking peoples of insular Southeast Asia and Oceania[END_REF][START_REF] Kayser | Independent histories of human Y chromosomes from Melanesia and Australia[END_REF] All but one of the Comorian O1 chromosomes are O1a-M50 (5.8%). The O1a-M50 Y chromosome has its highest incidence in SEA: Borneo (10-20%), Sulawesi (4%), Taiwanese aborigines (0-59%, mean 14%) and the Philippines (3-12%). [START_REF] Hurles | The dual origin of the malagasy in island southeast Asia and East Africa: evidence from maternal and paternal lineages[END_REF][START_REF] Kayser | The impact of the Austronesian expansion: evidence from mtDNA and Y chromosome diversity in the Admiralty Islands of Melanesia[END_REF][START_REF] Li | Paternal genetic affinity between Western Austronesians and Daic populations[END_REF] It has not been detected in the Middle East or the Indian subcontinent. [START_REF] Hurles | The dual origin of the malagasy in island southeast Asia and East Africa: evidence from maternal and paternal lineages[END_REF][START_REF] Regueiro | Iran: tricontinental nexus for Y-chromosome driven migration[END_REF]39,[START_REF] Sengupta | Polarity and temporality of highresolution Y-chromosome distributions in India identify both indigenous and exogenous expansions and reveal minor genetic influence of Central Asian pastoralists[END_REF] We performed an MDS with our STR data for the Y haplogroups O, C* and K* together with available STR data from candidate SEA populations (Figure 2b). The Comoros show a low affinity to the populations selected, even when C* and K* are not included (not shown), suggesting that these populations are not the source of SEA chromosomes on the Comoros.

Mitochondrial diversity on the Comoros

We have tested 577 Comorian samples for mitochondrial SNPs, and we define 9 distinct haplogroups (Figure 3). As for the Y chromosome, the majority of mitochondrial haplogroups on the Comoros, are of African origin. The haplogroups L0, L1, L2 and L3'4(xMN) compose 84.7% of the mitochondria in the Comoros sample, and their relative proportions, are most similar to profiles found in East and South East Africa. [START_REF] Behar | The dawn of human matrilineal diversity[END_REF][START_REF] Gonder | Whole-mtDNA genome sequence analysis of ancient African lineages[END_REF] The higher affinity with sub-Saharan East African populations is also evident in the MDS analysis (Figure 4a andb).

The remaining 15.3% of the Comoros sample is composed almost exclusively of haplogroups that can either be unambiguously identified as SEA (B4a1a1-PM, F3b, and M7c1c -10.6%), [START_REF] Hill | A mitochondrial stratigraphy for island southeast Asia[END_REF] or fall into the paragroup M(xD,E,M1,M2,M7) (4%) (Figure 3). The latter haplogroups are probably also originally from Southeast Asia, but of the 12 different M* HVS-I sequences on the Comoros, only two match published sequences: two M(xM7) mitochondria found on Madagascar. [START_REF] Tofanelli | On the origins and admixture of Malagasy: new evidence from high resolution analyses of paternal and maternal lineages[END_REF] We found no haplogroups that could be assigned to the Middle East.

Southeast Asian mitochondria

Of the SEA mitochondrial haplogroups present on the Comoros, F3b and M7c1c, like the Y-Hg O1-M50, each define an area of distribution that extends from Taiwan through the Philippines to Borneo 25 (Supplemental Figure 3). The MDS analysis with SEA populations shows greatest affinity to the Philippines and Borneo, although affinity is relatively weak (Figure 4c). Linguistic studies indicate Southeast Borneo to be the probable origin of the migration from SEA to Madagascar. [START_REF] Gray | Language trees support the express-train sequence of Austronesian expansion[END_REF] B4a1a1-PM (0.7%) is the major haplogroup throughout Polynesia (78%), [START_REF] Kayser | peopling of Indonesia: genetic contributions from pre-neolithic hunter-gatherers, Austronesian farmers, and Indian traders[END_REF] and on Madagascar (25%), but, within island SEA, it has not been found further West than South Borneo (1%). [START_REF] Hurles | The dual origin of the malagasy in island southeast Asia and East Africa: evidence from maternal and paternal lineages[END_REF][START_REF] Soodyall | Polynesian' mtDNA in the Malagasy[END_REF][START_REF] Tofanelli | On the origins and admixture of Malagasy: new evidence from high resolution analyses of paternal and maternal lineages[END_REF]

Male-biased gene flow from the Middle East

There are no mitochondrial lineages on the Comoros that are frequent in the Middle East (Figure 3). We have tested for, but did not find, the R haplogroups, H, J, T, U and V, or N(xR) that represent 80% of the mitochondria in Iran. 56 There is therefore striking evidence for male-biased gene flow from the Middle East to the Comoros, even if the unassigned mt-Hg M* and R* are designated as western Asian: 103/381 Y vs 27/577 mitochondria -Fisher's exact test, one sided, p<10 -22 . This is entirely consistent with male-dominated trade and religious proselytisation being the forces that drove the Middle Eastern gene flow to the Comoros. For African and SEA contributions, if Y haplogroups C* and K* are counted as SEA, the under representation of the male lineages are similar (Y to mt ratio: Africa 0.69, SEA 0.66).

An opposite female gene flow from Africa to the Middle East, is clearly evident in Yemen (34% mt-Hg L; 4% Y-Hg E-M2), Iraq and the Levant. [START_REF] Richards | Extensive female-mediated gene flow from sub-Saharan Africa into near eastern Arab populations[END_REF] However, no mt-Hg L has been found in Iran (n=712), 56 despite the presence of Y-Hg E-M2 (1.7%), [START_REF] Regueiro | Iran: tricontinental nexus for Y-chromosome driven migration[END_REF] supporting the idea that the elevated mt-Hg L frequency in the western Middle East is not exclusively a consequence of the Arab slave trade, but also of geography. [START_REF] Kivisild | Ethiopian mitochondrial DNA heritage: tracking gene flow across and around the gate of tears[END_REF]

Discussion

We reveal the Comoros population to be a genetic mosaic, the result of tripartite gene flow from the North, the East and the West. Admixture analysis of the maternal and paternal contributions indicates the gene pool to be predominately African (72%), with significant contributions from Western Asia (17%) and Southeast Asia (11%). Our study therefore provides the first unequivocal evidence that the Middle Eastern trade routes that developed along the East African coast, during the last 2 000 years, have left a genetic trace. Male and female SEA gene flow has already been described on Madagascar, in populations that speak Austronesian languages, [START_REF] Hurles | The dual origin of the malagasy in island southeast Asia and East Africa: evidence from maternal and paternal lineages[END_REF][START_REF] Tofanelli | On the origins and admixture of Malagasy: new evidence from high resolution analyses of paternal and maternal lineages[END_REF] but here we show that this extends beyond Madagascar, into African populations speaking languages from the Bantu family. This raises the question of whether the demic migration from SEA reached the East African mainland.

The frequencies of Y-Hg E-V22, E-M123, G2a, J1, J2, R1a1 and R2 in the Comorian sample are compatible with gene flow from Iran. [START_REF] Regueiro | Iran: tricontinental nexus for Y-chromosome driven migration[END_REF][START_REF] Sengupta | Polarity and temporality of highresolution Y-chromosome distributions in India identify both indigenous and exogenous expansions and reveal minor genetic influence of Central Asian pastoralists[END_REF] This concords with historical data which attests to the presence of traders from Shiraz in Iran on the Comoros, and also the Comorian's own oral traditions which recount that Shirazi princes came in ships and established colonies on the islands. On the island of Anjouan the term "Shirazi" is used to refer to someone of Middle Eastern appearance. There is historical evidence that 1 000 YBP Persian traders played an important role in trade along the East coast, and we therefore predict that an Iranian genetic signal will be detected among Swahili speakers at former Middle Eastern trading centres on the sub-Saharan East coast, such as the islands of Zanzibar and Kilwa off the coast of Tanzania.

Interestingly, there are a number of similarities between the genetic profile of the Comoros islanders and the Lemba of South Africa, a Bantu speaking people whose Semitic origins are evident at both the cultural and genetic level. [START_REF] Spurdle | The origins of the Lemba "Black Jews" of southern Africa: evidence from p12F2 and other Y-chromosome markers[END_REF][START_REF] Thomas | Y chromosomes traveling south: the cohen modal haplotype and the origins of the Lemba--the "Black Jews of Southern Africa[END_REF] The Lemba have from the centre and South of Madagascar have been studied to date. [START_REF] Hurles | The dual origin of the malagasy in island southeast Asia and East Africa: evidence from maternal and paternal lineages[END_REF][START_REF] Tofanelli | On the origins and admixture of Malagasy: new evidence from high resolution analyses of paternal and maternal lineages[END_REF] The low frequencies of E-M293 and A-M91, on both the Comoros and Madagascar, contrasts with the high frequency found in inland populations from Tanzania and Kenya, [START_REF] Luis | The Levant versus the Horn of Africa: Evidence for Bidirectional Corridors of Human Migrations[END_REF][START_REF] Henn | Y-chromosomal evidence of a pastoralist migration through Tanzania to southern Africa[END_REF] and could be characteristics of a genetic profile specific to sub-Saharan coastal East Africa.

The SEA haplogroups, shared with Madagascar, on the Comoros, are O1a-M50 for the Y chromosome and M7c1c, F3b and B4a1a1-PM for the mitochondria.

Consistent with their transit West across the Indian Ocean to the East African coast, O1a-M50, M7c1c and F3b are linked to maritime colonisation within island SEA. In contrast, B4a1a1-PM has not been found in island SEA further West than Southeast Borneo (1%) and has expanded mainly East into Polynesia but also West to Madagascar where it predominates. [START_REF] Hill | A mitochondrial stratigraphy for island southeast Asia[END_REF][START_REF] Kayser | The impact of the Austronesian expansion: evidence from mtDNA and Y chromosome diversity in the Admiralty Islands of Melanesia[END_REF] There are nevertheless several indicators that the Comoros' history of gene flow from SEA is distinct from Madagascar's: the absence of Y HgO2-M95 and the very low frequency of B4a1a1-PM, on the Comoros, the higher frequency of F3b (Comoros 8%; Madagascar 3.7%), the dissimilarity of M* HVS-I sequences and the low affinity between the Comoros' O1a-M50 chromosomes and those of Madagascar. The fact that O1a-M50, M7c1c, F3b and B4a1a1-PM have not been found at sites around the Indian Ocean, [START_REF] Metspalu | Most of the extant mtDNA boundaries in south and southwest Asia were likely shaped during the initial settlement of Eurasia by anatomically modern humans[END_REF][START_REF] Regueiro | Iran: tricontinental nexus for Y-chromosome driven migration[END_REF]39,[START_REF] Sengupta | Polarity and temporality of highresolution Y-chromosome distributions in India identify both indigenous and exogenous expansions and reveal minor genetic influence of Central Asian pastoralists[END_REF]56 outside SEA and now East Africa, is consistent with a colonising migration from SEA to East Africa directly across the Indian Ocean.

The Comoros population represents an exceptionally diverse genetic mosaic created by the complex process of human settlement around the Indian Ocean. The genealogy of our large sample is well-documented and will provide a solid base from which to explore human diversity in Madagascar and coastal sub-Saharan East Africa.

Supplemental information is available at the European Journal of Human Genetics

V22 V12 SRY4064 C3 C4 C5 C1 8a E* b1* E1b1b1a I H1 G2a G G1 b1c O1a 2 J2b2 J2b J2a2 J2a J1 O1 K* M R * R2 Q1 a1 Q1 a3 Q1 a2 O3 R1 b1 R1a 1c R1a

  high frequencies of the Middle Eastern Y chromosome HgJ-12f2a (25%), a potentially SEA Y, Hg-K(xPQR) (32%) and a Bantu Y, E-PN1 (30%) (similar to E-M2), raising the possibility that the Lemba and Comorian populations are consequences of similar demographic processes. The high-resolution genotyping of the Lemba Y chromosomes and mitochondria will elucidate this question. The Comoros and Madagascar show similarities in the paternal and maternal contribution from SEA and Africa. The absence of a strong Middle Eastern signal on Madagascar could be due to sampling bias, since Arab or Persian traders are known to have established posts on the Northwest coast of Madagascar, whereas only populations

Table 1 ). For 293 Y chromosomes, alleles of 17 short tandem repeats (STR) polymorphisms on the Y chromosome were amplified with the AmpFlSTR Yfiler PCR

 1 supplementary

	Amplification Kit (Applied Biosystems). Y STR haplotypes were determined for 15 of
	38 E-M2(xM191,U209) and 19 of 84 E-M191 chromosomes, and for all other
	chromosomes.

Table 1

 1 Variance of the principal Y haplogroups (n ≥5) on the Comoros based on 15 Ymicrosatellite loci. The markers DYS385a and b were not used in the calculation. The variance is calculated from STR data for all individuals of each haplogroup (n), except E1b1a7 (19 of 84 men) and E1b1a* (15 of 38 men).

	Haplogroup Mutation n Variance
	B2a	M150	6	1.01
	C*	M216	6	0.44
	E*	SRY4064	5	1.03
	E1b1a*	M2	15	0.46
	E1b1a7	M191	19	1.13
	E1b1a8a	U209	36	0.66
	E2	M75	14	0.75
	G2a	P15	7	0.10
	J1	M267	20	0.54
	J2a	M410	28	0.60
	J2b	M12	6	0.57
	K*	M9	10	1.16
	L1	M76	9	1.05
	O1a2	M50	22	0.30
	R1a1	M17	12	0.77
	R2	M124	9	1.18

Acknowledgements

We express our sincere gratitude to the people of the Comoros who participated in this study, and the medical personnel of the Comoros and the World Health Organisation for sample collection, and in particular Dr Said Ahamada Fazul, Dr Islam Abdallah, Mr Ismail Msaidié, Mrs Wardat Abdoulkarim, Mme Oumrati Haribou, Mr A Djazza and Mr A.I Abdallah. We thank the EFS (French Blood Transfusion Service) Alpesméditerranée for their support. We further thank Professor Nicolas Levy, Dr Catherine Badens, Caroline Lacoste, Julie di Cristofaro, Marie-Claude Bonino Fabien Ciné and Nathalie Eudes for help with genotyping. A.D. was supported by a MENRT post-graduate studentship from the French Government and a Fyssen post-doctoral fellowship. S.M. was supported by a studentship "Coopération France-Comores" given jointly from the French and the Comoros governments. This work was supported by the CNRS grant "OHLL -Origines de l'Homme, des langues et du language" and core funding from Inserm to MJM.

Titles and Legends to Figures

 [START_REF] Karafet | New binary polymorphisms reshape and increase resolution of the human Y chromosomal haplogroup tree[END_REF]Branches are labelled with the binary markers tested. Numbers without a letter represent "M" prefixed Y markers (e.g; 50=M50). [START_REF] Underhill | The phylogeography of Y chromosome binary haplotypes and the origins of modern human populations[END_REF] Putative geographic origin is indicated for each haplogroup: Af -sub-Saharan Africa, WSA -West and Southwest Asia, SEA -Southeast Asia and ? -uncertain. Frequencies of less than 5% have been rounded up or down to the nearest unit. [START_REF] Tofanelli | On the origins and admixture of Malagasy: new evidence from high resolution analyses of paternal and maternal lineages[END_REF] Oman (OMA), [START_REF] Luis | The Levant versus the Horn of Africa: Evidence for Bidirectional Corridors of Human Migrations[END_REF][START_REF] Alshamali | Local population structure in Arabian Peninsula revealed by Y-STR diversity[END_REF] Turkey (TUR), [START_REF] Cinnioglu | Excavating Y-chromosome haplotype strata in Anatolia[END_REF] North Pakistan (N-PAK), South Pakistan (S-PAK), North India (N-IND), South India (S-IND), [START_REF] Sengupta | Polarity and temporality of highresolution Y-chromosome distributions in India identify both indigenous and exogenous expansions and reveal minor genetic influence of Central Asian pastoralists[END_REF] Yemen (YEM), United Arab Emirates (UAE), Saudi Arabia (SAU), [START_REF] Alshamali | Local population structure in Arabian Peninsula revealed by Y-STR diversity[END_REF] , North Iran (N-IR), [START_REF] Roewer | A Y-STR database of Iranian and Azerbaijanian minority populations[END_REF] South Iran (S-IR), [START_REF] Alshamali | Local population structure in Arabian Peninsula revealed by Y-STR diversity[END_REF][START_REF] Roewer | A Y-STR database of Iranian and Azerbaijanian minority populations[END_REF] Malaysia (MAL), [START_REF] Chang | Haplotype diversity of 16 Ychromosomal STRs in three main ethnic populations (Malays, Chinese and Indians) in Malaysia[END_REF] Taiwan (Paiwan) (TAI), [START_REF] Wu | Y-chromosomal STRs haplotypes in the Taiwanese Paiwan population[END_REF] West Borneo (East Malaysia) (W-BOR) [START_REF] Chang | Haplotype diversity of 17 Y-chromosomal STRs in three native Sarawak populations (Iban, Bidayuh and Melanau) in East Malaysia[END_REF] and Bangladesh (BAN). 64 Populations: 1a), 1b) Comoros -this study, 2a) Tanzania, 2b) Kenya, [START_REF] Luis | The Levant versus the Horn of Africa: Evidence for Bidirectional Corridors of Human Migrations[END_REF] 3) Somalia, 66 4) Yemen, 39 5) Oman, 6) Egypt (*G-M201), [START_REF] Luis | The Levant versus the Horn of Africa: Evidence for Bidirectional Corridors of Human Migrations[END_REF] 7) Turkey (*L(xM349), [START_REF] Cinnioglu | Excavating Y-chromosome haplotype strata in Anatolia[END_REF] Madagascar, [START_REF] Hurles | The dual origin of the malagasy in island southeast Asia and East Africa: evidence from maternal and paternal lineages[END_REF][START_REF] Tofanelli | On the origins and admixture of Malagasy: new evidence from high resolution analyses of paternal and maternal lineages[END_REF] and 3) Sumatra, 4) Borneo, 5) Sulawesi, 6) Malaysia, 7) Philippines, 8) Taiwan. [START_REF] Hurles | The dual origin of the malagasy in island southeast Asia and East Africa: evidence from maternal and paternal lineages[END_REF][START_REF] Hill | A mitochondrial stratigraphy for island southeast Asia[END_REF][START_REF] Kayser | The impact of the Austronesian expansion: evidence from mtDNA and Y chromosome diversity in the Admiralty Islands of Melanesia[END_REF][START_REF] Li | Paternal genetic affinity between Western Austronesians and Daic populations[END_REF] Supplemental Table 1. Y chromosome SNP markers.

Supplemental Table 2. Y chromosome haplogroups and STR haplotypes.

Supplemental Table 3. Mitochondrial SNP markers. 4. Mitochondrial haplogroups and HVS-I sequence.

Supplemental Table