A reversible allelic partition process and Pitman sampling formula

TL;DR

This paper introduces a new continuous-time Markov chain model for dynamic allelic partitions that extends existing sampling formulas, maintains mutual dependence over time, and exhibits reversible behavior with a stationary distribution.

Contribution

It presents a novel Markov chain model extending the Pitman sampling formula, capturing reversible dynamics and asymptotic behavior in allelic partition processes.

Findings

The model maintains mutual dependence among multiplicities over time.

When death rate exceeds birth rate, the system has a reversible distribution.

The population converges to a stationary configuration with finite families and individuals.

Abstract

We introduce a continuous-time Markov chain describing dynamic allelic partitions which extends the branching process construction of the Pitman sampling formula in Pitman (2006) and the birth-and-death process with immigration studied in Karlin and McGregor (1967), in turn related to the celebrated Ewens sampling formula. A biological basis for the scheme is provided in terms of a population of individuals grouped into families, that evolves according to a sequence of births, deaths and immigrations. We investigate the asymptotic behaviour of the chain and show that, as opposed to the birth-and-death process with immigration, this construction maintains in the temporal limit the mutual dependence among the multiplicities. When the death rate exceeds the birth rate, the system is shown to have reversible distribution identified as a mixture of Pitman sampling formulae, with negative binomial mixing distribution on the population size. The population therefore converges to a stationary random configuration, characterised by a finite number of families and individuals.

Equations

$$\mathbb{A}_{n}=\bigg{\{}\mathbb{m}\in\mathbb{Z}^{n}_{+},\ \sum_{i=1}^{n}im_{i}=n\bigg{\}}$$

$$s(m):=\sum _{i\ge 1}i{m}_i,k(m):=\sum _{i\ge 1}{m}_i,$$

$$\mathbb{A}:=\bigcup _{n\ge 1}{\mathbb{A}}_n=\{m\in {\mathbb{Z}}_{+}^{\infty },s(m)<\infty \},$$

$$\begin{array}{r}{\text{ESF}}_n^{\theta }(m)=\frac{n!}{{\theta }_{(n)}}\prod _{i\ge 1}{\left(\frac{\theta }{i}\right)}^{m_i}\frac{1}{m_i!}\text{mathbbm}{1}_{\{s(m)=n\}},\end{array}$$

$$\begin{array}{r}{\text{PSF}}_n^{\alpha ,\theta }(m)=\frac{n!}{{\theta }_{(n)}}\prod _{i=1}^{k(m)}(\theta +i\alpha )\prod _{j\ge 1}{\left[\frac{(1-\alpha {)}^{j-1}}{j!}\right]}^{m_j}\frac{1}{m_j!}\text{mathbbm}{1}_{\{s(m)=n\}}.\end{array}$$

$$\begin{array}{r}p(m^{\prime }|m)\propto \{\begin{array}{ll}\theta +\alpha k(m),& m^{\prime }=m+e_1,\\ (i-\alpha )m_i,& m^{\prime }=m-e_i+e_{i+1},i\ge 1,\\ 0,& \text{else},\end{array}\end{array}$$

$$\begin{array}{r}\underset{n\to \infty }{\lim }\frac{K_n^{0,\theta }}{\log n}=\theta \text{a.s.},\end{array}$$

$$\begin{array}{r}\underset{n\to \infty }{\lim }\frac{K_n^{\alpha ,\theta }}{\log n}=S^{\alpha ,\theta }\text{a.s.},\end{array}$$

$$\underset{h\to 0}{\lim }\frac{1}{h}\mathbb{P}(I(t+h)-I(t)=1|I(t)=k)=\theta +\alpha k,k\ge 0,$$

$$\begin{array}{r}{q}_{i,i+1}=i-\alpha ,i\ge 1.\end{array}$$

$$\begin{array}{r}\mathbb{M}(t)\sim \prod _{i\ge 1}\text{Po}(\theta b_t^i/i),b_t=\{\begin{array}{ll}[e^{(1-\mu )t}-1]/[e^{(1-\mu )t}-\mu ],& \mu \ne 1\\ t/(1+t),& \mu =1.\end{array}\end{array}$$

$$\begin{array}{r}S(t)\sim \text{N-Bin}(\theta ,b_t),\mathbb{P}(S(t)=n)=\frac{{\theta }_{(n)}}{n!}(1-b_t{)}^{\theta }{b}_t^n,n=0,1,\dots \end{array}$$

$$\begin{array}{r}\mathbb{P}(\mathbb{M}(t)=m|S(t)=n)={\text{ESF}}_n^{\theta }(m),m\in \mathbb{A},n\ge 1.\end{array}$$

$$\begin{array}{r}\mathbb{M}(t)\stackrel{d}{\to }(X_1,X_2,\dots ),t\to \infty ,\end{array}$$

$$\begin{array}{r}\pi (m)=C(\theta /\alpha {)}_{(k(m))}\prod _{i\ge 1}\mathrm{Po}(m_i;{\alpha }_i{\mu }^{-i}),{\alpha }_i=\frac{\alpha (1-\alpha {)}_{(i-1)}}{i!},m\in \mathbb{A},\end{array}$$

$$\begin{array}{r}C=e^{1-(1-1/\mu {)}^{\alpha }}(1-1/\mu {)}^{\theta }.\end{array}$$

$$\begin{array}{rl}q(m^{\prime }|m)& =\{\begin{array}{ll}\theta +\alpha k(m),& m^{\prime }=m+e_1,\\ (i-\alpha )m_i,& m^{\prime }=m-e_i+e_{i+1},i\ge 1,\\ \mu i{m}_i,& m^{\prime }=m-e_i+e_{i-1},i\ge 1,\\ -\theta -(1+\mu )s(m),& m^{\prime }=m,\\ 0,& \mathrm{else}.\end{array}\end{array}$$

$$\begin{array}{rl}S(t)& :=s(\mathbb{M}(t))=\sum _{i\ge 1}i{M}_i(t)\end{array}$$

$$\begin{array}{r}r(n^{\prime }|n)=\{\begin{array}{ll}\theta +n,& n^{\prime }=n+1,\\ \mu n,& n^{\prime }=n-1,\\ -\theta -(1+\mu )n,& n^{\prime }=n,\\ 0,& \text{otherwise}.\end{array}\end{array}$$

$$\theta +\alpha k(m)+\sum _{i\ge 1}(i-\alpha )m_i=\theta +n.$$

$$\sum _{i\ge 1}\mu i{m}_i=\mu n.$$

$$q_{i,i+1}=i-\alpha ,q_{i,i-1}=\mu i,i\ge 1$$

$$\underset{t\to \infty }{\lim }\mathbb{P}(S(t)=n|S(0)=0)=\lambda (n),\forall n\in {\mathbb{Z}}_{+}$$

$$\sum _{n\in {\mathbb{Z}}_{+}}\lambda (n)r(n^{\prime }|n)=0,\forall n^{\prime }\in {\mathbb{Z}}_{+}.$$

$$\begin{array}{r}\lambda (n)=\frac{{\theta }_{(n)}}{n!}{\mu }^{-n}(1-1/\mu {)}^{\theta },n\ge 0,\end{array}$$

$$\begin{array}{r}\lambda (n)r(n^{\prime }|n)=\lambda (n^{\prime })r(n|n^{\prime }),n\ne n^{\prime }.\end{array}$$

$$\begin{array}{rl}\lambda (n)r(n+1|n)& =\frac{{\theta }_{(n)}}{n!}{\mu }^{-n}(1-1/\mu {)}^{\theta }(\theta +n)\\ & =\frac{{\theta }_{(n+1)}}{(n+1)!}{\mu }^{-(n+1)}(1-1/\mu {)}^{\theta }\mu (n+1)\\ & =\lambda (n+1)r(n|n+1);\end{array}$$

$$\begin{array}{rl}\lambda (n)r(n-1|n)& =\frac{{\theta }_{(n)}}{n!}{\mu }^{-n}(1-1/\mu {)}^{\theta }\mu n\\ & =\frac{{\theta }_{(n-1)}}{(n-1)!}{\mu }^{-(n-1)}(1-1/\mu {)}^{\theta }(\theta +n)\\ & =\lambda (n-1)r(n|n-1).\end{array}$$

$$q(m^{(1)}|m)q(m^{(2)}|m^{(1)})\dots q(m^{\prime }|m^{(l)})>0.$$

$$S(t)=\sum _{i\ge 1}i{M}_i(t)=0⟺\mathbb{M}(t)=e_0=(0,0,\dots ),t\ge 0,$$

Full text

Keywords: birth-and-death process, branching process, immigration, stationary process, Pólya urn, population dynamics.

MSC Primary: 60G10, 60J10. Secondary: 92D25, 60J80.

1 Allelic partitions and sampling formulae

A partition of $n\in\mathbb{N}$ is an unordered collection $\pi=\{n_{1},\dots,n_{k}\}$ of $k\leq n$ positive integers whose sum equals $n$. A common equivalent way to describe $\pi$ is by means of the so-called allelic partition, which groups the partition sets by size. We denote by $\mathbb{m}=(m_{1},\dots,m_{n})\in\mathbb{Z}^{n}_{+}$ the associated vector of multiplicities, where $m_{i}$ counts the number of repetitions of $i$ in $\pi$, and let

[TABLE]

be the finite set of all allelic partitions of $n$. It will be useful to embed $\mathbb{A}_{n}$ into $\mathbb{Z}^{\infty}_{+}$ by considering infinite vectors $\mathbb{m}=(m_{1},\dots,m_{n},0,0,\dots)$, defining also

[TABLE]

which are respectively the number of items and groups (or positive multiplicities) in $\mathbb{m}$, and

[TABLE]

the countable set of all allelic partitions. A random allelic partition $\mathbb{M}=(M_{i})_{i\geq 1}$ can then be defined as a random variable taking values in $\mathbb{A}$.

The literature on the distributional properties of random allelic partition is rich and well established. Ewens ​(1972) derived the distribution of the random allelic partition induced by a sample of $n$ genes drawn from a selectively neutral population at equilibrium (see also Karlin & McGregor, 1972). This is described by the celebrated Ewens sampling formula (ESF), which assigns probability

[TABLE]

to the configuration $\mathbb{m}\in\mathbb{A}$ with $m_{i}$ alleles appearing exactly $i$ times for each $i\geq 1$. Here $\theta>0$ and $\theta_{(n)}=\theta(\theta+1)\dots(\theta+n-1)$ is the ascending factorial. The impact of the ESF has been significant in a number of fields beyond population genetics. For example, Antoniak ​(1974) derived it independently from a sample from a Dirichlet process, a cornerstone nonparametric prior distribution in Bayesian nonparametric statistics; Hoppe ​(1984) recovered the formula as the marginal distribution of a Markov chain generated by a Pólya-like urn model comprising an infinite number of colors, analogous to that introduced in Blackwell & MacQueen ​(1973) and described equivalently in Aldous ​(1985) with the famous metaphor of the Chinese restaurant process. See Johnson et al. ​(1997) and Crane ​(2016) for reviews and Feng ​(2010) for connections with population genetics.

Pitman ​(1995) introduced a two-parameter generalization of (1), often referred to as the Pitman sampling formula (PSF), whereby for parameters $0\leq\alpha<1$ and $\theta>-\alpha$, the probability assigned to a random allelic partition $\mathbb{m}\in\mathbb{A}$ is

[TABLE]

For $\alpha=0$, equation (2) immediately recovers (1), so that $\textnormal{PSF}_{n}^{0,\theta}\equiv\textnormal{ESF}_{n}^{\theta}$. For $0<\alpha<1$, the distribution arises for example in the study of stable processes with index $\alpha$, the case $\alpha=1/2$ being related to the zeros of Brownian motion (see also Pitman, 1997), and in the study of the partition structures induced by a sample from a two-parameter Poisson–Dirichlet distribution ​(Perman et al., 1992; Perman, 1993; Pitman & Yor, 1997). The corresponding generalization of Hoppe’s urn is defined in Pitman ​(1995, 1996b) as a Markov chain of allelic partitions with initial state $\mathbb{e}_{0}$ and transition probabilities

[TABLE]

where $\mathbb{e}_{i}=(\delta_{ij})_{j\geq 1},\ i\geq 0$, and the normalising constant is given by $\theta+s(\mathbb{m})$. The first transition in (3) can be understood as the addition to the gene sample of a new gene of a previously unobserved allele, while the transition from $\mathbb{m}$ to $\mathbb{m}-\mathbb{e}_{i}+\mathbb{e}_{i+1}$ corresponds to adding a gene whose allele had been previously observed $i$ times. For $\alpha=0$, (3) coincides with Hoppe’s urn, while for any $0\leq\alpha<1$ the marginal distribution of the state of the chain after $n$ transitions is given by (2).

Among other relevant contributions to the related theory, the asymptotic behavior of the number $K^{\alpha,\theta}_{n}$ of different alleles was derived by Korwar & Hollander ​(1973) in the case $\alpha=0$, for which

[TABLE]

and by Pitman ​(1997) for $0<\alpha<1$, where we have that

[TABLE]

$S^{\alpha,\theta}$ being an absolutely continuous random variable on $(0,\infty)$ whose law is related to the Mittag–Leffler distribution. Concerning the PSF, we further refer the reader to Kerov ​(2006); Pitman ​(1996a) for various characterizations thereof, James et al. ​(2008) for connections with Bayesian nonparametric statistics and again Feng ​(2010) for connections with population genetics. Note that similar asymptotic results to the above are available for Gibbs-type models ​(Gnedin & Pitman, 2006; De Blasi et al., 2015), which include the one- and two-parameter models recalled above, as well as a model in Gnedin ​(2010) which produces a finite but random number of families in the limit.

In this paper we are interested in the connection between Ewens–Pitman sampling formulae and the theory of stochastic processes for population growth. In particular, we construct a birth-and-death process with immigration whose dynamics modifies the one arising as after (3) by allowing the removal of items from the system; and we show, in a particular regime, that the process is reversible with respect to a certain mixture of Pitman sampling formulae. The rest of the paper is organised as follows. In Section 2 we present our contribution and its connections with past related work. Section 3 describes in detail the birth-and-death population model with immigration, and how it can be constructed via branching processes. Finally, Section 4 analyses the reversible regime, identifying explicitly the reversible distribution and its connections with Pitman sampling formulae.

2 Main contribution and related work

The seminal work of Karlin & McGregor ​(1967) describes a population wherein new families are initiated at a sequence of random times generated by a stochastic process $I=\{I(t),t\geq 0\}$, which can be thought of as immigration events, and then evolve, independently of one another, according to the law of a common continuous-time Markov chain on $\mathbb{Z}_{+}$, whose infinitesimal rates are denoted $q_{ij}$. As pointed out by Tavaré ​(1989), if new families are interpreted as novel mutant alleles of a given gene, then the scheme may be regarded as a version of the infinitely-many-neutral-alleles model and can provide, under specific choices of its probabilistic components, a generating mechanism for the sampling formulae (1) and (2). In particular, let $I$ be a pure-birth processes started at [math] with birth rates

[TABLE]

and assume that the process describing the evolution of each family is started at 0 and has rates

[TABLE]

If now $M_{i}(t)$ is the number of families of size $i$ at time $t$, and $\mathbb{M}(t)=(M_{i}(t))_{i\geq 1}$ is the induced allelic partition, then $\mathbb{M}=\{\mathbb{M}(t),t\geq 0\}$ defines a continuous-time Markov chain on $\mathbb{A}$ whose embedded jump chain has transition probabilities (3) and marginal distributions described by the PSF; see Feng & Hoppe ​(1998) and Pitman ​(2006). When $\alpha=0$, this construction reduces to the birth process with immigration considered in Tavaré ​(1987), which provides the corresponding embedding of Hoppe’s urn leading to the ESF (cf. Section 2.7.2 in Feng, 2010).

In the above schemes, each transition entails the addition of one individual to the population, and both the population size and the number of families diverge almost surely as time increases. This is a convenient mathematical simplification but may be undesirable or result in a lack of flexibility in different contexts, and is arguably unrealistic when modelling the evolution of populations. Therefore, it can be of interest to account for the death of individuals and to study the related implications on the family structure.

Motivated by the above considerations, we consider here a modified version of the population model underlying the PSF, by assuming that the individuals are endowed with independent random exponential lifetimes of parameter $\mu>0$, with $\mu=0$ recovering the original construction, thereby introducing uniform death events in the scheme. See Definition 1 below for the details. By leaving unchanged the other rules governing the evolution of the family structure, we thus obtain more flexible dynamics whereby families are started, fluctuate in size and possibly become extinct with the passage of time, resulting in a population with varying size and number of families. Note that, because of exchangeability, the PSF is consistent with respect to uniform deletion, hence a death event induces a random partition still following the PSF. For consistency with respect to deletion of an entire family, a property known in the literature as regeneration, see Gnedin & Pitman ​(2005).

For $\alpha=0$, our construction recovers the detailed version of the birth-and-death process with immigration (BDI) in Karlin & McGregor ​(1967), object of thorough study in Kendall ​(1975) and Tavaré ​(1989). In this case, the new families are started at the times of a Poisson point process of intensity $\theta$, and each family subsequently evolves according to an independent linear birth-and-death process. Theorem 2.1 in Karlin & McGregor ​(1967) derives the marginal distribution of the induced allelic partition $\mathbb{M}(t)$, which is that of a sequence of independent Poisson random variables given by

[TABLE]

Anticipated by earlier results in Watterson ​(1974) and Kendall ​(1975), Tavaré ​(1989) reformulated the above distribution as a mixture of Ewens sampling formulae in (1) with a negative binomial mixing distribution on $n$. Specifically, letting $S(t)$ being the number of alive individuals at time $t$, the population size process $S=\{S(t),t\geq 0\}$ defines by construction a BDI process, with negative binomial marginal distribution

[TABLE]

Then, by combining (7) and (8), the ESF is recovered as the conditional distribution of the actual allelic partition given the population size, i.e.,

[TABLE]

Within the study of the BDI process, a great attention has been dedicated to investigating the long-run behavior of the model. By taking the limit as $t\to\infty$ in (7), we have the convergence in law

[TABLE]

where, for $0<\mu\leq 1$, $X_{i}\overset{\textnormal{ind}}{\sim}\textnormal{Po}(\theta/i)$, while for $\mu>1$, $X_{i}\overset{\textnormal{ind}}{\sim}\textnormal{Po}(\theta\mu^{-i}/i)$. In the first case the death rate does not exceed the birth rate, and the asymptotic regime corresponds to the weak limit as $n\to\infty$ of the Ewens partition structure (1) derived in Arratia et al. ​(1992), characterized by an infinite sample size and a logarithmic growth of the underlying number of groups. Instead, for $\mu>1$, the death events occur at a faster rate than births, causing each family to eventually become extinct almost surely (see Kendall, 1975) and preventing the indefinite growth observed previously. In fact, in the second regime the population size process $S$ is easily seen in Lemma 2 below to have $\textnormal{N-Bin}(\theta,\mu^{-1})$ reversible distribution. In turn, Kendall ​(1975) showed the reversibility of the induced allelic partition process, and shed light on the representation of the reversible distribution as an analogous mixture as that arising from (9), characterized by a $\textnormal{N-Bin}(\theta,\mu^{-1})$ mixing on the population size.

In this paper, we encode the description of the evolving family structure directly in terms of the induced allelic partition process, and then focus on investigating the long-run behavior of the model. Because of the more involved rules governing the formation of new families, the sharp distributional result derived in Karlin & McGregor ​(1967) is not accessible in our construction. Nonetheless, by observing that the overall dynamics of the total population size is left unchanged, we show that for $\mu>1$ the model is reversible, and identify the reversible distribution as a mixture of Pitman sampling formulae with respect to the same negative binomial mixing on the population size appearing for $\alpha=0$. For $0<\alpha<1$, the mixture can be written in closed form as

[TABLE]

where

[TABLE]

Thus, in the reversible regime, the system will converge to a stationary random configuration that features almost surely finitely many families and individuals. However, as opposed to the previous model, the mutual dependence among the multiplicities is seen in (11) to be preserved in the limit.

3 A birth-and-death process with immigration

3.1 Definition

Definition 1**.**

For $0\leq\alpha<1$, $\theta>-\alpha$ and $\mu>0$, let $\mathbb{M}=\{\mathbb{M}(t),\ t\geq 0\},\ \mathbb{M}(t)=(M_{i}(t))_{i\geq 1}$, be a continuous-time Markov chain with state space $\mathbb{A}$, initial state $\mathbb{M}(0)=\mathbb{e}_{0}$ and temporally homogeneous transition rates

[TABLE]

We notice that the first two lines in (13) coincide precisely with the transitions in the sequential construction by Pitman ​(1995) displayed in (3). These describe the initiation of a new family of size one and the addition of one individual to a pre-existing family of size $i$ respectively. The third transition instead represents, for $i\geq 2$, the death of a member of a family of size $i$, while for $i=1$ it corresponds to the extinction of a family comprising a single individual. The choice of the initial state $\mathbb{M}(0)=(0,0,\dots)$ is rather common across the literature, and can be interpreted as the population being initially empty, as in a territory about to be colonised.

The infinitesimal generator corresponding to the transition rates in (13) are easily seen to be stable, as $q(\mathbb{m}|\mathbb{m})$ is finite for all $\mathbb{m}\in\mathbb{A}$, and conservative, since each row has zero sum (cf. Norris, 1997). Furthermore, by adapting the argument in Section 4 in Kendall ​(1975), we can deduce from the properties of the associated population size process that $\mathbb{M}$ is also regular, i.e. it performs almost surely only a finite number of jumps in every finite time interval. In particular, for any $t\geq 0$, let

[TABLE]

be the random number of items underlying the allelic partition $\mathbb{M}(t)$. Then we have from Definition 1 that $S=\{S(t),t\geq 0\}$ defines a continuous-time Markov chain on $\mathbb{Z}_{+}$, with initial state $S(0)=0$ and transition rates

[TABLE]

Indeed, by properly aggregating the rates in (13), we see that conditionally given the population size $S(t)=n$ (that is to say $\mathbb{M}(t)=\mathbb{m}$ for any $\mathbb{m}\in\mathbb{A}$ with $s(\mathbb{m})=n$) one of the following events occur. A new individual is introduced in the population (either initiating a new family or joining a pre-existing one) increasing its size by one unit, and this happens with overall rate

[TABLE]

Alternatively, the death of one of the $n$ individuals occurs, decreasing the population size by one unit, which happens at the total rate

[TABLE]

Thus, irrespective of the parameter $\alpha$, $S$ defines a classical BDI process, with unit birth rate, death rate $\mu$ and immigration rate $\theta$. The law of $S$ is seen in (15) to be insensitive to the choice of $\alpha$, so that we conclude that the more refined rules governing the formation of new families do not impact the dynamics of the total population size, but rather influence the detailed probabilistic fluctuations in the family configuration.

Finally, we notice that $\mathbb{M}$ and $S$ share the same random jump time, as each transition of the allelic partition process involves either the introduction or the death of an individual. Therefore, as the generator of the birth-and-death process with immigration is known to be regular (see, e.g., Section 2 in Kendall ​(1975)), we can directly conclude about the regularity of $\mathbb{M}$, which accordingly identifies (jointly with the initial state) a well-defined chain on $\mathbb{A}$.

3.2 Branching process construction

A biological basis for the three-parameter process in Definition 1 can be provided in terms of a population of individuals grouped into families evolving according to a sequence of births, deaths and immigrations, by introducing uniform death events in the models studied by Feng & Hoppe ​(1998) and Pitman ​(2006).

First assume that $\theta>0$, and consider a population whose members have independent exponential lifetimes with parameter $\mu$, and autonomously reproduce according to independent unit-rate Poisson processes (taken without loss of generality, upon appropriate time rescaling). At a birth, any offspring produced by the oldest member of each family will either join its parent’s family, with probability $1-\alpha$, or start a new family, with probability $\alpha$. All other family members give birth to individuals that will remain within their parent’s family. Finally, an independent Poisson process of rate $\theta$ is assumed to bring immigrants into the population that, at the time of entry, will initiate a new family each.

For $-\alpha<\theta\leq 0$, the above rules can be modified by assuming that no immigrant can enter the population, but rather that the overall oldest member in the population produces at rate $1-\alpha$ offspring that will join its family, and at rate $\alpha+\theta$ individuals that will initiate new families.

Thus, by exploiting basic properties of independent Poisson processes (e.g., a Poissonization argument as in Athreya & Karlin, 1968), it follows that, if we encode the family structure into the corresponding sequence of multiplicities, the dynamics of the above population generates a continuous-time Markov chain of allelic partitions, whose transition rates recover precisely those of the process $\mathbb{M}$ in (13). In particular, each family evolves according to an independent birth-and-death process with initial state 1 and transition rates

[TABLE]

directly generalizing (6). On the other hand, the introduction of deaths changes substantially the probabilistic laws that govern the input process $I=\{I(t),\ t\geq 0\}$ (counting the number of new families formed up to time $t$). In particular, as opposed to Karlin & McGregor ​(1967), here $I$ does not define a Markov process, since the probability of introducing a new family at each transition depends on the actual number of families $k(\mathbf{M}(t))$ which, due to families possibly becoming extinct, is in general different from $I(t)$. This in turn represents a serious obstruction in deriving distributional results for quantities of interest related to the model, in particular for the total number of families. Still, as resulting from (15), the dynamics of the total population size is tractable, and this fact will be exploited in the next section to study a reversibility regime for the process.

4 The reversible regime

As seen above, the population size process $S$ defines a classical BDI process. Hence, from the transition rates in (15) we immediately deduce that $S$ is irreducible, as each pair of states are mutually accessible. Moreover, the process is known to be transient for $0<\mu<1$, null recurrent for $\mu=1$ and positive recurrent for $\mu>1$ (cf. Section 2 in Kendall, 1975). In the latter case, $S$ is ergodic, implying the convergence in law

[TABLE]

for a unique stationary distribution $\lambda=(\lambda(n))_{n\geq 0}$ on $\mathbb{Z}_{+}$ satisfying

[TABLE]

The above sequence of equations can be solved recursively to obtain

[TABLE]

so that $\lambda$ results to be a negative binomial distribution $\textrm{N-bin}(\theta,\mu^{-1})$ with parameters $\theta$ and $\mu^{-1}$. The following lemma shows that, under the condition $\mu>1$, the population size process in fact enjoys the stronger property of being reversible with respect to $\lambda$.

Lemma 2**.**

Let $S$ be the population size process associated to $\mathbb{M}$ in Definition 1 through (14). Then, if $\mu>1$, $S$ is reversible with respect to its stationary distribution $\lambda$ in (16).

Proof.

As $S$ is regular and irreducible, the reversibility of the process is equivalent to the detailed balance condition

[TABLE]

In view of (15), we only need to show (17) for the cases $n\geq 0,\ n^{\prime}=n+1$, and $n\geq 1,\ n^{\prime}=n-1$, as the transition rates corresponding to all other transitions are equal to 0. For $n\geq 0$, $n^{\prime}=n+1$, the left hand side of (17) becomes

[TABLE]

and similarly, for $n\geq 1$, $n^{\prime}=n-1$ we have

[TABLE]

∎

We now move on to study the properties of $\mathbb{M}$. We first notice that $\mathbb{M}$ is irreducible, as we see from (13) that the state space $\mathbb{A}$ constitutes a single communication class, i.e., for any $\mathbb{m},\mathbb{m}^{\prime}\in\mathbb{A}$ we can find a finite path $\mathbb{m}^{(1)},\dots,\mathbb{m}^{(l)}\in\mathbb{A}$ for which

[TABLE]

Thus, we can argue as in Section 5 in Kendall ​(1975) to draw some preliminary conclusions concerning the influence of the death rate on the dynamics of the model, based upon the properties of the sample size process $S$. In particular, in view of the one-to-one correspondence

[TABLE]

we deduce that the initial condition $\mathbb{e}_{0}$ and, by irreducibility, the entire state space $\mathbb{A}$ itself, are characterized by the same recurrence and transience properties of the state [math] for $S$. Accordingly, we conclude that $\mathbb{M}$ is transient for $0<\mu<1$, null recurrent for $\mu=1$ and ergodic for $\mu>1$.

We shall henceforth focus on the latter case, wherein the underlying population almost surely undergoes a sequence of extinctions, separated by intervals of times of finite expected length during which its size and family configuration fluctuate according to the description in Section 3. Such dynamics were shown in Lemma 2 to generate a reversible population size process. The following theorem, which states our main result, shows that the same property holds for $\mathbb{M}$.

Theorem 3**.**

Let $\mathbb{M}$ be as in Definition 1, with $0<\alpha<1$ and $\mu>1$. Then, $\mathbb{M}$ is reversible with respect to the distribution $\pi$ with weights given by (11).

Proof.

That $\pi$ is a well-defined probability distribution on $\mathbb{A}$, as well as the explicit expression of the constant $C$ in (12), will follow from the mixture representation below. We then proceed to show the detailed balance condition

[TABLE]

which, since $\mathbb{M}$ is regular and irreducible, is equivalent to the reversibility.

In view of the expression for the transition rates (13) of $\mathbb{M}$, we only need to show (18) for the cases in which, for given $\mathbb{m}\in\mathbb{A}$, we have $\mathbb{m}^{\prime}=\mathbb{m}+\mathbb{e}_{1}$, $\mathbb{m}^{\prime}=\mathbb{m}-\mathbb{e}_{i}+\mathbb{e}_{i+1}$ for $i\geq 1$, or finally $\mathbb{m}^{\prime}=\mathbb{m}+\mathbb{e}_{i}-\mathbb{e}_{i+1}$ for some $i\geq 0$. The rate associated to the other transitions are indeed null.

For $\mathbb{m}\in\mathbb{A}$ and $\mathbb{m}^{\prime}=\mathbb{m}+\mathbb{e}_{1}$, denoting $k=k(\mathbb{m}),\ s=s(\mathbb{m}),$ so that $k(\mathbb{m}+\mathbb{e}_{1})=k+1$ and $s(\mathbb{m}+\mathbb{e}_{1})=s+1$, the left hand side of (18) becomes

[TABLE]

Here we have used the identity

[TABLE]

and the fact that

[TABLE]

with the convention $(1-\alpha)_{(0)}=1$, yielding $\alpha_{1}=\alpha$.

Instead, for $\mathbb{m}^{\prime}=\mathbb{m}-\mathbb{e}_{i}+\mathbb{e}_{i+1}$, denoting as before $k=k(\mathbb{m})$ and $s=s(\mathbb{m}),$ we have $k(\mathbb{m}-\mathbb{e}_{i}+\mathbb{e}_{i+1})=k$ and $s(\mathbb{m}-\mathbb{e}_{i}+\mathbb{e}_{i+1})=s+1$. Thus, proceeding similarly

[TABLE]

where we have used

[TABLE]

Finally, consider the case $\mathbb{m}^{\prime}=\mathbb{m}+\mathbb{e}_{i}-\mathbb{e}_{i+1}$. Using the same notation as above, we have $k(\mathbb{m}+\mathbb{e}_{i}-\mathbb{e}_{i+1})=k$ and $s(\mathbb{m}+\mathbb{e}_{i}-\mathbb{e}_{i+1})=s-1$. A similar computation yields

[TABLE]

concluding the proof. ∎

The reversible distribution $\pi$ in of Theorem 3 can be represented as a mixture of Pitman sampling formulae in (2) with respect to the reversible $\textrm{N-Bin}(\theta,\mu^{-1})$ distribution of the population size process $S$. Specifically, we have

[TABLE]

Indeed, for $\alpha>0$, the PSF (2) can be rewritten as

[TABLE]

where $\alpha_{i}$ is as in (11). Then, for $\mathbb{m}\in\mathbb{A}$, it follows that

[TABLE]

The series appearing as the argument of the exponential can be simplified by noticing that

[TABLE]

Replaced into the last line of the previous display, this recovers the expression for the reversible distribution $\pi$ with the constant $C$ as in (12), implying (19).

Theorem 3 yields interesting insights in the long-run behavior of the three-parameter process: in particular, in view of the ergodicity, $\mathbb{M}(t)$ will converge in distribution to a stationary random allelic partition whose distribution admits the mixture representation (19), and is characterized by a random, but almost surely finite number of items and groups. Thus, the reversible regime of $\mathbb{M}$ reveals clear analogies with the results in Kendall ​(1975) for the BDI process, essentially due to the equivalence of the induced population size processes. However, as opposed to the case $\alpha=0$, $\mathbb{M}$ in the limit to maintains the mutual dependence among the multiplicities, when $\alpha>0$.

Concerning the range $0<\mu\leq 1$, we deduce at once from the absence of positive recurrence that the process is necessarily non-stationary, and in fact, because of the slower death rate, the associated population size process $S$ is seen in (8) to diverge as $t\to\infty$. On the other hand, for $\theta>0$, we can deduce from the description of the population model in Section 3, that the number of families $K(t)$ alive at time $t$ stochastically dominates the number of families in the corresponding model with $\alpha=0$, as it can be rewritten as

[TABLE]

where $K_{1}(t)$ counts the family founded by immigrants, and $K_{2}(t)$ is the number of those started by the newborns leaving their family. Now $K_{1}(t)$ is evidently equal in distribution to the number of families comprised in a population governed by the BDI process, which, in turn, exhibits the logarithmic growth (4). Accordingly, we conclude that $\mathbb{M}$ will evolve towards infinite structures characterized by an infinite number of individuals and families. Whether a precise description of the limiting behavior similar to (10) can be achieved within this formulation, remains for the moment an open question.

Acknowledgements

The authors are grateful to the Associate Editor and two referees for helpful comments. The first author is supported by the ERC grant No. 647812 (UQMSI) and partially by the EPSRC grant EP/L016516/1 for the Cambridge Centre for Analysis. The second and third authors are partially supported by the Italian Ministry of Education, University and Research (MIUR), through PRIN 2015SNS29B and through “Dipartimenti di Eccellenza” grant 2018-2022.