domain: stat sdv

text:
Jean Peyhardi

Fabien Laroche

Frédéric Mortier

Pólya-splitting distributions as stationary solutions of multivariate birth-death processes under extended neutral theory

Keywords: Species diversity, Ecological communities, Joint Species Distribution Model, Neutral theory, Multivariate birth-death jump processes, Stationary Distributions 62H05, 60J74, 92D40

Multivariate count distributions are crucial for the inference of ecological processes underpinning biodiversity. In particular, neutral theory provides useful null distributions allowing the evaluation of adaptation or natural selection. In this paper, we build a broader family of multivariate distributions: the Polya-splitting distributions. We show that they emerge naturally as stationary distributions of a multivariate birthdeath process. This family of distributions is a consistent extension of non-zero sum neutral models based on a master equation approach. It allows considering both total abundance of the community and relative abundances of species. We emphasize that this family is large enough to encompass various dependence structures among species. We also introduce the strong closure under addition property that can be useful to generate nested multi-level dependence structures. Although all Pólya splitting distributions do not share this property, we provide numerous example verifying it. They include the previously known example with independent species, and also new ones with alternative dependence structures. Overall, we advocate that Polya-splitting distribution should become a part of the classic toolbox for the analysis of multivariate count data in ecology, providing alternative approaches to joint species distribution framework. Comparatively, our approach allows to model dependencies between species at the observation level, while the classical JSDM's model dependencies at the latent process strata.

Introduction

Understanding the processes that shape the biodiversity of ecological communities is a major question in the context of global changes. A large fraction of empirical research about understanding metacommunities relies on pattern-to-process approaches, i.e. detecting the signature of processes through a statistical analysis of the observed distribution of species in space and time. The success of pattern-to-process approaches relies on building an appropriate null hypothesis, where some target process is nullified, and testing whether observations deviate from it.

For 20 years, a family of neutral models, inspired from population genetics, have been proposed as a baseline generating null hypotheses to investigate the effect of species ecological niches. The neutrality assumption consists in assuming that all individuals are ecologically equivalent irrespective of their species, genotypes, etc. hence cancelling any effect associated to ecological niches. One of the most famous example of neutral model is the model introduced by [START_REF] Hubbell | The unified neutral theory of biodiversity and biogeography[END_REF]. It is a zero-sum game: the total number of individuals in a community is assumed constant, and dead individuals are immediately replaced by offspring of the remaining individuals with equal chance to reproduce for any of them. Neutral model has challenged former pattern-to-process approaches of ecological niches, which were based on permutational apporaches, by showing that neutrality itself could generate non-random structure and thus should be filtered out of patterns using specific models [START_REF] Bell | The co-distribution of species in relation to the neutral theory of community ecology[END_REF][START_REF] Canard | Emergence of structural patterns in neutral trophic networks[END_REF]. More precisely, it has been shown that this model yields a local distribution of species abundances within communities that follows a Dirichlet-multinomial distribution [START_REF] Donnelly | Likelihoods and Simulation Methods for a Class of Nonneutral Population Genetics Models[END_REF][START_REF] Etienne | A dispersal-limited sampling theory for species and alleles[END_REF][START_REF] Harris | Linking Statistical and Ecological Theory: Hubbell's Unified Neutral Theory of Biodiversity as a Hierarchical Dirichlet Process[END_REF]. In these studies, the Dirichlet-multinomial distribution has sometimes been called a 'dispersal-limited' multinomial distribution [START_REF] Etienne | A dispersal-limited sampling theory for species and alleles[END_REF].

The Dirichlet-multinomial distribution has several practical interests for pattern-to-process analyses of empirical communities. First, it satisfies useful property denoted weak closure under addition implying in particular that (i) if two species are lumped together, the multivariate distribution of the resulting distribution is still a Dirichlet-multinomial distribution with a natural adaptation of parameters [START_REF] Laroche | Analyzing snapshot diversity patterns with the neutral theory can show functional groups' effects on community assembly[END_REF]; (ii) considering a subgroup of species, the distribution of species abundances conditionnally to the size of the subgroup is also Dirichlet-multinomial, with a natural adaptation of parameters [START_REF] Laroche | A neutral theory for interpreting correlations between species and genetic diversity in communities[END_REF]. These properties have been used to design several tests in empirical ecological studies. Second, it is readily implemented as a hierarchical process [START_REF] Harris | Linking Statistical and Ecological Theory: Hubbell's Unified Neutral Theory of Biodiversity as a Hierarchical Dirichlet Process[END_REF], hence facilitating the computational aspects or the use of a Bayesian framework for inference purposes.

Importantly, the Dirichlet-multinomial distribution is robust to relaxing the unrealistic zero-sum assumption, by modelling the stochastic dynamics of species abundances within a community as a multivariate jump process. Community size then becomes a random variable fluctuating in time. Using this approach, [START_REF] Haegeman | Relaxing the zero-sum assumption in neutral biodiversity theory[END_REF] identified a class of neutral models for which decomposing species abundances stationary distribution within a community into a sum distribution (i.e. total number of individuals) and a conditional split distribution (i.e. species labels of individuals) led to a split distribution that is still a Dirichlet multinomial.

Although those results about Dirichlet-multinomial split have triggered an important and fruitful area of empirical research, they still represent a fraction of the possibilities offered by the neutral theory. Most of the litterature about detecting deviation from neutrality in species abundance data has focused on the split distribution conditionally to the total number of individuals, but to our knowledge only few studies simultaneously discuss the sum and the split (i.e. the full multivariate abundance distribution) as a tool to evaluate community composition [START_REF] Etienne | The zero-sum assumption in neutral biodiversity theory[END_REF].However this field has been mostly abandoned in favour to other statistical frameworks based for instance on the multivariate Poisson -lognormal distribution in the context of the joint species distribution model (JSDM) [START_REF] Aïtchison | The multivariate poisson-log normal distribution[END_REF][START_REF] Warton | So many variables: Joint modeling in community ecology[END_REF][START_REF] Ovaskainen | Joint Species Distribution Modelling: With Applications in R. Ecology, Biodiversity and Conservation[END_REF]. In addition, studies focusing on the split distribution have mostly focused on formulations of the neutrality assumption and ancillary hypotheses that necessarily lead to Dirichlet-multinomial distributions. If these assumptions are relaxed or modified, we expect that new split distributions can be obtained in models that still arguably remain neutral.

Our aim here is to show how relaxing some assumptions of neutral models can generalize the sum-split decomposition with Dirichlet-multinomial split to the more general Polyasplitting distributions family. We show that our framework covers classical distributions but also lead to new ones. Hence our work promotes a unified neutral-based statistical framework able to tackle the full multivariate abundance distribution of species, thus making an interesting mechanistic alternative to current phenomenological JSDMs framework. We also present a new general property denoted the strong closure under addition. This property allows in particular proposing approaches based on recursive application of splitting distributions to generate communities mixing dependent or independent species or group of species simultaneously Section 2 describes the family of multivariate Pólya splitting distributions. We specifically focus on nine examples of such distributions sharing the property of strong closure under addition. Section 3 shows that these distributions are stationary solutions of the master equation under a specific parametric hypothesis on the ratio between birth and death rates. Section 4 shows that this parametric assumption corresponds to a mildly extended version of the neutral theory of biodiversity where the immigration rate of a species can depend on its local abundance following a relationship identical across all species.

Polya splitting distributions

The first subsection recall the definition of multivariate Pólya distributions as urn models with n random drawing. Then, assuming that n is a random number, the second subsection presents the larger family of Pólya splitting distributions, introduced by Peyhardi and Fernique (2017) and more generally studied by [START_REF] Peyhardi | Splitting models for multivariate count data[END_REF]; [START_REF] Peyhardi | On quasi pólya thinning operator[END_REF]. Nine examples of Pólya splitting distributions with remarkable properties are presented, which are repeatedly referred to in the rest of the article.

Multivariate Pólya distributions

The Pólya urn model is generally presented in terms of n random drawings of balls from an urn, that initially contains θ j ∈ N * balls of the j th color. One ball is drawn at random and then replaced with c ∈ Z additional balls of the same color. A negative value for c means that balls are removed from the urn. This procedure is repeated n times and focus is made on the count N = (N 1 , . . . , N J ) of drawn balls for the J ≥ 2 different colors. Let |N | = J j=1 N j denotes the sum of the vector N and ∆ n = n ∈ N J : |n| = n (resp. n = n ∈ N J : |n| ≤ n ) the discrete simplex (resp. the discrete corner of the hypercube). The multivariate count distribution for N is known as the multivariate Pólya distribution and will be denoted by P

P |N |=n (N = n) = J j=1 n j -1 k=0 r [c] θ j (k j ) n-1 k=0 r [c] |θ| (k) , where r [c] θ (k) = θ+ck k+1 1 θ+ck≥0 , θ = (θ 1 , . . . , θ J ) ∈ Θ J c . The indicator function ensures that r [c] θ (k) ≥ 0 even if c < 0. Let us define R [c] θ (n) = n-1 k=0 r [c]
θ (k), then the pmf becomes

P |N |=n (N = n) = J j=1 R [c] θ j (n j ) R [c] |θ| (n)
.

The multivariate Pólya distribution turns out to be the multivariate hyergeometric distribution when c = -1, the multinomial distribution when c = 0 and the multivariate negative hyergeometric distribution when c = 1. Reasoning by equivalence on the pmf, it can be shown that these three distributions are the representative elements of their equivalence classes:

{P [c] ∆n (θ) : c < 0}, {P [c] ∆n (θ) : c = 0} and {P [c]
∆n (θ) : c > 0}. Therefore, in the following, a focus will be made only on the three cases c ∈ {-1, 0, 1}. Let us note that the pmf for the last two cases can be extended for continuous values of θ j ∈ R * + for j = 1, . . . , J. The multivariate negative hyergeometric distribution is thereby extended to the Dirichlet multinomial distribution. In summary, the three cases c ∈ {-1, 0, 1} respectively correspond to the

• multivariate hypergeometric distribution, denoted by H ∆n (θ) with Θ -1 = N * and R [-1] θ (n) = θ n ,
• multinomial distribution, denoted by M ∆n (π) where π = θ/|θ|,

Θ 0 = R * + and R [0] θ (n) = θ n ,
• Dirichlet-multinomial distribution, denoted by DM ∆n (θ) where Θ 1 = R * + and R

[1]

θ (n) = n+θ-1 n ,
The pmf of these three distributions are presented in Table A.2 of Appendix A. The support is ∆ n when c = 0 or c = 1 and is ∆ n ∩ θ when c = -1, i.e., the intersection between the simplex ∆ n and the hyper-rectangle θ = {n ∈ N J : n 1 ≤ θ 1 , . . . , n J ≤ θ J }. It should be noted that some authors refers to the Dirichlet multinomial distribution as the multivariate Pólya distribution. All along the paper, the multivariate Pólya distribution will refer to the general case that encompasses the three cases c ∈ {-1, 0, 1}.

Specific Pólya splitting distributions with remarkable properties

Pólya distributions cannot be considered as a sensu stricto multivariate distribution. Indeed, the sum of the random vector N is fixed and only J -1 elements over J are free. This kind of distribution, supported on ∆ n , is said to be singular. Another way to relax the fixed sum assumption, is considering the sum |N | as a random variable. The sum |N | then follows an univariate count distribution L(ψ) and the vector N = (N 1 , . . . , N J ) given the sum |N | = n follows a multivariate Pólya distribution P

[c] ∆n (θ 1 , . . . , θ J ). We thus obtain a Polya splitting distribution, that can be viewed as a compound distribution denoted as follows:

N ∼ P [c] ∆n (θ) ∧ n L(ψ),
where ψ is an univariate or multivariate set of unknown parameters (e.g. univariate for the Poisson distribution and bivariate for the negative binomial distribution). To be more explicit the pmf of the Pólya splitting distribution is given by P

(N = n) = P |N |=|n| (N = n)P (|N | = |n|).
According to L(ψ), several multivariate distributions can be defined sharing interesting properties.

The weak closure under addition The weak closure under addition is the stability under convolution of R

[c] θ : for all θ, θ ∈ Θ c R [c] θ * R [c] θ = R [c] θ+θ , (1)
where * denotes the discrete convolution. This identity plays a central role among the family of Polya splitting distributions and corresponds to the Vandermonde's identity, the Newton's binomial theorem and the Hagen-Rothe's identity when c = -1, 0 and 1 respectively. It could be shown that equation (1) implies the stability of the multivariate distribution under marginalization; see Theorem 1 of [START_REF] Peyhardi | Splitting models for multivariate count data[END_REF] for details. Let focus on two such stability properties.

(i) if two species are lumped together, the multivariate distribution of the resulting distribution is still a Pólya splitting distribution with a natural adaptation of parameters [START_REF] Laroche | Analyzing snapshot diversity patterns with the neutral theory can show functional groups' effects on community assembly[END_REF], i.e., we have

(N 1 + N 2 , N 3 , . . . , N J ) ∼ P [c] ∆n (θ 1 + θ 2 , θ 3 , . . . , θ J ) ∧ n L(ψ),
(ii) considering a subgroup of species, the distribution of species abundances conditionally to the size of the subgroup is also Pólya distribution with a natural adaptation of parameters [START_REF] Laroche | A neutral theory for interpreting correlations between species and genetic diversity in communities[END_REF], i.e., we have

(N 1 , N 2 , N 3 )|N 1 + N 2 + N 3 = n ∼ P [c] ∆n (θ 1 , θ 2 , θ 3 ).
It should be noted that all Pólya splitting distributions share this property, i.e., it holds for any sum distribution L(ψ).

The strong closure under addition A Pólya splitting distribution is said to be strongly closed under addition if the sum distribution and all the marginal distributions belong to the same family of parametric distributions, i.e., if we have ∀j ∈ {1, . . . , J} ∃ψ j : N j ∼ L(ψ j ).

While weak property is share by all Pólya splitting distributions, strong property holds for only specific sum distribution see [START_REF] Peyhardi | On quasi pólya thinning operator[END_REF] for details about this closure property.

In the following, we propose nine sum distributions verifying strong closure property. The first three distributions naturally extend singular Pólya distributions to their non-singular version. The three others leads to independence splitting distributions. Finally, the three last distributions allows the generalization of classical univariate count processes to their multivariate version.

Canonical cases

In the canonical case, the sum distribution is defined as the univariate version of the nonsingular Pólya distribution P

[c] m (θ, γ) . The sum distribution is then denoted by P

[c] m (θ, γ); see Table A.1 for details about its pmf and support. Theorem 4 of [START_REF] Peyhardi | Splitting models for multivariate count data[END_REF] showed that we have the following distribution identity

P [c] ∆n (θ) ∧ n P [c] m (|θ|, γ) = P [c] m (θ, γ).
See Table 1 to write this identity in the three cases c ∈ {-1, 0, 1}. It could be shown that marginals belongs to the same family of distribution, more precisely N j ∼ P

[c] m (θ j , |θ -j | + γ) (the strong closure under addition holds). The non-singular version has the advantage that its support n has a dimension equal to J (whereas the support of the singular version ∆ n has a dimension equal to J -1). Therefore the J variables N 1 , . . . , N J are free (not related by linear function) even if they are not independent. The graphical model of independence of such a distribution is complete [START_REF] Peyhardi | On quasi pólya thinning operator[END_REF]. In summary for the canonical case, the variables are free and not independent but the support is bounded.

Independence cases

There exists a sum distribution L * such that all the variables N 1 , . . . , N J are mutually independent. Moreover, all the margins belong to the same family L * , i.e. the strong closure under addition holds. It could be shown that the pmf of this distribution is

P (|N | = n) = R [c] |θ| (n)α n ∞ m=0 R [c] |θ| (m)α m , (2)
for some α > 0 and the marginals are given by

P (N j = n) = R [c] θ j (n)α n ∞ m=0 R [c] θ j (m)α m . ( 3
)
The distribution L * belongs to the family of power series distributions. It turns out to be the binomial distribution (c = -1), the Poisson distribution (c = 0) and the negative binomial distribution (c = 1) respectively; see Table 1 for details.

Dependent non-canonical cases

The last line of Table 1 presents three Pólya splitting distributions that share the strong closure under addition without independence assumption and different of the canonical case (demonstrations are given in Appendix B). Those distributions are obtained from the independent case assuming the parameter α (see eq. ( 2)) is a random variable. For instance when c = -1, the sum distribution is a binomial distribution compound by a beta distribution, i.e., a beta binomial distribution. According to (3), it is easily seen that marginals also follow beta binomial distributions. This particular distribution can be viewed as a new multivariate extension of the beta binomial distribution different from the more usual Dirichlet multinomial (non singular version). It is supported on the hyper-rectangle θ . For c = 0, the Poisson distribution is coumpound by a gamma distribution to obtain a negative binomial distribution. The corresponding Pólya splitting distribution turns out to be the natural multivariate extension, i.e., the negative multinomial distribution. Finally, for c = 1, the negative binomial is compound by a beta distribution to obtain a negative beta binomial distribution, also known as the univariate generalized Waring distribution [START_REF] Irwin | The generalized Waring distribution applied to accident theory[END_REF]. The corresponding Pólya splitting distribution turns out to be its natural multivariate extension, known as the multivariate generalized Waring distribution (MGWD) [START_REF] Xekalaki | The multivariate generalized Waring distribution[END_REF]. This dependent and non-canonical case has been recently formalized by [START_REF] Peyhardi | On quasi pólya thinning operator[END_REF] as the inverse Pólya distribution.

Stationary distributions of multivariate birth-death processes

We show that the class of Pólya splitting distributions introduced in former section exactly corresponds to the stationary distributions of multivariate birth-death processes under specific parametric assumptions on jumping rates. We further provide more precise conditions

c = -1 c = 0 c = 1 Canonical H ∆n (θ) ∧ n H m (|θ|, γ) M ∆n (π) ∧ n B m (p) DM ∆n (θ) ∧ n βB m (|θ|, γ) cases = = = H m (θ, γ) M m (p • π) DM m (θ, γ) θ ∈ N * J , γ ∈ N * , m ∈ N * , m ≤ |θ| + γ π ∈ ∆, m ∈ N * θ ∈ R * J + , γ ∈ R * + , m ∈ N * (support = ( m \ m-γ ) ∩ θ ) (support = m ) (support = m ) Independent H ∆n (θ) ∧ n B |θ| (p) M ∆n (π) ∧ n P(λ) DM ∆n (θ) ∧ n N B(|θ|, p) cases = = = J j=1 B θ j (p) J j=1 P(π j λ) J j=1 N B(θ j , p) θ ∈ N * J , p ∈ (0, 1) π ∈ ∆, λ ∈ R * + θ ∈ R * J + , p ∈ (0, 1) (support = θ ) (support = N J ) (support = N J ) Dependent H ∆n (θ) ∧ n βB |θ| (a, b) M ∆n (π) ∧ n N B(a, p) DM ∆n (θ) ∧ n βN B(|θ|, a, b) non-canonical = = = cases Mβ θ (a, b) N M(a, p • π) MGWD(b, θ, a) θ ∈ N * J , a ∈ R * + , b ∈ R * + π ∈ ∆, a ∈ R * + , b ∈ R * + θ ∈ R * J + , a ∈ R * + , b ∈ R * + (support = θ ) (support = N J ) (support = N J )
Table 1: Nine Pólya splitting distributions that share the strong closure under addition. The notation N ∼ J j=1 L j means that each N j follows the distribution L j independently of all other j = j. on rates that lead to the nine Pólya splitting distributions previously described (see Table 1). Let us start with the univariate case (J = 1) to recall classical results that will be useful to explicitly describe the stationary distribution of the sum in the multivariate case (J ≥ 2).

The univariate case

Let N (t) denote an univariate birth/death process with q + (n) (resp. q -(n)) denoting the birth rate (resp. death rate) for a population of size n. Let p n = P (N = n) denote the pmf of the stationary distribution. It can be shown, by induction on n, that solutions of the master equation at stationary state are solutions of the detailed balance equation:

q -(n + 1)p n+1 = q + (n)p n
Remark that for every count distribution (p n ) n∈N , it is possible to find two sequences of birth and death rates such that the detailed balance holds. Moreover, assuming that the death rate is positive, i.e., q -(n) > 0 for all n ≥ 1, then the detailed balance becomes

p n+1 = p n q(n),
where q(n) = q + (n)/q -(n + 1). The support of such a distribution is necessary of the connected form {0, . . . , m} with m ∈ N ∪ {∞} and we have:

p n = Q n k≥0 Q k ,
where Q n = n-1 k=0 q(k) for n ≥ 1 and Q 0 = 1. Specific stationary distribution are obtain according to specific parametric assumption on q(n). For instance if q(n) = 1 n+1 α for some α ∈ R * + then the stationary distribution is a Poisson distribution with parameter α. See Appendix C for several examples of parametric assumption on q(n) that lead to usual univariate distributions (e.g. binomial, negative binomial).

The multivariate case

Here we now describe the multivariate jump process N (t) = {N 1 (t), . . . , N J (t)} depicting species abundances within a community. We then focus on sufficient conditions on jumping rates to ensure the existence of a stationary distribution with detailed balance. In this specific case, it is straightforward to derive a closed form of the stationary distribution. In the following let p n (t) = P {N (t) = n} denote the pmf at time t and p n = P (N = n) denote the pmf at stationary state.

The master equation

∂p n (t) ∂t = J j=1
p n-e j (t)q - j (n -e j ) + p n+e j (t)q + j (n + e j ) -p n+e j (t){q - j (n) + q + j (n)} (4)

where q - j (n) (resp. q + j (n)) denotes the jumping rate from n to ne j (resp. to n + e j ) and e i is a vector of size J where the elements are all equal to zero except the i th equal to one. It is usual to assume that q - j (n) > 0 for all n ∈ N J such that n j > 0 (i.e., any individual is mortal). Moreover it is assumed that q + j (0) > 0 for all j = 1, . . . , J (where 0 = (0, . . . , 0) denotes the null vector) in order to avoid the case of non observed species (i.e., P (N j = 0) = 1).

The Detailed balance equation is given, for all j ∈ {1, . . . , J} and all n ∈ N J , by

p n+e j q - j (n + e j ) = p n q + j (n). ( 5
)
Since it is assumed that q - j (n) > 0 for all n ∈ N J such that n j > 0 then the detailed balance becomes p n+e j = p n q j (n)

where q j (n) = q + j (n)/q - j (n + e j ). The idea is then to recursively use the equality (6) in order to express p n according to p 0 and thus derive a closed analytical formula for the stationary distribution pmf of the multivariate jump process.

Necessary and sufficient conditions for detailed balance Note that the detailed balance holds if and only if the product of the quantitites q j along every increasing path between 0 and n is the same. This is equivalent to assume that this product is constant along every path between n and n + e i + e j with i = j. There are only two such path: n → n + e i → n + e i + e j and n n + e j n + e j + e i . For the first path the detailed balance equation gives p n+e i +e j = p n q i (n) q j (n + e i ). For the second path the detailed balance equation gives p n+e i +e j = p n q j (n) q i (n + e j ). Therefore a necessary and sufficient condition to the existence of a solution is the equality

q i (n) q j (n + e i ) = q j (n) q i (n + e j ), (7)
for all i = j ∈ {1, . . . , J} and all n ∈ N J such that p n = 0. Equation ( 7) corresponds to the Kolmogorov's criterion in the case of a multivariate birth-death process with positive death rates. This criterion is a necessary and sufficient condition for the reversibility of the process.

Parametric assumption on birth and death rates Assume that there exists some parameters c ∈ {-1, 0, 1} and θ = (θ 1 , . . . , θ J ) ∈ Θ J c and two non-negative functions s + and s -such that the birth and death rates have the following form

q + j (n) = s + (|n|)(θ j + cn j )1 θ j +cn j ≥0 , q - j (n) = s -(|n|)n j . (8)
The birth-deah rate thus becomes

q j (n) = s(|n|) r [c] θ j (n j ). ( 9
)
where s(n) = s + (n) s -(n+1) for all n ∈ N and r

[c] θ (n) = θ+cn n+1 1 θ+cn≥0 . It is easily seen that this parametric assumption (9) respects the Kolmogorov's criterion (7) and thus the detailed balance equation (5). In order to obtain a well defined stationary distribution we add the following assumption on s(n):

n≥0 R [c] |θ| (n) n-1 k=0 s(k) < ∞. ( 10
)
Theorem 1 Assume that the hypothesis (9) and (10) hold then

• the stationary distribution for N = (N 1 , . . . , N J ) is the Polya splitting distribution

P [c] ∆n (θ) ∧ n L
• L is the stationary distribution of a univariate process with birth/death ratio equal to

q(n) = s(n)r [c]
|θ| (n) (its support is necessary of the form {0, . . . , m} where m ∈ N∪{∞}).

Proof: Under assumption (9), using recursively (6), we obtain the pmf of the stationary distribution as follows:

p n = p 0 |n|-1 m=0 s(m) J j=1 n j -1 k=0 r [c] θ j (k).
Using notations of Section 2 we obtain

p n = p 0 |n|-1 m=0 s(m) J j=1 R [c] θ j (n j ).
Since the pmf is written as a recurrent product, it can be remarked that the support of the stationary distribution is exactly the connexe m when c = 0 or c = 1 and the connexe m ∩ θ when c = -1, where m is the smaller integer such that s(m) = 0. Note that if m = +∞ then the support becomes N J when c = 0 or c = 1 and θ when c = -1. Otherwise, We know that the probability of n can be conditioned by the sum as follows:

p n = P (|N | = |n|) P |N |=|n| (N = n).
By identifiability between the two previous equalities (on the support) we obtain the sum distribution and the split distribution (repartition into components given the sum). More precisely we have

P |N |=n (N = n) = J j=1 R [c] θ j (n j ) C(n) ,
for all n ∈ N, where C(n) is the normalizing constant

C(n) = n∈∆n J j=1 R [c] θ j (n j ),
and

P (|N | = n) = C(n) n-1 k=0 s(k) m≥0 C(m) m-1 k=0 s(k)
.

Note that C(n) is positive since the identification is made on the support. Now remark that C(n) turns out to be the convolution

C = R [c] θ 1 * • • • * R [c]
θ J . According to the Newton's binomial theorem, for c = 0, (respectively the Vandermonde's identity for c = -1 and the Hagen-Rothe identity for c = 1) we have

R [c] θ 1 * • • • * R [c] θ J = R [c]
|θ| , and thus the pmf of N given the sum |N | = n is

P |N |=n (N = n) = J j=1 R [c] θ j (n j ) R [c] |θ| (n) , i.e., given |N | = n we have N ∼ P [c]
∆n (θ) (multivariate Polya distribution). The sum distribution is now given by

P (|N | = n) = R [c] |θ| (n) n-1 k=0 s(k) m≥0 R [c] |θ| (m) m-1 k=0 s(k) . ( 11
)
This is a proper distribution according to assumption (10). Moreover, by definition we have R

[c] |θ| (n) = n-1 k=0 r [c] |θ| (k) (with convention R [c]
|θ| (0) = 1) and so the sum distribution can be viewed as the stationary distribution of an univariate process whose the ratio of birth/death is q(n) = s(n)r |θ| (n).

Parametric hypothesis on s(n)

For any univariate count distribution L with a support of the form {0, . . . , m} where m ∈ N ∪ {∞}, there exists a birth/death rate q(n) (and thus a function s(n)) such that the corresponding birth-death process converges toward the stationary distribution L. Therefore, for any Pólya splitting distribution P where m ∈ N ∪ {∞}) we are able to write the birth and death rates (through r θ (n) and s(n)) such that the multivariate jump process converges toward the stationary distribution P

[c] ∆n (θ) ∧ n L.
Let us illustrate this fact with nine parametric assumptions on s(n) that lead to the nine Pólya splitting distributions of Table 1. According to Theorem 1 the sum is a univariate birth-death process driven by the birth/death rate q(n) = s(n)r

[c] |θ| (n). Appendix C details the parametric form for q(n) that leads to specific univariate distributions (binomial, Poisson, etc . . . ). By identification it is possible to find the parametric form of s(n) in each case c ∈ {-1, 0, 1} and deduce the sum distribution as compiled in table 2. It is also possible to find the parametric form of s(n) in a general way, i.e., for any c ∈ {-1, 0, 1}; see Appendix D for details about the canonical cases.

Setting apart the canonical case with c = 1 and γ < 1, which will be discussed later, we observe clear qualitative differences among the nine examples. Each case corresponds to a particular variation profile of s(n). The function s decreases in the canonical case, is constant in the independent case and increases in the dependent non-canonical case (see Appendix F for details). Moreover, as c increases in {-1, 0, 1}, the convexity of s(n) shows consistent changes which do not depend on the considered case. When c = -1, s(n) is convex, when c = 0, s(n) is linear and, when c = 1, s(n) is concave.

The exception to above pattern is the canonical case with c = 1 when γ < 1. In this case, s(n) is increasing and convex, thus resembling the dependent non-canonical case with c = -1. This exception is obtained when the sum follows the beta binomial distribution βB m (|θ|, γ) with γ < 1 inducing a peak at n = m, hence promoting saturation of community size at m.

c = -1 c = 0 c = 1 Canonical s(n) = m -n γ -m + n + 1 1 m-γ≤n<m s(n) = m -n γ 1 n<m s(n) = m -n γ + m -n -1 1 n<m cases γ ∈ N * , m ∈ N * , m ≤ |θ| + γ γ ∈ R * + , m ∈ N * γ ∈ R * + , m ∈ N * Independent s(n) = α s(n) = α s(n) = α cases α ∈ R * + α ∈ R * + α ∈ (0, 1) Dependent s(n) = a + n |θ| + b -n -1 1 n<|θ| s(n) = a + n |θ| + b s(n) = b + n |θ| + a + b + n non-canonical cases a ∈ R * + and b ∈ R * + a ∈ R * + and b ∈ R * + a ∈ R * + and b ∈ R * +
Table 2: Parametric hypothesis on s(n) that lead to the nine Pólya splitting distributions of Table 1 4 Biological interpretations of the Pólya splitting distributions

From a biological perspective, the birth and death rates of the multivariate jump process defined by the master equation (see eq. ( 4)) are generally assumed to have the following form:

q + j (n) = m j (n) + n j b j (n) q - j (n) = n j d j (n) where m j (n)
is the immigration rate of species j, b j (n) is the per-capita local reproduction rate of species j and d j (n) is the per-capita local death-or-emigration rate of species j. Classical neutral models make two additional assumptions. The first corresponds to the strong neutrality:

q + j (n) = m(|n|)π j + n j b(|n|) q - j (n) = n j d(|n|) (12)
The second assumption is detailed balance. Given the strong neutrality assumption and discarding degenerated cases m(0) = 0 or m(1) = 0, detailed balance occurs if and only if there exists a constant Ĩ ≥ 0 such that

b(n) = Ĩm(n) (13)
for all n ∈ N * (see Appendix E). The parameter I = 1/ Ĩ his known as the effective number of migrants [START_REF] Etienne | A novel genealogical approach to neutral biodiversity theory[END_REF]). Here we focus on a generalization of equations ( 12), considering new expressions of birth and death rates as follows:

q + j (n) = m(|n|) π j -Kn j + n j b(|n|) 1 π j -Kn j ≥0 q - j (n) = n j d(|n|) (14)
with K ∈ R + . The master equation with assumption ( 14) can be called a neutral model with density-dependent immigration, because the rate of immigration of species j is now m(|n|)(π j -Kn j ), which depends on species j and on the local population size n j . Assuming K = 0 leads to the usual equations ( 12). We still assume (13), which is a necessary and sufficient condition to obtain detailed balance (see Appendix E.2 for details about the necessity). As corollary of Theorem 1, we obtain the link between the neutral model with density-dependent immigration ( 14) and assumption (8) that leads to Pólya splitting distributions.

Corollary 1 Any Polya splitting distribution (with sum distribution support of the form {0, . . . , m} where m ∈ N ∪ {∞}) can be obtained as a stationnary distribution of a jumping process with rates (14) verifying detailed balance condition (13). Reciprocally, let a jumping process with rates (14) with detail balance condition (13). If at least K or Ĩ is null then the stationary distribution is a Polya splitting distribution.

Proof: Let P According to Theorem 1 this is the stationary distribution of a jumping process satisfying assumption (9) on q j (n). Let

π j = θ j /|θ| and m(n), d(n) such that m(n) d(n + 1) = |θ|s(n). (15)
Moreover, let assume that

K = |θ| -1 if c = -1, 0 if c = 0 or c = 1, and Ĩ = 0 if c = -1 or c = 0, |θ| -1 if c = 1. Then q j (n) = m(|n|) d(|n| + 1) π j + ( Ĩ -K)n j n j + 1 1 θ j +cn j ≥0 , q j (n) = m(|n|) π j -Kn j + b(|n|)n j d(|n| + 1)(n j + 1) 1 θ j +cn j ≥0 ,
where b(n) = Ĩm(n) and 1 θ j +cn j ≥0 = 1 π j -Kn j ≥0 and thus the desired result. Reciprocally, let a jumping process with rates ( 14) with detail balance condition ( 13). If at least K or Ĩ is null then the parametrization is reversible.

In the previous section, we derived results about the variation of s(n) among the nine Pólya splitting distributions presented in table 1. Migration m(n) and death d(n) rates directly relate to s(n) through s(n) = |θ| -1 m(n) d(n+1) . Therefore, when assuming that one of these rates is constant with community size n, one obtains biologically interpretable results about the density-dependence effect on the other rate. For instance, if one assumes that d(n) = 1 (i.e., constant per-capita local death or emigration rate), m(n) = |θ|s(n). Here again discarding the canonical case with c = 1 and γ < 1, we obtain that m(n) always decreases in canonical cases, is constant in independent cases and increases in dependent noncanonical cases. From a biological perspective, this suggests that, as community size increases and before hitting a potential regulation threshold m < ∞, immigration and reproduction of within the community become harder in the canonical case (negative density dependence), remain unaffected in the independent case and become easier in the other case (positive density dependence). We also obtain that m(n) is convex when c = -1, which suggest that the marginal variation in reproduction and immigration increases as community density increases. Migration m(n) is linear when c = 0 and is concave when c = 1, which means that the marginal variation in reproduction and immigration increases as community density decreases (see Table 3).

Discussion and perspectives

We presented the Polya-splitting distributions, a set of multivariate distributions with two key properties : the sum is a positive random variable on N and the split conditionally to the sum is a Polya distribution. We recalled that Polya-distribution (the split) can be classified

c = -1 c = 0 c = 1 All K = |θ| -1 K = 0 K = 0 cases Ĩ = 0 Ĩ = 0 Ĩ = |θ| -1 d(|n|) = 1 d(|n|) = 1 d(|n|) = 1 Canonical m(|n|) = |θ| (m-|n|) γ+1-(m-|n|) 1 |n|≤m m(|n|) = |θ| γ (m -|n|)1 |n|≤m m(|n|) = |θ| m-|n| γ-1+m-|n| 1 |n|≤m cases b(|n|) = 0 b(|n|) = 0 b(|n|) = m-|n| γ-1+m-|n| 1 |n|≤m Independent m(|n|) = α|θ| m(|n|) = α|θ| m(|n|) = α|θ| cases b(|n|) = 0 b(|n|) = 0 b(|n|) = α Dependent m(|n|) = |θ| a+|n| |θ|+b-|n|-1 m(|n|) = |θ| a+|n| |θ|+b m(|n|) = |θ| b+|n| |θ|+a+b+|n| non-canonical cases b(|n|) = 0 b(|n|) = 0 b(|n|) = b+|n| |θ|+a+b+|n|
Table 3: Parametric hypothesis on extended neutral models leading to the nine Pólya splitting distributions presented in table 1. Dark gray indicates cases where the split distribution has been characterized [START_REF] Haegeman | Relaxing the zero-sum assumption in neutral biodiversity theory[END_REF] while light gray stands for cases where both split and sum distributions have been characterized [START_REF] Etienne | The zero-sum assumption in neutral biodiversity theory[END_REF]. Other cases (white boxes) are new.

in three categories, depending on a parameter c = -1, 0 or 1. The case c=1 corresponds to a Dirichlet-multinomial split, the case c=0 to a multinomial split and the case c=-1 to a hypergeometric split. Our main contribution is to connect those distributions to the neutral theory of biodiversity in ecology, a useful null model allowing the evaluation of non-neutral processes such as adaptation or natural selection [START_REF] Alonso | The merits of neutral theory[END_REF]. We found that for any Polya-splitting distribution, irrespective of the value of c in the split but with a sum support of the form {0, . . . , m} where m ∈ N ∪ {∞}, there exists a multivariate jump process of neutral species with such stationnary distribution. However, staying at the very general level for the sum distribution, the associated transition rates may not have a straightforward biological interpretation. We therefore exhibited, nine transition rates parametrization with meaningful biological interpretation leading to usual parametric distributions. Reciprocally, if multivariate jump process of neutral species follows the detailed balance assumption and has a well-defined stationary distribution of the sum, then the multivariate distribution of species counts is a Polya-splitting distribution. [START_REF] Kadmon | Integrating the effects of area, isolation, and habitat heterogeneity on species diversity: A unification of island biogeography and niche theory[END_REF] had already shown that a variety of neutral jump processes with detailed balance assumption (formalized in eq. 2 in their work) could generate Multinomial or Dirichlet-Multinomial split distributions that correspond to Polya distributions with c = 0 or 1 respectively. However, they did not identify neutral jump processes that could generate c = -1 because they made the classic assumption that the positive jump rate of a species j linearly increased with the local population size of species j with a positive or null slope b that corresponds to a per-capita bitrh rate. Under this assumption, only Multinomial and Dirichlet-multinomial split can be obtained (c = 0 or 1). Because we started from the description of the full family of Polya-distribution, we were able to ask the question of whether Polya-splitting distribution with hypergeometric split (c = -1) could also be obtained as stationnary distributions of neutral jump processes -which we showed to be true -and what were the peculiarities of these processes that depart from e.g. models condidered in [START_REF] Kadmon | Integrating the effects of area, isolation, and habitat heterogeneity on species diversity: A unification of island biogeography and niche theory[END_REF]. We evidenced that the key difference is that hypergeometric split when the positive jump rate of species j linearly decreases with its local population size. This would correspond to a negative per capita birth rate in [START_REF] Kadmon | Integrating the effects of area, isolation, and habitat heterogeneity on species diversity: A unification of island biogeography and niche theory[END_REF], hence explaining why they did not explore this track. Here, we showed that this negative relationship between positive jump rate and local population size can emerge in a well-defined case: when a limited quantity of available propagules can immigrate from the regional pool to a local site where local birth rate is zero, i.e. the community is a pure sink.

We focused our study on the stationary state of a biological community following a neutral jump process. In practice, only a fraction of the community is observed, through a sampling process. In the case c = 0 or 1, [START_REF] Etienne | A dispersal-limited sampling theory for species and alleles[END_REF] noticed that the Dirichletmultinomial split of the sum among species verify a "subsampling property": when applying a hypergeometric sampling process with fixed sample size over the community, the resulting subsample still followed a Dirichlet-multinomial distribution with the same parameters. We conjecture that this property still holds for the case c = -1, although this remains to be properly shown. This property remarkably simplifies the statistical study of species relative abundances within the community. By contrast, the hypergeometric sampling with fixed size does not allow studying the sum distribution, because the sample size is artificially controlled, independently from the real community size. This is quite limiting: our results emphasized the importance of studying the total sum abundance as a random variable, because of its links with the dependence structures among species. In particular, we show that independence is a consequence of parametric assumption made on birth and death rates and not a necessary assumption per se, contrary to what was posited by other authors [START_REF] Etienne | The zero-sum assumption in neutral biodiversity theory[END_REF]. Therefore we suggest that a stronger focus should now be given on sampling processes that preserves information about the community size, like process that controls the distance covered or the time spent during sampling rather than the number of individuals. This requires to explicit sampling models accounting for the spatio-temporal distribution of studied organisms [START_REF] Jousimo | A spatio-temporally explicit random encounter model for large-scale population surveys[END_REF] and to study stability properties of associated thinning operators [START_REF] Peyhardi | On quasi pólya thinning operator[END_REF].

Pólya splitting distributions induce only two types of dependence structures: either all species are independent or fully dependent with homogenized correlation sign. To extend this binary setting towards more complex nested dependence structures between species or communities, we suggest the use of recursive application of splitting distributions. From this perspective, the strong closure under addition property plays an essential role by preserving distributions across levels, hence allowing a full control of generated dependencies. This emphasizes the importance of using appropriate choices of Pólya splitting distributions at each level to ensure strong closure under addition.

For instance, let us assume that we aim at simulating species communities composed of five species s 1 , . . . , s 5 . The first three species (s 1 , s 2 and s 3 ) belongs to a community C 1 and are all fully positively correlated, corresponding to mutualist species context or species sharing underlying environmental factors. The two others species (s 4 and s 5 ) belongs to a second community, C 2 correlated to C 1 . These two species s 4 and s 5 are assumed to be independent. Such situation can be easily obtained combining negative binomial for the sum distributions and multinomial or Dirichlet multinomial distributions for split components (with specific constraint on parameters); see figure 1). Another example relies on the simulation of negatively correlated species within the first community C 1 (exclusive species context) and non-dependent species within C 2 and assuming independence between C 1 and C 2 . Such simulation can be performed combining binomial distributions for sum parts with multinomial and/or hypergeometric distribution for splits. However, does such sequential structure is solution of the master equation? It remains an open question and should be carefully study.

While multinomial or Dirichlet multinomial regression has already been used in neutral community analysis [START_REF] Jabot | Reconciling neutral community models and environmental filtering: theory and an empirical test[END_REF][START_REF] Jabot | A stochastic dispersal-limited trait-based model of community dynamics[END_REF], the inclusion of environmental factors in the Pólya splitting distributions is a natural extension. It could be achieved assuming parameters varied according to covariates as follows:

P [c] ∆n {θ(x)} ∧ n L{ψ(x)},
where x = (x 1 , . . . , x p ) denotes the vector explanatory variables (see [START_REF] Peyhardi | Splitting models for multivariate count data[END_REF] for more details in the multinomial splitting regression context). Note that parameters could be constrained to be the same in the split and sum parts.

Finally, combining graph hierarchical approach with the inclusion of environmental covariates at each node leads to propose nested multi-level inhomogenous splitting models. Such models should be interesting alternatives to classical approaches used in joint species distribution contexts mainly based on conditional Independence's [START_REF] Warton | So many variables: Joint modeling in community ecology[END_REF][START_REF] Ovaskainen | Joint Species Distribution Modelling: With Applications in R. Ecology, Biodiversity and Conservation[END_REF] and the use of the multivariate Poisson log-Normal distribution. Comparatively, our approach allows to model dependencies between species at the observation level, while the classical JSDM's model dependencies at the latent process strata. While correlation relationships estimated at the latent processes inform correlations between observations, it does not allow to deduce dependencies structures at the observation scales [START_REF] Aïtchison | The multivariate poisson-log normal distribution[END_REF][START_REF] Chiquet | The poisson-lognormal model as a versatile framework for the joint analysis of species abundances[END_REF]. A null correlation does not imply independence in the multivariate Poisson context.

N B(θ

1 + θ 2 , p) M ∆n (π 1 , π 2 ) N B (θ 1 + θ 2 , p 1 ) M ∆n (π 1 , π 2 , π 3 ) N B(θ 1 + θ 2 , p 1 ) N B(θ 1 + θ 2 , p 2 ) N B(θ 1 + θ 2 , p 3 ) N B(θ 1 + θ 2 , p 2 ) DM ∆n (θ 1 , θ 2 ) N B(θ 1 , p 2 ) N B(θ 2 , p 2 ) Dependent communities
Dependent species Independent species

[c] m (θ, γ) θ ∈ Θ c , γ ∈ Θ c , m ∈ N * R [c] θ (n)R [c] γ (m -n) R [c] θ+γ (m) hypergeometric (c = -1) H m (k, l) k ∈ N * , l ∈ N * ,m ∈ N * , m ≤ k + l {0, . . . , m} ∩ {m -l, . . . , k} k n l m-n k+l m binomial (c = 0) B m (p) p ∈ (0, 1), m ∈ N * {0, . . . , m} m n p n (1 -p) m-n beta binomial (c = 1) βB m (a, b) a ∈ R * + , b ∈ R * + ,m ∈ N * {0, . . . , m} n+a-1 n m-n+b-1 m-n m+a+b-1 m
[c] ∆n (θ) n ∈ N * , θ ∈ Θ J c J j=1 R [c] θ j (n j ) R [c] |θ| (n) multivariate hypergeometric (c = -1) H ∆n (k) n ∈ N * , k ∈ N * J ∆ n ∩ k J j=1 k j n j |k| n multinomial (c = 0) M ∆n (π) n ∈ N * , π ∈ ∆ ∆ n n n J j=1 π n j j Dirichlet multinomial (c = 1) DM ∆n (α) n ∈ N * , α ∈ R * J + ∆ n J j=1 n j +α j -1 n j n+|α|-1 n Non-singular version multivariate Pólya P [c] m (θ, γ) m ∈ N * , θ ∈ Θ J c , γ ∈ Θ c R [c] γ (m -|n|) J j=1 R [c] θ j (n j ) R [c] |θ|+γ (m) multivariate hypergeometric (c = -1) H m (k, l) m ∈ N * , k ∈ N * J , l ∈ N * ( m \ m-l ) ∩ k l m-|n| J j=1 k j n j |k|+l m m ≤ |k| + l multinomial (c = 0) M m (π * ) m ∈ N * , π * ∈ m m n (1 -|π * |) m-|n| J j=1 π n j j Dirichlet multinomial (c = 1) DM m (α, β) m ∈ N * , α ∈ R * J + , β ∈ R * + m m-|n|+β-1 m-|n| J j=1 n j +α j -1 n j m+|α|+β-1 m
Table A.2: Notations and pmf of multivariate Pólya distributions (singular and non-singular versions)

• c = 1 According to Theorem 1 of [START_REF] Peyhardi | Splitting models for multivariate count data[END_REF] the marginals are given by the beta-binomial damage distribution

βB n (θ j , |θ -j |) ∧ n βN B(|θ|, a, b) = βB n (θ j , |θ -j |) ∧ n N B(|θ|, p) ∧ p β(a, b) = βB n (θ j , |θ -j |) ∧ n N B(|θ|, p) ∧ p β(a, b) = N B(θ j , p) ∧ p β(a, b) βB n (θ j , |θ -j |) ∧ n βN B(|θ|, a, b) = βN B(θ j , a, b)
This demonstration follows the same ways as in case c = -1.

C Some stationary distributions of univariate birthdeath process

According to different assumptions on the ratio q

(n) = q + (n) q -(n+1) , we find different distributions (p n ) n≥0 . Recall that p n = Q n m≥0 Q m , with Q n = n-1 k=0 q(k).
C.1 Univariate Pólya distributions

• Hypergeometric distribution: if q(n) = (k-n)(m-n) (n+1)(l-m+n+1) 1 max(0,m-l)≤n<min(m,k) with k ∈ N * , l ∈ N * and m ≤ k + l then Q n = ( k n )( l m-n ) ( l m ) and p n = k n l m-n k+l m , max(0, m -l) ≤ n ≤ min(m, k), i.e., N ∼ H m (k, l). • Binomial distribution: if q(n) = (M -n)π (n+1)(1-π) 1 n<m with π ∈ (0, 1) and m ∈ N * then Q n = m n π n (1 -π) n and p n = m n π n (1 -π) m-n , 0 ≤ n ≤ m i.e., N ∼ B m (π).

D Parametric hypothesis on s(n) for the canonical case

Assume that s

(n) = 1/r [c]
γ (m -n -1) for some γ ∈ Θ and m ∈ N. Then we have

n-1 k=0 s(k) = 1 n-1 k=0 r [c] γ (m -k -1) = 1 r [c] γ (m -1) × • • • × r [c] γ (m -n) = r [c] γ (m -n -1) × • • • × r [c] γ (0) r [c] γ (m -1) × • • • × r [c] γ (0) n-1 k=0 s(k) = R [c] γ (m -n) R [c] γ (m) Therefore m n=0 R [c] |θ| (n) n-1 k=0 s(k) = 1 R [c] γ (m) m n=0 R [c] |θ| (n)R [c] γ (m -n) = 1 R [c] γ (m) (R [c] |θ| * R [c] γ )(m) m n=0 R [c] |θ| (n) n-1 k=0 s(k) = 1 R [c] γ (m) R [c] |θ|+γ (m)
Finally the pmf of the sum given by (11) becomes

P (|N | = n) = R [c] |θ| (n)R [c] γ (m -n) R [c] |θ|+γ (m)
, for all n ≤ m and zero otherwise. This is the pmf of the univariate Pólya distribution P Recall the assumptions of classic neutral theory regarding the jumping rates (12):

q + j (n) = m(|n|)π j + n j b(|n|) q - j (n) = n j d(|n|) where we assume that n > 0 =⇒ d(n) > 0 (i.e. no individual is immortal).

We assume that the above process has a stationnary distribution with support S ⊂ N J . Because no individual is immortal, 0 ∈ S and if n ∈ S and n ∈ N J | ∀i ∈ {1, ..., J}, n i ≤ n i then n ∈ S.

Define We seek for necessary conditions to obtain detailed balance of the stationary distribution which is depicted by the Kolmogrorov criterion (7): ∀n ∈ K , ∀(i, j) ∈ {1, ..., J} 2 , q i (n) q j (n + e i ) = q j (n) q i (n + e j ), where q j (n) = q + j (n) q - j (n+e j ) . Recall that q j is well defined because no individual is immortal. Using the expression of jumping rates, the Kolmogorov criterion becomes :

  ). Its probability mass function (pmf) is given by

  a sum support of the form {0, . . . , m}

L

  be a Pólya splitting distribution with sum distribution support of the form {0, . . . , m} where m ∈ N ∪ {∞}.

Figure 1 :

 1 Figure 1: Example of possible simulation schemes combining multi dependent levels

  , γ). Therefore the multivariate stationary distribution is the non-singular version of the Pólya distribution thanks to the identity P According to Theorem 1 of[START_REF] Peyhardi | Splitting models for multivariate count data[END_REF] the marginals are given by the damage distributionsP [c] n (θ j , |θ -j |) ∧ n P [c] m (|θ|, γ) = P [c] m (θ j , |θ| + γ),for all j = 1, . . . , J.E Necessary and sufficient condition for detailed balance in the classic neutral theory and with densitydependent migration E.1 Classic neutral theory

  K m = min{n ∈ N|m(n) = 0} and K = min{n ∈ N|m(n) = b(n) = 0}. By definition, K m ≤ K ≤ +∞. If K m = 0 then S = {0}. If K m > 0 then S = K .In what follows we assume that K m > 0 and S = K .

  [m(|n|)π i + n i b(|n|)] [m(|n| + 1)π j + n j b(|n| + 1)] = [m(|n|)π j + n j b(|n|)] [m(|n| + 1)π i + n i b(|n| + 1)] ,which can be simplified as :(π i n j -π j n i ) [m(|n|)b(|n| + 1) -m(|n| + 1)b(|n|)] = 0 which implies in turn that : ∀n ∈ {1, ..., K -1}, m(n)b(n + 1) = m(n + 1)b(n) (16)Declining constraint (16) along possible initializations of m(n) and b(n) yields:• if m(1) = b(1) = 0 : K = 1 and constraint (16) disappears. • if m(1) = 0 and b(1) > 0 : ∀n ∈ {1, ..., K -1}, m(n) = 0 • if m(1) > 0 and b(1) = 0 : ∀n ∈ {1, ..., K -1}, b(n) = 0 • if m(1)> 0 and b(1) > 0 : ∀n ∈ {1, ..., K -1} : m(n) > 0, b(n) > 0 and : • if m(1) = b(1) = 0 : K = 1 and constraint (17) disappears. • if m(1) = 0 and b(1) > 0 : ∀n ∈ {1, ..., K -1}, m(n) = 0 • if m(1) > 0 and b(1) = 0 : ∀n ∈ {1, ..., K -1}, b(n) = 0 • if m(1) > 0 and b(1) > 0 : ∀n ∈ {1, ..., K -1} : m(n) > 0, b(n) > 0 and : b(n) = Ĩm(n)where Ĩ = b(1) m(1) . In summary, we have shown that the stationnary distribution verifies detailed balance only if one of the following conditions holds• K m = 0 • K m = 1 and ∀n ∈ {1, ..., K -1}, m(n) = 0 • K m > 1 and ∃ Ĩ ≥ 0 | ∀n ∈ {1, ..., K }, b(n) = Ĩm(n)Reciprocally, it is straightforward to show that each of these conditions is sufficient to obtain detailed balance of the stationary distribution.

F

  Variation and convexity ofn → s(n) m -n -1) 2 1 n≤m cases γ ∈ N * , m ∈ N * , m ≤ |θ| + γ γ ∈ R * + , m ∈ N * γ ∈ R * + , m ∈ N * =⇒ s decreases on (m -γ, m) =⇒ s decreases on (0, m) γ < 1 =⇒ s increases on (0, m) γ > 1 =⇒ s decreases on (0, m)=⇒ s increases on (m -γ, m) =⇒ s constant on (0, m) s increases on (0, m -1) 1 2 < γ < 1 =⇒ s increases on (0, m) γ > 1 =⇒ s decreases on (0, on (0, |θ|) =⇒ s constant on (0, +∞) =⇒ s increases on (0, +∞) =⇒ s increases on (0, |θ|) =⇒ s constant on (0, +∞) =⇒ s decreases on (0, +∞)

  Table A.1: Notations and pmf of univariate Pólya distributions

	Name	Notation	space parameter	support	pmf (p n )
	Singular version
	multivariate Pólya	P

Table F .

 F 1: Variation and convexity of s for the nine Pólya splitting distributions of Table1

Acknowledgements

This research was supported by the GAMBAS project funded by the French National Research Agency (ANR-18-CE02-0025). F. Laroche was supported by the BloBiForM project funded by French National Research Agency (ANR-19-CE32-0002-01).

Contribution

F. Laroche: Conceptualization, Formal analysis, Writing -Original Draft F. Mortier: Conceptualization, Writing -Original Draft, Funding acquisition. J Peyhardi: Methodology, Formal analysis, Writing -Original Draft.

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Declarations

For the purpose of Open Access, a CC-BY 4.0 public copyright licence (https://creativecommons. org/licenses/by/4.0/) has been applied by the author to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission.

Name

Notation

space parameter support pmf (p n ) Let us show the strong closure under addition for the three Pólya splitting distributions presented in the third line of Table 1. We have to show that marginal distributions and sum distribution belong to the same family.

• c = -1 According to Theorem 1 of [START_REF] Peyhardi | Splitting models for multivariate count data[END_REF] the marginals are given by the hypergeometric damage distribution

The first equality uses the definition of the beta-binomial distribution, the second one uses the Fubini theorem (inversion of sum on n and integral on p), the third one uses the stability of the binomial distribution under hypergeometric damage process (obtained in the case of independence) and the last one uses again the definition of the beta-binomial distribution.

• c = 0 According to Theorem 1 of [START_REF] Peyhardi | Splitting models for multivariate count data[END_REF] the marginals are given by the binomial damage distribution B n (π j ) ∧ n N B(r, p). Theorem 6 of [START_REF] Peyhardi | Splitting models for multivariate count data[END_REF] showed the stability of the negative binomial distribution under the binomial damage process, i.e., we have

where p = π j p π j p+1-p . The demonstration is based on the generative function of a binomial damage distribution. This result can also be obtained by following the way of the demonstartion of the previous case, recalling that a negative binomial is a Poisson mixed by a gamma distribution.

and

C.2 Other univariate distributions

i.e., N ∼ P(α).

• Negative binomial distribution: if q(n) = n+α n+1 π with π ∈ (0, 1) and α ∈ (0, ∞) then

i.e., N ∼ N B(α, π).

Geometric distribution: if q(n) = π with π ∈ (0, 1) then Q n = π n and

Remark that the geometric distribution is a sub-case of the negative-binomial distribution, more precisely we have G(π) = N B(1, π). α+a,b) and

where B(a, b) = Γ(a)Γ(b) Γ(a+b) , i.e., N ∼ βN B(α, a, b). Remark that if the parameters a and b are positive integers, then the beta-binomial distribution turns out to be the negative hypergeometric distribution. Otherwise, the beta negative binomial distribution is also called the generalized waring distribution [START_REF] Irwin | The generalized Waring distribution applied to accident theory[END_REF].

m(1) . In summary, we have shown that the stationnary distribution verifies detailed balance only if one of the following conditions holds

Reciprocally, it is straightforward to show that each of these conditions is sufficient to obtain detailed balance of the stationary distribution.

E.2 Neutral theory with density-dependent migration

We now turn to the extension of neutral theory including density-dependent immigration as defined in equation ( 14), which we recall here :

with K ∈ R + . We define K j = max{n j ∈ N|π j -Kn j > 0} + 1 and K the hypercube {0, ..., K 1 } × {0, ..., K 2 } × ... × {0, ..., K J }.

We also define K m and K like in previous section and we assume that K m > 0. Then the support of the stationary distribution is S = K ∩ K

We seek for necessary conditions to obtain detailed balance of the stationary distribution which is depicted by the Kolmogrorov criterion (7). Using the expression of jumping rates in ( 14