- Research
- Open access
- Published:

# Theoretical size distribution of fossil taxa: analysis of a null model

*Theoretical Biology and Medical Modelling*
**volumeÂ 4**, ArticleÂ number:Â 12 (2007)

## Abstract

### Background

This article deals with the theoretical size distribution (of number of sub-taxa) of a fossil taxon arising from a simple null model of macroevolution.

### Model

New species arise through speciations occurring independently and at random at a fixed probability rate, while extinctions either occur independently and at random (background extinctions) or cataclysmically. In addition new genera are assumed to arise through speciations of a very radical nature, again assumed to occur independently and at random at a fixed probability rate.

### Conclusion

The size distributions of the pioneering genus (following a cataclysm) and of derived genera are determined. Also the distribution of the number of genera is considered along with a comparison of the probability of a monospecific genus with that of a monogeneric family.

## Background

Mathematical modelling of the evolution of lineages goes back at least to Yule[1] who developed the eponymous *Yule process* (homogeneous pure birth process) in which speciations occur independently and at random. Yule's model did not include extinctions *per se*, because he believed that they resulted only from cataclysmic events. This issue was discussed at greater length by Raup[2], who distinguished between background and episodic extinctions. Raup started from a homomogeneous birth-and-death process model (in which background extinctions occur, like speciations, independently and at random) for which he presented mathematical results, and described more complex models of extinction including episodic extinctions and a mixture of episodic and background extinctions. However he gave no mathematical results for these models. Stoyan[3] considered a time in-homogeneous birth-and death process, in which speciation and background extinction rates varied with time, based on the idea that younger paraclades have higher speciation rates, while older ones have higher background extinction rates.

There has been considerable discussion (*e.g*. Raup[2]; Patzkowsky[4]; Przeworski and Wall[5]) about the suitability of the null birth-and-death process model (with constant birth and death rates) as a macroevolutionary model of species diversification. In order to truly assess the validity of such a model it is necessary to have a full understanding of its properties which can then be compared with the fossil record. Specifically analysis is needed to generate hypotheses, which can be tested against available data. To date such an analysis is incomplete, relying on the partial analytic results of Raup[2] and the simulation results of Patzkowsky[4] and Przeworski and Wall[5].

Analytic results are clearly superior to simulation ones. In particular with analytic results for the size distribution of a clade one can fit the model via a multinomial likelihood, using observed size distributions, and thence test the adequacy of the underlying birth-and-death model using a statistical goodness-of-fit test. In addition analytic results are preferable to simulation ones, in that it is much easier to interpret a parametric formula than a collection of simulation results; and one does not have to distinguish between sampling variation due to a finite number of runs (noise) and signal.

It is the purpose of this paper to conduct a more thorough analysis of the birth-and-death model than that previosly carried out by Raup[2]. In particular we obtain results for size distributions of taxa and probabilities of monotypic taxa. In this paper we confine attention to obtaining analytic results and defer actual fitting and testing of the fit, using observed fossil data, to a future paper.

We develop the mathematical model presented by Raup[2] (and used in simulations by the above authors) to include the possibility of episodic, cataclysmic extinctions in which complete lineages are destroyed. We consider a hiearchy of models, which can include both cataclysmic and background extinctions of species and examine the resulting size distributions of extinct genera. We start (following section), as did Yule, by considering cataclysmic extinction only. Furthermore like Patzkowsky[4] and Przeworski and Wall [5], we assume that at any time an existing species can split, yielding a new species so radically different from existing ones that it becomes the founding member of a new genus. Thus we assume that the probability of a new genus being formed in an infinitesimal interval (*t*, *t* + *dt*) is proportional to the total number of species in existence at time *t*. We derive results for the size distribution of extinct genera.

In the third and fourth sections we do the same assuming only background extinctions (but no cataclysmic extinction); and both cataclysmic and background extinctions (although the results here are limited). The fifth section is devoted to the distribution of the number of genera derived from the pioneering species and in the final section the probability of a monotypic genus is compared with that of a monogeneric family.

## Cataclysmic extinctions only

Yule[1] considered the evolution of a genus begining with one species at time *t* = 0, which thenceforth evolves as a homogeneous pure birth process (Yule process) with speciation rate (birth parameter) *Î»*. He then showed that *N*_{
t
}, the number of species alive at time *t*, follows a geometric distribution with probability mass function (pmf)

*p*_{
n
}(*t*; 1) = Pr{*N*_{
t
}= *n*|*N*_{0} = 1} = *e*^{-Î»t}(1 - *e*^{-Î»t})^{n - 1}Â Â Â (1)

for *n* = 1,2,.... If instead there are initially *n*_{0} species then from standard results (*e.g*. Bailey, 1964) the distribution of *N*_{
t
}is negative binomial with pmf

{p}_{n}(t;{n}_{0})=\mathrm{Pr}\{{N}_{t}=n|{N}_{0}={n}_{0}\}=\left(\begin{array}{c}n\xe2\u02c6\u20191\\ {n}_{0}\xe2\u02c6\u20191\end{array}\right){e}^{\xe2\u02c6\u2019\mathrm{\xce\xbb}{n}_{0}t}{(1\xe2\u02c6\u2019{e}^{\xe2\u02c6\u2019\mathrm{\xce\xbb}t})}^{n\xe2\u02c6\u2019{n}_{0}}\left(2\right)

for *n* = *n*_{0}, *n*_{0} + 1,....

We now consider evolution of genera, and of species within genera, over an epoch between cataclysmic events. Let the time origin be the time of the previous cataclysm, and suppose only a single genus (containing *n*_{0} species) survived that cataclysm. Let *Ï„* be the time of the succeeding cataclysm. Yule assumed that new genera were formed from old in a process analogous to that of speciation, thereby establishing that the time in existence of any genus would follow a truncated exponential distribution, with parameter equal to the rate at which new genera are formed from old. But it is more realistic to assume that a new genus is formed when a speciation within an existing genus is of such a radical form as to qualify the new species as belonging to a completely new genus. Thus the probabilty of a new genus being formed in an infinitesimal interval (*t*, *t* + *dt*) should be proportional to *the existing number of species in all existing genera in the family* (and not to the existing number of genera in the family). We let

*K*_{
t
}denote the number of genera at time *t*, evolved from the pioneeering *n*_{0} species;

*L*_{
t
}denote the number of species at time *t* in all genera, evolved from the pioneeering *n*_{0} species; and

*N*_{
t
}denote the number of species in the pioneering genus at time *t*.

We assume that speciations (within a genus) occur at the rate *Î»* and new genera are formed from existing species at the rate *Î³*. Then to order *o*(*dt*) the following state transitions (of *K*_{
t
}, *L*_{
t
}, *N*_{
t
}) can occur in (*t*, *t* + *dt*):

(*k*, *l* - 1, *n* - 1) â†’ (*k*, *l*, *n*) with probability *Î»*(*n* - 1)*dt*

(*k*, *l* - 1, *n*) â†’ (*k*, *l*, *n*) with probability *Î»*(*l* - 1 - *n*)*dt*

(*k* - 1, *l* - 1, *n*) â†’ (*k*, *l*, *n*) with probability *Î³*(*l* - 1)*dt*

(*k*, *l*, *n*) â†’ (*k*, *l*, *n*) with probability 1 - (*Î»* + *Î³*)*ldt*.

Letting *p*_{k, l, n}(*t*) = P(*K*_{
t
}= *k*, *L*_{
t
}= *l*, *N*_{
t
}= *n*), the following differential-difference equations can be established from the above:

\begin{array}{ccc}\frac{d}{dt}{p}_{k,l,n}(t)& =& \mathrm{\xce\xbb}(n\xe2\u02c6\u20191){p}_{k,l\xe2\u02c6\u20191,n\xe2\u02c6\u20191}(t)+\mathrm{\xce\xbb}(l\xe2\u02c6\u20191\xe2\u02c6\u2019n){p}_{k,l\xe2\u02c6\u20191,n}(t)\\ +\mathrm{\xce\xb3}(l\xe2\u02c6\u20191){p}_{k\xe2\u02c6\u20191,l\xe2\u02c6\u20191,n}(t)\xe2\u02c6\u2019(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})l{p}_{k,l,n}(t).\end{array}\left(3\right)

Using the generating function

\mathrm{\xce\xa6}(x,y,z;t)={\displaystyle \underset{k=1}{\overset{\mathrm{\xe2\u02c6\u017e}}{\xe2\u02c6\u2018}}{\displaystyle \underset{l=1}{\overset{\mathrm{\xe2\u02c6\u017e}}{\xe2\u02c6\u2018}}{\displaystyle \underset{n=1}{\overset{\mathrm{\xe2\u02c6\u017e}}{\xe2\u02c6\u2018}}{p}_{k,l,n}}(t){x}^{k}{y}^{l}{z}^{n}}},\left(4\right)

multiplying (3) by *x*^{k}*y*^{l}*z*^{n}and summing yields the following partial differential equation

Î¦_{
t
}= *y*(*Î»y* + *Î³xy* - (*Î»* + *Î³*)) Î¦_{
y
}+ *Î»yz*(*z* - 1) Î¦_{
z
}, Â Â Â (5)

which can be solved by the method of characteristics (*e.g*. Bailey,[6]) with initial condition *Ï•*(*x*, *y*, *z*; 0) = x{y}^{{n}_{0}}{z}^{{n}_{0}}. From the solution the generating functions of *K*_{
t
}, *L*_{
t
}and *N*_{
t
}can be derived. They are

{\mathrm{\xce\xa6}}_{K}(x,t)=E({x}^{{K}_{t}})=x{\left\{\frac{p(t)}{1\xe2\u02c6\u2019x[1\xe2\u02c6\u2019p(t)]}\right\}}^{{n}_{0}},\left(6\right)

{\mathrm{\xce\xa6}}_{L}(y,t)=E({y}^{{L}_{t}})={\left\{\frac{y{e}^{\xe2\u02c6\u2019(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})t}}{1\xe2\u02c6\u2019y[1\xe2\u02c6\u2019{e}^{\xe2\u02c6\u2019(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})t}]}\right\}}^{{n}_{0}},\left(7\right)

{\mathrm{\xce\xa6}}_{N}(z,t)=E({z}^{{N}_{t}})={\left\{\frac{z{e}^{\xe2\u02c6\u2019\mathrm{\xce\xbb}t}}{1\xe2\u02c6\u2019z(1\xe2\u02c6\u2019{e}^{\xe2\u02c6\u2019\mathrm{\xce\xbb}t})}\right\}}^{{n}_{0}},\left(8\right)

where

p(t)=\frac{(\mathrm{\xce\xbb}+\mathrm{\xce\xb3}){e}^{\xe2\u02c6\u2019(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})t}}{\mathrm{\xce\xb3}+\mathrm{\xce\xbb}{e}^{\xe2\u02c6\u2019(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})t}}.\left(9\right)

From this it is clear that both the total number of species, *L*_{
t
}, and the number of species in the pioneering genus, *N*_{
t
}, have negative binomial distributions (with parameters *n*_{0} and *e*^{-(Î»+ Î³)t}and n_{0} and *e*^{-Î»t}respectively); while the number of genera *K*_{
t
}has a distribution related to the negative binomial â€“ precisely *K*_{
t
}+ *n*_{0} - 1 has a negative binomial distribution with parameters *n*_{0} and *p*(*t*). The expected number of genera at time *t* is

E({K}_{t})=1+\frac{{n}_{0}\mathrm{\xce\xb3}}{\mathrm{\xce\xbb}+\mathrm{\xce\xb3}}\left[{e}^{(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})t}\xe2\u02c6\u20191\right].\left(10\right)

It can be shown (see Appendix) that the times of formation of derived genera constitute an *order statistic process*. This means that they can be considered as the order statisics of a collection of independent, identically distributed (iid) random variables. From this it is shown that at any fixed time *Ï„*, the times *t*_{1}, *t*_{2},...,*t*_{
k
}that the derived genera have been in existence are iid random variables with probability density function (pdf)

\begin{array}{cc}{f}_{k}(t)=\frac{(\mathrm{\xce\xbb}+\mathrm{\xce\xb3}){e}^{\xe2\u02c6\u2019(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})t}}{1\xe2\u02c6\u2019{e}^{\xe2\u02c6\u2019(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})\mathrm{\xcf\u201e}}},& 0<t<\mathrm{\xcf\u201e}\end{array}.\left(11\right)

By summing (3) over *k* and *l* one can show that *N*_{
t
}is a pure birth process with birthrate *Î»*; and by summing over *k* and *n* that *L*_{
t
}is a pure birth process with birthrate *Î»* + *Î³*. From the fact that a pure birth process is an order statistic process it can be shown (see Appendix) that at time *Ï„* the times since establishment of all non-pioneering species in the pioneering *genus* are independently distributed random variables, with a truncated exponential distribution with pdf

\begin{array}{cc}{f}_{N}(t)=\frac{\mathrm{\xce\xbb}{e}^{\xe2\u02c6\u2019\mathrm{\xce\xbb}t}}{1\xe2\u02c6\u2019{e}^{\xe2\u02c6\u2019\mathrm{\xce\xbb}\mathrm{\xcf\u201e}}},& 0<t<\mathrm{\xcf\u201e};\end{array}\left(12\right)

and that the times since establishment of all non-pioneering species in the pioneering *family* are independently distributed random variables, with a truncated exponential distribution with pdf

\begin{array}{cc}{f}_{L}(t)=\frac{(\mathrm{\xce\xbb}+\mathrm{\xce\xb3}){e}^{\xe2\u02c6\u2019(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})t}}{1\xe2\u02c6\u2019{e}^{\xe2\u02c6\u2019(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})\mathrm{\xcf\u201e}}},& 0<t<\mathrm{\xcf\u201e}\end{array}.\left(13\right)

Note the fact that *f*_{
L
}(*t*) â‰¡ *f*_{
K
}(*t*) *i.e*. the marginal distribution of the time since establishment of a derived genus in the family is the same as that of a derived species in the family.

Consider now the case when *Ï„* is the time of the first cataclysm since the appearance of the pioneering genus. The size distribution of all derived (non-pioneering) genera at the time of the cataclysm can be obtained by integrating the geometric pmf *p*_{
n
}(*t*; 1) in (1) with respect to the truncated exponential distribution *f*_{
K
}(*t*) between 0 and *Ï„*. This yields the pmf

\begin{array}{lll}{q}_{n}^{deriv}\hfill & =\hfill & {\displaystyle {\xe2\u02c6\xab}_{0}^{\mathrm{\xcf\u201e}}{p}_{n}(t;1){f}_{K}(t)dt}\hfill \\ =\hfill & \frac{1+\mathrm{\xce\xb3}/\mathrm{\xce\xbb}}{1+{e}^{\xe2\u02c6\u2019(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})\mathrm{\xcf\u201e}}}[B(2+\mathrm{\xce\xb3}/\mathrm{\xce\xbb},n)\xe2\u02c6\u2019{B}_{{e}^{\xe2\u02c6\u2019\mathrm{\xce\xbb}\mathrm{\xcf\u201e}}}(2+\mathrm{\xce\xb3}/\mathrm{\xce\xbb},n)],\hfill \end{array}\left(14\right)

where

\begin{array}{cc}B(a,b)=\frac{\mathrm{\xce\u201c}(a)\mathrm{\xce\u201c}(b)}{\mathrm{\xce\u201c}(a+b)},& {B}_{x}(a,b)={\displaystyle {\xe2\u02c6\xab}_{0}^{x}{z}^{a\xe2\u02c6\u20191}{(1\xe2\u02c6\u2019z)}^{b\xe2\u02c6\u20191}dz}\end{array}

are the *beta function* and *incomplete beta functions*, respectively. Alternatively the term in square brackets can be expressed in terms of the cumulative distribution function (cdf) *F*(*x*; *a*, *b*) of the *beta distribution* with parameters *a* and *b* leading to

{q}_{n}^{deriv}=\frac{(1+\mathrm{\xce\xb3}/\mathrm{\xce\xbb})B(2+\mathrm{\xce\xb3}/\mathrm{\xce\xbb},n)}{1\xe2\u02c6\u2019{e}^{\xe2\u02c6\u2019(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})\mathrm{\xcf\u201e}}}\left[1\xe2\u02c6\u2019F({e}^{\xe2\u02c6\u2019\mathrm{\xce\xbb}\mathrm{\xcf\u201e}};2+\mathrm{\xce\xb3}/\mathrm{\xce\xbb},n)\right].\left(15\right)

This can be readily computed using standard statistical software.

The distribution of the size of the pioneering genus at time *Ï„* has pmf {q}_{n}^{pion} = *p*_{
n
}(*Ï„*; *n*_{0}) where *p*_{
n
}is negative binomial pmf given by (2). The distribution of the size of all existing genera at time *Ï„* is simply a mixture of {q}_{n}^{pion} and {q}_{n}^{deriv}. Precisely

{q}_{n}={\mathrm{\xcf\u20ac}}_{K}(\mathrm{\xcf\u201e}){q}_{n}^{pion}+[1\xe2\u02c6\u2019{\mathrm{\xcf\u20ac}}_{K}(\mathrm{\xcf\u201e})]{q}_{n}^{deriv},\left(16\right)

where *Ï€*_{
K
}(*Ï„*) is the probability that a genus in existence at time *Ï„* is the pioneering genus, *i.e*.

{\mathrm{\xcf\u20ac}}_{K}(\mathrm{\xcf\u201e})=E\left(\frac{1}{{K}_{\mathrm{\xcf\u201e}}}\right)={\displaystyle {\xe2\u02c6\xab}_{0}^{1}\frac{{\mathrm{\xce\xa6}}_{K}(s,\mathrm{\xcf\u201e})}{s}}ds,\left(17\right)

which can be evaluated as

{\mathrm{\xcf\u20ac}}_{K}(\mathrm{\xcf\u201e})=\frac{(\mathrm{\xce\xbb}+\mathrm{\xce\xb3}){e}^{\xe2\u02c6\u2019(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})\mathrm{\xcf\u201e}}}{\mathrm{\xce\xb3}[1\xe2\u02c6\u2019{e}^{\xe2\u02c6\u2019(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})\mathrm{\xcf\u201e}}]}\mathrm{log}\left(\frac{\mathrm{\xce\xb3}+\mathrm{\xce\xbb}{e}^{\xe2\u02c6\u2019(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})\mathrm{\xcf\u201e}}}{(\mathrm{\xce\xbb}+\mathrm{\xce\xb3}){e}^{\xe2\u02c6\u2019(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})\mathrm{\xcf\u201e}}}\right).\left(18\right)

Note that as *Ï„* â†’ âˆž, *Ï€*_{
K
}(*Ï„*) â†’ 0 and

{q}_{n}\xe2\u2020\u2019\frac{(\mathrm{\xce\xb3}/\mathrm{\xce\xbb}+1)\mathrm{\xce\u201c}(\mathrm{\xce\xb3}/\mathrm{\xce\xbb}+2)\mathrm{\xce\u201c}(n)}{\mathrm{\xce\u201c}(\mathrm{\xce\xb3}/\mathrm{\xce\xbb}+n+2)}.\left(19\right)

This distribution was obtained by Yule[1] and is now known as the *Yule distribution*; for this distribution *q*_{
n
}behaves asymptotically like a power-law, *i.e*.,

*q*_{
n
}~ (*Î³*/*Î»* + 1)*Î“*(*Î³*/*Î»* + 2) Ã— *n*^{-(2 + Î³/Î»)}

as *n* â†’ âˆž, yielding the asymptotic straight line when *q*_{
n
}is plotted against *n* on logarithmic axes. We note in passing that setting *Î³* = 0 in (19) does *not* yield the size distribution (as *Ï„* â†’ âˆž) of a single genus, since when *Î³* = 0, *Ï€*_{
K
}â‰¡ 1. In this case *N*_{
Ï„
}â†’ âˆž with probability one.

Figure 1 shows the size distribution of pioneering and derived genera, along with the mixed distribution of all genera, calculated from the above formulae, for different values of *n*_{0} and *Ï„*. They show how the results of Yule [1] need to be modified to take into account the effects of: (a) the evolution of new genera ; (b) pioneering genera of size (*n*_{0}) greater than one; and (c) the time, *Ï„*, until cataclysmic extinction. Large values of *Ï„* (right-hand panels), resulting in straight-line plots on the log-log scale, correspond most closely to the situation considered initially by Yule. In this case approximate power-law (fractal) distributions occur. The deviations from such a power-law distribution are greatest when cataclysmic extinction occurs earlier (smaller *Ï„*) and when the number of species in the pioneering genus (*n*_{0}) differs greatly from one (lower panels). The distribution of derived genera (dotted lines) is unaffected by the initial size (*n*_{0}) of the pioneering genus. However the overall size distribution is affected (especially at values immediately above *n*_{0}) because of the fact that the pioneering genus size has support on {*n*_{0}, *n*_{0} + 1,...} while that of derived genera is on {1, 2,...}. This effect becomes less important when a long time elapses before the cataclysmic extinction event (because when *Ï„* is large, *Ï€*_{
K
}(*Ï„*) is smallâ€“derived genera will in probability outnumber the pioneering one).

## Background extinctions only

In this section we consider the size distribution of a fossil genus, starting with a single species (the case of a genus beginning with *n*_{0} species is considered later in this section), subject to speciations at rate *Î»* and background (individual) extinctions occurring independently and at random, at rate *Î¼*.

Thus *N*_{
t
}, the number of species alive *t* time units after the origin of the genus, follows a homogeneous birth and death process. Let *M*_{
t
}denote the total number of species in the genus that have existed by time *t* (*i.e*. *M*_{
t
}= 1 + number of speciations). The size of an extinct genus is a random variable *M*_{
T
}, where *T* itself is a random variable, denoting the time of extinction. Since no speciations can occur in a genus once it is extinct, we have that for *t* â‰¥ *T*, *M*_{
t
}â‰¡ *M*_{
T
}. However *T* may not be finite (*N*_{
t
}> 0 for all *t*). Thus finding the distribution of the size of an extinct genus will involve conditioning on *T* < âˆž (or *N*_{âˆž} = 0). Clearly it is given by the distribution of *M*_{âˆž} conditional on *N*_{âˆž} = 0.

Now let

*p*_{m, n}(*t*) = Pr(*M*_{
t
}= *m*, *N*_{
t
}= *n*). Â Â Â (20)

It was shown by Kendall[7] that *p*_{m, n}satisfies the differential-difference equations

\frac{d}{dt}{p}_{m,n}(t)=\xe2\u02c6\u2019(\mathrm{\xce\xbb}+\mathrm{\xce\xbc})n{p}_{m,n}(t)+\mathrm{\xce\xbb}(n\xe2\u02c6\u20191){p}_{m\xe2\u02c6\u20191,n\xe2\u02c6\u20191}(t)+\mathrm{\xce\xbc}(n+1){p}_{m,n+1}(t)\left(21\right)

with initial condition

*p*_{m, n}(0) = 1 if *m* = *n* = 1; *p*_{m, n}(0) = 0 otherwise.

Let

\mathrm{\xce\xa8}(s,z;t)={\displaystyle \underset{m=1}{\overset{\mathrm{\xe2\u02c6\u017e}}{\xe2\u02c6\u2018}}{\displaystyle \underset{n=0}{\overset{\mathrm{\xe2\u02c6\u017e}}{\xe2\u02c6\u2018}}{p}_{m,n}(t){s}^{m}{z}^{n}}}\left(22\right)

be the generating function for *M*_{
t
}, *N*_{
t
}. Muliplying both sides of (21) by *s*^{m}*z*^{n}and summing over *m* = l,... âˆž; *n* = 0,...,âˆž yields the partial differential equation

Î¨_{
t
}= (*sz*^{2}*Î»* - (*Î»* + *Î¼*)*z* + *Î¼*)Î¨_{
z
}. Â Â Â (23)

This equation was derived and solved by Kendall[7], using the method of characteristics. The solution is (for *Î»* â‰ *Î¼*)

\mathrm{\xce\xa8}(s,z;t)=\frac{\mathrm{\xce\xb2}(sz\xe2\u02c6\u2019\mathrm{\xce\pm})\mathrm{exp}(\mathrm{\xce\xbb}\mathrm{\xce\pm}t)+\mathrm{\xce\pm}(\mathrm{\xce\xb2}\xe2\u02c6\u2019sz)\mathrm{exp}(\mathrm{\xce\xbb}\mathrm{\xce\xb2}t)}{(sz\xe2\u02c6\u2019\mathrm{\xce\pm})\mathrm{exp}(\mathrm{\xce\xbb}\mathrm{\xce\pm}t)+(\mathrm{\xce\xb2}\xe2\u02c6\u2019sz)\mathrm{exp}(\mathrm{\xce\xbb}\mathrm{\xce\xb2}t)},\left(24\right)

where *Î±* = *Î±*(*s*), *Î²* = *Î²*(*s*) are the two (positive) roots of the quadratic equation

*Î»* *x*^{2} - (*Î»* + *Î¼*)*x* + *Î¼s* = 0. Â Â Â (25)

These roots are distinct for 0 â‰¤ *s* â‰¤ 1, except when *Î»* = *Î¼*, where the roots are distinct for 0 â‰¤ *s* â‰¤ 1, but coincide for *s* = 1. We select *Î²*(*s*) to be the smaller root, so that

\mathrm{\xce\xb2}(s)=\frac{1}{2}\left\{1+\frac{\mathrm{\xce\xbc}}{\mathrm{\xce\xbb}}\xe2\u02c6\u2019\sqrt{{\left(1+\frac{\mathrm{\xce\xbc}}{\mathrm{\xce\xbb}}\right)}^{2}\xe2\u02c6\u2019\frac{4\mathrm{\xce\xbc}s}{\mathrm{\xce\xbb}}}\right\}\left(26\right)

and note that *Î±*(1) = max{*Î»*, *Î¼*}/*Î»*, *Î²*(1) = min{*Î»*, *Î¼*}/*Î»* and *Î»*[*Î±*(1) - *Î²*(1)] = |*Î»* - *Î¼*|.

From (24) the individual generating function *Ïˆ*_{
M
}(*s*; *t*) = *E*({s}^{{M}_{t}}) of *M*_{
t
}(and similarly that of *N*_{
t
}) can be derived. Specifically

{\mathrm{\xce\xa8}}_{M}(s;t)=E({s}^{{M}_{t}})=\frac{\mathrm{\xce\xb2}(s\xe2\u02c6\u2019\mathrm{\xce\pm})+\mathrm{\xce\pm}(\mathrm{\xce\xb2}\xe2\u02c6\u2019s){e}^{\xe2\u02c6\u2019\mathrm{\xce\xbb}(\mathrm{\xce\pm}\xe2\u02c6\u2019\mathrm{\xce\xb2})t}}{(s\xe2\u02c6\u2019\mathrm{\xce\pm})+(\mathrm{\xce\xb2}\xe2\u02c6\u2019s){e}^{\xe2\u02c6\u2019\mathrm{\xce\xbb}(\mathrm{\xce\pm}\xe2\u02c6\u2019\mathrm{\xce\xb2})t}}.\left(27\right)

Expanding this in a power-series expansion will yield the size distribution of the number of species which have existed by a finite time *t*. Simple closed-form expressions are not obtainable, but the expansion can be done numerically for specified parameter values using a computer mathematics program such as Maple VII[8]. It is easy to show that

E({M}_{t})={{\mathrm{\xce\xa8}}^{\xe2\u20ac\xb2}}_{M}(1)=1+\frac{\mathrm{\xce\xbb}}{\mathrm{\xce\xbb}\xe2\u02c6\u2019\mathrm{\xce\xbc}}\left[{e}^{(\mathrm{\xce\xbb}\xe2\u02c6\u2019\mathrm{\xce\xbc})t}\xe2\u02c6\u20191\right].\left(28\right)

Note that for *Î»* > *Î¼*, *E*(*M*_{
t
}) â†’ âˆž as *t* â†’ âˆž; while for *Î»* <*Î¼*, *E*(*M*_{
t
}) â†’ *Î¼*/(*Î¼* - *Î»*).

To find the distribution of the size of an extinct genus we consider the distribution of *M*_{
t
}conditional on *N*(*t*) = 0. This has generating function Î©(*s*; *t*) = *E*({s}^{{M}_{t}}|*N*_{
t
}= 0) given by

\mathrm{\xce\copyright}(s;t)={\displaystyle \underset{m=1}{\overset{\mathrm{\xe2\u02c6\u017e}}{\xe2\u02c6\u2018}}pr({M}_{t}=m|{N}_{t}=0){s}^{m}}=\frac{\mathrm{\xce\xa8}(s,0;t)}{pr({N}_{t}=0)}.\left(29\right)

The probabilty of extinction by time *t* in the denominator can be evaluated as Î¨ (1, 0; *t*) (or from standard results on birth and death processes) yielding

\mathrm{\xce\copyright}(s;t)=\left[\frac{\mathrm{\xce\pm}\mathrm{\xce\xb2}(1\xe2\u02c6\u2019{e}^{\xe2\u02c6\u2019\mathrm{\xce\xbb}(\mathrm{\xce\pm}\xe2\u02c6\u2019\mathrm{\xce\xb2})t})}{\mathrm{\xce\pm}\xe2\u02c6\u2019\mathrm{\xce\xb2}{e}^{\xe2\u02c6\u2019\mathrm{\xce\xbb}(\mathrm{\xce\pm}\xe2\u02c6\u2019\mathrm{\xce\xb2})t}}\right]\left[\frac{\mathrm{max}\{\mathrm{\xce\xbb},\mathrm{\xce\xbc}\}\xe2\u02c6\u2019\mathrm{min}\{\mathrm{\xce\xbb},\mathrm{\xce\xbc}\}{e}^{\xe2\u02c6\u2019|\mathrm{\xce\xbb}\xe2\u02c6\u2019\mathrm{\xce\xbc}|t}}{\mathrm{\xce\xbc}(1\xe2\u02c6\u2019{e}^{\xe2\u02c6\u2019|\mathrm{\xce\xbb}\xe2\u02c6\u2019\mathrm{\xce\xbc}|t})}\right],\left(30\right)

for *Î»* â‰ *Î¼*, and

\mathrm{\xce\copyright}(s;t)=\left[\frac{\mathrm{\xce\pm}\mathrm{\xce\xb2}(1\xe2\u02c6\u2019{e}^{\xe2\u02c6\u2019\mathrm{\xce\xbb}(\mathrm{\xce\pm}\xe2\u02c6\u2019\mathrm{\xce\xb2})t})}{\mathrm{\xce\pm}\xe2\u02c6\u2019\mathrm{\xce\xb2}{e}^{\xe2\u02c6\u2019\mathrm{\xce\xbb}(\mathrm{\xce\pm}\xe2\u02c6\u2019\mathrm{\xce\xb2})t}}\right]\left[\frac{\mathrm{\xce\xbb}t+1}{\mathrm{\xce\xbb}t}\right]\left(31\right)

when *Î»* = *Î¼*.

Since once a genus is extinct it remains extinct forever, the size distribution

{q}_{m}^{\xe2\u20ac}\begin{array}{c}\underset{\xc2\xaf}{\underset{\xc2\xaf}{\text{def}}}\end{array}\mathrm{Pr}\{{M}_{\mathrm{\xe2\u02c6\u017e}}=m|{N}_{\mathrm{\xe2\u02c6\u017e}}=0\}\left(32\right)

of an extinct fossil genus can be found by letting *t* â†’ âˆž in the generating function Î©(*s*; *t*) above. Since *Î±*(*s*) â‰¥ *Î²*(*s*), with the inequality strict for 0 â‰¤ *s* < 1, we have *e*^{-Î»(Î±-Î²)t}â†’ 0 as *t* â†’ âˆž. Thus if we let *t* â†’ âˆž in the generating function above, we deduce that for all *Î»* > 0 and *Î¼* > 0,

\begin{array}{llll}\mathrm{\xce\copyright}(s;\mathrm{\xe2\u02c6\u017e})\hfill & =\hfill & {\displaystyle \underset{m=1}{\overset{\mathrm{\xe2\u02c6\u017e}}{\xe2\u02c6\u2018}}{q}_{m}^{\xe2\u20ac}{s}^{m}}=\frac{\mathrm{max}\{\mathrm{\xce\xbb},\mathrm{\xce\xbc}\}\mathrm{\xce\xb2}(s)}{\mathrm{\xce\xbc}}\hfill & \left(33\right)\hfill \\ =\hfill & \frac{(\mathrm{\xce\xbb}+\mathrm{\xce\xbc})}{2\mathrm{min}\{\mathrm{\xce\xbb},\mathrm{\xce\xbc}\}}\left\{1\xe2\u02c6\u2019\sqrt{1\xe2\u02c6\u2019\frac{4\mathrm{\xce\xbb}\mathrm{\xce\xbc}s}{{(\mathrm{\xce\xbb}+\mathrm{\xce\xbc})}^{2}}}\right\}.\hfill & \left(34\right)\hfill \end{array}

Using the binomial theorem to expand the square root in (34) yields the pmf {q}_{m}^{\xe2\u20ac} for the size of an extinct fossil genus. Where *m* â‰¥ *n*_{0} = 1,

\begin{array}{llll}{q}_{m}^{\xe2\u20ac}\hfill & =\hfill & \frac{(\mathrm{\xce\xbb}+\mathrm{\xce\xbc})}{\mathrm{min}\{\mathrm{\xce\xbb},\mathrm{\xce\xbc}\}}\frac{(2m\xe2\u02c6\u20192)!}{(m\xe2\u02c6\u20191)!m!}\frac{{(\mathrm{\xce\xbb}\mathrm{\xce\xbc})}^{m}}{{(\mathrm{\xce\xbb}+\mathrm{\xce\xbc})}^{2m}}\hfill & \left(35\right)\hfill \\ ~\hfill & \frac{(\mathrm{\xce\xbb}+\mathrm{\xce\xbc})}{4{\mathrm{\xcf\u20ac}}^{1/2}\mathrm{min}\{\mathrm{\xce\xbb},\mathrm{\xce\xbc}\}{m}^{3/2}}{\left[\frac{4\mathrm{\xce\xbb}\mathrm{\xce\xbc}}{{(\mathrm{\xce\xbb}+\mathrm{\xce\xbc})}^{2}}\right]}^{m}.\hfill & \left(36\right)\hfill \end{array}

We observe that asymptotically *q*_{
m
}decays faster than a power-law, except in the case when *Î»* = *Î¼* when it follows a power law with exponent -3/2.

The expected size of an extinct genus can be found by evaluating the derivative Î©_{
s
}(1; âˆž), yielding

E({M}_{\mathrm{\xe2\u02c6\u017e}}|{N}_{\mathrm{\xe2\u02c6\u017e}}=0)=\{\begin{array}{ll}\mathrm{\xce\xbb}/(\mathrm{\xce\xbb}\xe2\u02c6\u2019\mathrm{\xce\xbc}),\hfill & \mathrm{\xce\xbb}>\mathrm{\xce\xbc};\hfill \\ \mathrm{\xe2\u02c6\u017e}\hfill & \mathrm{\xce\xbb}=\mathrm{\xce\xbc};\hfill \\ \mathrm{\xce\xbc}/(\mathrm{\xce\xbc}\xe2\u02c6\u2019\mathrm{\xce\xbb}),\hfill & \mathrm{\xce\xbb}<\mathrm{\xce\xbc}.\hfill \end{array}\left(37\right)

The case *Î»* = *Î¼* represents a phase transition analogous to the percolation phase transition (Hughes[9], Grimmett[10]). For this case although with probability one the genus goes extinct (*i.e*. *N*_{âˆž} = 0, w.p.1), the expected time for this to happen is infinite.

If there were initially *n*_{0} species in the genus, the expressions for the generating functions (24), (27) and (34) need to be modified by raising the expressions on the right-hand side to the *n*_{0}th power. In particular, if we denote the pmf for the size of an extinct genus by {q}_{m}^{\xe2\u20ac}(*n*_{0}) we have

{\displaystyle \underset{m={n}_{0}}{\overset{\mathrm{\xe2\u02c6\u017e}}{\xe2\u02c6\u2018}}{q}_{m}^{\xe2\u20ac}}({n}_{0}){s}^{n}={\left\{\frac{(\mathrm{\xce\xbb}+\mathrm{\xce\xbc})}{2\mathrm{min}\{\mathrm{\xce\xbb},\mathrm{\xce\xbc}\}}\left[1\xe2\u02c6\u2019\sqrt{1\xe2\u02c6\u2019\frac{4\mathrm{\xce\xbb}\mathrm{\xce\xbc}s}{{(\mathrm{\xce\xbb}+\mathrm{\xce\xbc})}^{2}}}\right]\right\}}^{{n}_{0}}.\left(38\right)

We deduce at once from Eq. 38) that

\begin{array}{cc}{q}_{m}^{\xe2\u20ac}({n}_{0})={\left[\frac{(\mathrm{\xce\xbb}+\mathrm{\xce\xbc})}{2\mathrm{min}\{\mathrm{\xce\xbb},\mathrm{\xce\xbc}\}}\right]}^{{n}_{0}}{\left[\frac{4\mathrm{\xce\xbb}\mathrm{\xce\xbc}}{{(\mathrm{\xce\xbb}+\mathrm{\xce\xbc})}^{2}}\right]}^{m}{Q}_{m}({n}_{0})& (m\xe2\u2030\yen {n}_{0}),\end{array}\left(39\right)

where

{\displaystyle \underset{m={n}_{0}}{\overset{\mathrm{\xe2\u02c6\u017e}}{\xe2\u02c6\u2018}}{Q}_{m}}({n}_{0}){z}^{m}={[1\xe2\u02c6\u2019{(1\xe2\u02c6\u2019z)}^{1/2}]}^{{n}_{0}}.\left(40\right)

The extraction of numerical values for the coefficients *Q*_{
m
}(*n*_{0}) for a modest fixed value of *n*_{0} is not difficult in practice. Alternatively, *Q*_{
m
}(*n*_{0}) can be found by a contour integral argument that we shall not write out here, leading to the formula

{Q}_{m}({n}_{0})=\frac{1}{\mathrm{\xcf\u20ac}}{\displaystyle \underset{\underset{jodd}{\underset{(}{j=1}}}{\overset{{n}_{0}}{\xe2\u02c6\u2018}}\left(\begin{array}{c}{n}_{0}\\ j\end{array}\right)}\mathrm{sin}(j\mathrm{\xcf\u20ac}/2)\frac{\mathrm{\xce\u201c}(j/2+1)\mathrm{\xce\u201c}(m\xe2\u02c6\u2019j/2)}{\mathrm{\xce\u201c}(m+1)}(m\xe2\u2030\yen {n}_{0}).\left(41\right)

In particular, the following simple formula holds for *n*_{0} = 1, 2, 3 or 4:

\begin{array}{cc}{Q}_{m}({n}_{0})=\frac{{n}_{0}(2m\xe2\u02c6\u20192)!}{{2}^{2m\xe2\u02c6\u20191}(m\xe2\u02c6\u20191)!m!}\{1\xe2\u02c6\u2019\frac{({n}_{0}\xe2\u02c6\u20191)({n}_{0}\xe2\u02c6\u20192)}{4(m\xe2\u02c6\u20193/2)}\},& m\xe2\u2030\yen {n}_{0}\end{array}.

From Eqs (39) and (41) we see that for arbitrary fixed *n*_{0} â‰¥ 1,

{q}_{m}^{\xe2\u20ac}({n}_{0})~{\left[\frac{\mathrm{\xce\xbb}+\mathrm{\xce\xbc}}{2\mathrm{min}\{\mathrm{\xce\xbb},\mathrm{\xce\xbc}\}}\right]}^{{n}_{0}}{\left[\frac{4\mathrm{\xce\xbb}\mathrm{\xce\xbc}}{{(\mathrm{\xce\xbb}+\mathrm{\xce\xbc})}^{2}}\right]}^{m}\frac{{n}_{0}}{2{\mathrm{\xcf\u20ac}}^{1/2}{m}^{3/2}}

as *m* â†’ âˆž. The right-hand side of this differs from that of (36) only by a multiplicative constant, and for all *n*_{0} â‰¥ 1 asymptotically {q}_{m}^{\xe2\u20ac}(*n*_{0}) decays faster than a power law except in the case *Î»* = *Î¼*, when it follows a power law with exponent -3/2.

Fig. 2 shows the distribution of the size of an extinct genus plotted on logarithmic axes, for two values of *n*_{0} and three values of *Î¼* with *Î»* = 1. In the case *n*_{0} = 1 (left-hand panel), an approximate power-law distribution (straight-line plot) can be seen in the case of equal birth and death rates (*Î»* = *Î¼*, the solid line). When the birth and death rates differ (*Î»* â‰ *Î¼*) there is departure from the power-law with faster decay in probabilities as genus size increases both when *Î»* > *Î¼* and when *Î»* <*Î¼*. In the case when the initial size *n*_{0} of the pioneering genus exceeds one (right-hand panel), similar results pertain asymptotically (large genus sizes), but perturbations in the size distribution occur at the lower end (around *n*_{0}).

## Both background and cataclysmic extinctions

We have very limited results in the case. The difficulty lies in the fact that at the time (*Ï„*, say) at which the cataclysmic extinction event occurs, different genera will have been in existence for different lengths of time. Unlike the case discussed in an earlier (no background extinctions) where we established that the times of establishment of new genera formed an order-statistic process, whence it followed that at time *Ï„*, the times in existence of distinct genera constituted iid random variables with a truncated exponential distribution, in the present case (with background extinctions) we have not been able to establish that the times of establishment of new genera constitute an order-statistic process. Thus it has not been possible to determine the size distribution of derived genera, destroyed in the cataclysm, since their time in existence is unknown. This is particularly unfortunate, since it seems that in fact for many fossil families both background and cataclysmic extinctions have occurred (Raup and Sepkoski [11]).

The only genus for which the time in existence is known is the pioneering genus. The pgf of the size of this genus is given by {\left[{\mathrm{\xce\xa6}}_{M}\left(s;\mathrm{\xcf\u201e}\right)\right]}^{{n}_{0}} where Î¦_{
M
}is defined in (27). This cannot be expanded in terms of simple functions to obtain explicit probabilities for sizes, although of course it can always be done numerically for specific parameter values. The expected size of the pioneering genus is

E({M}_{pion})={n}_{0}\left(1+\frac{\mathrm{\xce\xbb}}{\mathrm{\xce\xbb}\xe2\u02c6\u2019\mathrm{\xce\xbc}}\left[{e}^{(\mathrm{\xce\xbb}\xe2\u02c6\u2019\mathrm{\xce\xbc})\mathrm{\xcf\u201e}}\xe2\u02c6\u20191\right]\right).\left(42\right)

## 1 Size distribution of families

In this section we consider the number of *genera* in the *family* derived from the pioneering species, assuming (as in the second section) that new genera are created by extreme speciations (at probabilistic rate *Î³*) and (as in the third section) that background extinctions occur at the rate *Î¼*.

It can be shown (see Appendix) that the number of genera, *G*_{
t
}, which have existed up to time *t* has a generating function Î¦_{
G
}(*s*; *t*) = *E*({s}^{{G}_{t}}) given by

{\mathrm{\xce\xa6}}_{G}(s;t)=s{\left[\frac{(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})}{\mathrm{\xce\xbb}+\mathrm{\xce\xb3}s}\stackrel{\xcb\u0153}{\mathrm{\xce\xa8}}\left(\frac{\mathrm{\xce\xbb}+\mathrm{\xce\xb3}s}{\mathrm{\xce\xbb}+\mathrm{\xce\xb3}};t\right)\right]}^{{n}_{0}},\left(43\right)

where \stackrel{\xcb\u0153}{\mathrm{\xce\xa8}} is the same as Î¨_{
M
}in (27), but with *Î»* replaced by *Î»* + *Î³*. This can be verified directly in the case *Î¼* = 0 (only cataclysmic extinctions) for which *G*_{
t
}â‰¡ *K*_{
t
}(see second section) with *G*_{
t
}+ *n*_{0} - 1 having a negative binomial distribution. In the more general case the proof is somewhat technical and is relegated to the Appendix. The expected number of genera in the family can easily be determined from (43) as

E({G}_{t})=1+{n}_{0}\frac{\mathrm{\xce\xb3}}{\mathrm{\xce\xbb}+\mathrm{\xce\xb3}\xe2\u02c6\u2019\mathrm{\xce\xbc}}\left({e}^{(\mathrm{\xce\xbb}+\mathrm{\xce\xb3}\xe2\u02c6\u2019\mathrm{\xce\xbc})t}\xe2\u02c6\u20191\right).\left(44\right)

If, following a cataclysmic event from which *n*_{0} species survived, a subsequent cataclysm occurred *Ï„* time units later, the size distribution of the family (number of genera) derived from these *n*_{0} pioneering species, would have pgf Î¦_{
G
}(*s*; *Ï„*). While no simple expansion of this is possible it can be done numerically. Some examples are shown in Fig. 3. The distributions show considerable deviation from a power law (straight line in logarithmic plots). They appear similar to the corresponding distributions of number of species in a genus (Fig. 1, top row) for smaller values of *Ï„*, but are further from the power-law form for larger *Ï„*. Thus it would appear that under the birth-and-death model power-law (fractal-like) size distributions are less likely to occur at higher taxonomic levels.

## Monotypic taxa

One characteristic of interest in the empirical study of lineages is the proportion of monotypic taxa. Przeworski and Wall[5] compared the proportions of monospecific genera and of monogeneric families observed in the fossil record with results from a simulation of a birth-and-death process model. In this section we compute probabilities of such monotypic taxa. We consider the cases of (1) only background extinctions; and (2) only cataclysmic extinctions.

### Only background extinctions

For a genus in existence for *t* time units, the probability of it having only ever contained one species by that time is

\mathrm{Pr}({M}_{t}=1)={{\mathrm{\xce\xa8}}^{\xe2\u20ac\xb2}}_{M}(0;t)=\underset{s\xe2\u2020\u20190}{\mathrm{lim}}\frac{{\mathrm{\xce\xa8}}_{M}(s;t)}{s}=\frac{\mathrm{\xce\xbb}{e}^{\xe2\u02c6\u2019(\mathrm{\xce\xbb}+\mathrm{\xce\xbc})t}+\mathrm{\xce\xbc}}{\mathrm{\xce\xbb}+\mathrm{\xce\xbc}}\left(45\right)

where Î¨_{
M
}is as in (27). Since all extinct fossil genera are finite in size, the probability of such a genus being monospecific is (from the results in fourth section)

Pr(monospecificgenus)=\text{P}r({M}_{\mathrm{\xe2\u02c6\u017e}}=1|{M}_{\mathrm{\xe2\u02c6\u017e}}<\mathrm{\xe2\u02c6\u017e})=\{\begin{array}{ll}\frac{\mathrm{\xce\xbc}}{\mathrm{\xce\xbb}+\mathrm{\xce\xbc}},\hfill & \mathrm{\xce\xbb}\xe2\u2030\xa4\mathrm{\xce\xbc},\hfill \\ \frac{\mathrm{\xce\xbb}}{\mathrm{\xce\xbb}+\mathrm{\xce\xbc}},\hfill & \mathrm{\xce\xbb}>\mathrm{\xce\xbc}.\hfill \end{array}\left(46\right)

Note that this is never less than one half (with this minimum value occurring when *Î»* = *Î¼*), so in the absence of cataclysmic extinctions, one should expect at least half of all extinct genera to be monospecific.

Consider now the distribution of the number of *genera* derived from a pioneering genus with *n*_{0} species. Again since all observed extinct families will be of finite size, the probability of such a fossil family being monogeneric is

\mathrm{Pr}({G}_{\mathrm{\xe2\u02c6\u017e}}=1|{G}_{\mathrm{\xe2\u02c6\u017e}}<\mathrm{\xe2\u02c6\u017e})=\{\begin{array}{ll}{\left(\frac{\mathrm{\xce\xbb}+\mathrm{\xce\xb3}}{\mathrm{\xce\xbc}}\right)}^{{n}_{0}}{{\mathrm{\xce\xa6}}^{\xe2\u20ac\xb2}}_{G}(0,\mathrm{\xe2\u02c6\u017e}),\hfill & \mathrm{\xce\xbb}+\mathrm{\xce\xb3}>\mathrm{\xce\xbc},\hfill \\ {{\mathrm{\xce\xa6}}^{\xe2\u20ac\xb2}}_{G}(0,\mathrm{\xe2\u02c6\u017e}),\hfill & \mathrm{\xce\xbb}+\mathrm{\xce\xb3}\xe2\u2030\xa4\mathrm{\xce\xbc}.\hfill \end{array}\left(47\right)

where

{{\mathrm{\xce\xa6}}^{\xe2\u20ac\xb2}}_{G}(0,\mathrm{\xe2\u02c6\u017e})=\frac{\xe2\u02c6\u201a}{\xe2\u02c6\u201as}{\mathrm{\xce\xa6}}_{G}(s,\mathrm{\xe2\u02c6\u017e}){|}_{s=0}={\left[\left(\frac{(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})}{\mathrm{\xce\xbb}}\right){\mathrm{\xce\xa8}}_{L}\left(\frac{\mathrm{\xce\xbb}}{\mathrm{\xce\xbb}+\mathrm{\xce\xb3}};\mathrm{\xe2\u02c6\u017e}\right)\right]}^{{n}_{0}}

using (43). Thus, using (34), when *Î»* + *Î³* > *Î¼*

Pr(monogenericfamily)={\left[\frac{\mathrm{\xce\xbb}+\mathrm{\xce\xb3}}{2\mathrm{\xce\xbb}\mathrm{\xce\xbc}}\left(\mathrm{\xce\xbb}+\mathrm{\xce\xb3}+\mathrm{\xce\xbc}\xe2\u02c6\u2019\sqrt{{(\mathrm{\xce\xbb}+\mathrm{\xce\xb3}+\mathrm{\xce\xbc})}^{2}\xe2\u02c6\u20194\mathrm{\xce\xbb}\mathrm{\xce\xbc}}\right)\right]}^{{n}_{0}};

and when *Î»* + *Î³* â‰¤ *Î¼*, the right hand side is modified by the fraction (*Î»* + *Î³*)/(2*Î»Î¼*) being replaced by 1/(2*Î»*).

Comparing the probability of a monospecific genus with that of a mono-generic family is complicated in general because of the number of parameters. But one can show that with *n*_{0} = 1, the probability of a monogeneric family always exceeds that of a monospecific genus if the rate of formation of new genera is suitably small - *i.e*. if 0 <*Î³* <*Î³*_{0}, for some positive *Î³*_{0} (depending on *Î»* and *Î¼*). In this case of course the probability of a monogeneric family will also exceed 0.5.

### Only cataclysmic extinctions

If a cataclysmic extinction event occurs at time *Ï„*, the probabilities of a monotypic genus and of a monogeneric family can be found easily from the results of the second section using the explicit expressions for the generating functions of the number of species *N*_{
Ï„
}, (8); and for the number of genera *L*_{
Ï„
}, (6). Specifically if there is initially a single species in the genus the probability that it is monospecific at the time of extinction is

Pr(monospecific genus) = Pr(*N*_{
Ï„
}= 1) = e^{-Î»Ï„}, Â Â Â (48)

which is simply the probabilty of no speciations in (0, *Ï„*). In contrast the probabilty of a monogeneric family is

\mathrm{Pr}(monogenericfamily)=Pr({K}_{\mathrm{\xcf\u201e}}=1)={[p(\mathrm{\xcf\u201e})]}^{{n}_{0}}={\left[\frac{(\mathrm{\xce\xbb}+\mathrm{\xce\xb3}){e}^{\xe2\u02c6\u2019(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})\mathrm{\xcf\u201e}}}{\mathrm{\xce\xb3}+\mathrm{\xce\xbb}{e}^{\xe2\u02c6\u2019(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})\mathrm{\xcf\u201e}}}\right]}^{{n}_{0}}.\left(49\right)

Comparing the right-hand sides of the above two equations, one can show that provided *Î³* <*Î»*/*n*_{0} then Pr(monogeneric family) > Pr(monospecific genus) for *Ï„* less than some threshold value *Ï„*_{0}, say; but for *Ï„* > *Ï„*_{0} the inequality is reversed. Thus as with the case of only background extinctions, monogeneric fossil families should be more common than monospecific fossil genera when the inter-cataclysm period is short. However if the inter-cataclysm period is longer the situation may be reversed.

## Concluding remarks

In the paper a number of analytic results on the size distributions of genera and families and on the probabilities of monospecific taxa have been derived under the assumption of a simple homogeneous birth-and-death model and various extinction scenarios. The results are incomplete due to the complexity of the analysis, especially in the case when both cataclysmic and background extinctions can occur. However it is hoped that there are sufficient results to enable testing of the birth-and-death model using empirical taxon size distributions obtained from the fossil record.

Undoubtedly more complex plausible extinction scenarios than the two extremes discussed in this paper could be considered. For example one could consider major extinction events which resulted in the destruction of a significant proportion (but not all) of species within a genus. However realistically formulating a model for this, not to mention its subsequent analysis, seems to present a formidable task.

One could also consider the size distribution of taxa existing over more than one inter-cataclysmic epoch. In this case one would need to consider mixtures of the distributions, using different (but assumed known) values of *Ï„*. In principle this is not difficult to do. If the durations of inter-cataclysmic epoch were not known one could consider *Ï„* as a random variable and consider the resulting infinite mixture. As a null model for catclysmic extinction events, it seems reasonable to assume that they occur independently at random, so that the time between two successive events would have an exponential distribution. An overall distribution for the size of a taxon could then be obtained by integrating the results obtained in the earlier sections with respect to an exponential density. This has been considered in another paper (Hughes and Reed[12]) where it is shown that, under certain conditions, the resulting size distributions exhibit fractal-like behaviour.

## Appendix

A point process {*X*_{
t
}, *t* â‰¥ 0} is said to be an *order statistic process* (Feigin[13]) if conditional on *X*_{
Ï„
}- *X*_{0} = *k* the successive jump times (times of events) *T*_{1}, *T*_{2},...,*T*_{
k
}are distributed as the order statistics of *k* independent, identically distributed random variables with support on [0, *Ï„*]. The simplest example is when {*X*_{
t
}} is a Poisson process, for which conditional on *X*_{
Ï„
}- *X*_{0} = *k*, it is well known that the event times *T*_{1}, *T*_{2} ,..., *T*_{
k
}have the same distribution as the order statistics of of *k* independent, uniformly distributed random variables on [0, *Ï„*].

For a given order statistic process the order statistic distribution can be shown (Feigin[13](Theorem 2)) to have cdf

F(t)=\frac{m(t)\xe2\u02c6\u2019m(0)}{m(\mathrm{\xcf\u201e})\xe2\u02c6\u2019m(0)},0\xe2\u2030\xa4y\xe2\u2030\xa4\mathrm{\xcf\u201e},\left(50\right)

where *m*(*t*) = *E*(*X*_{
t
}).

Puri[14] (Theorem 8) gives conditions for a non-homogeneous birth process, with birth rates *Î¸*_{
i
}(*t*), to be an order statistic process. For the process {*K*_{
t
}} (the number of genera) in second section, the birth rates *Î¸*_{
k
}(*t*) are given by

*Î¸*_{
k
}(*t*)*dt* = Pr (*K*(*t* + *dt*) = *k* + 1|*K*(*t*) = *k*)*dt* + *o*(*dt*). Â Â Â (51)

If we sum over *l* and *n* in (3) we find that with *p*_{
k
}(*t*) = Pr{*K*_{
t
}= *k*},

\begin{array}{llll}\frac{d}{dt}{p}_{k}(t)\hfill & =\hfill & \mathrm{\xce\xb3}E({L}_{t}|{K}_{t}=k\xe2\u02c6\u20191){p}_{k\xe2\u02c6\u20191}(t)\xe2\u02c6\u2019\mathrm{\xce\xb3}E({L}_{t}|{K}_{t}=k){p}_{k}(t)\hfill & \left(52\right)\hfill \\ =\hfill & {\mathrm{\xce\xb8}}_{k\xe2\u02c6\u20191}(t){p}_{k\xe2\u02c6\u20191}(t)\xe2\u02c6\u2019{\mathrm{\xce\xb8}}_{k}(t){p}_{k}(t)\hfill & \left(53\right)\hfill \end{array}

so that *K*_{
t
}does evolve under a non-homogeneous birth process, with birth rates

*Î¸*_{
k
}(*t*) = *Î³E*(*L*_{
t
}|*K*_{
t
}= *k*). Â Â Â (54)

We now calculate *Î¸*_{
k
}(*t*) explicitly. From Eq. (6),

{p}_{k}(t)=\frac{{({n}_{0})}_{k\xe2\u02c6\u20191}p{(t)}^{{n}_{0}}}{(k\xe2\u02c6\u20191)!}{[1\xe2\u02c6\u2019p(t)]}^{k\xe2\u02c6\u20191}\left(55\right)

with *p*(*t*) = [(*Î»* + *Î³*)*e*^{-(Î» + Î³)t}]/[*Î³* + *Î»e*^{-(Î» + Î³)t}] and we note for later use that

\frac{{p}^{\xe2\u20ac\xb2}(t)}{p(t)}=\frac{\mathrm{\xce\xb3}(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})}{\mathrm{\xce\xb3}+\mathrm{\xce\xbb}{e}^{\xe2\u02c6\u2019(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})t}}.

Since *p*_{0}(*t*) = 0, we have

{\mathrm{\xce\xb8}}_{1}(t)=\xe2\u02c6\u2019\frac{{{p}^{\xe2\u20ac\xb2}}_{1}(t)}{{p}_{1}(t)}=\xe2\u02c6\u2019\frac{{n}_{0}{p}^{\xe2\u20ac\xb2}(t)}{p(t)}=\frac{\mathrm{\xce\xb3}(\mathrm{\xce\xbb}+\mathrm{\xce\xb3}){n}_{0}}{\mathrm{\xce\xb3}+\mathrm{\xce\xbb}{e}^{\xe2\u02c6\u2019(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})t}}.\left(56\right)

For *k* â‰¥ 1 we have from (53) a difference equation to solve for *Î¸*_{
k
}(*t*):

(*k* - 1)*Î¸*_{k - 1}(*t*) - [1 - *p*(*t*)](*n*_{0} + *k* - 2)*Î¸*_{
k
}(*t*) = (*n*_{0} + *k* -2){*n*_{0} [1 - *p*(*t*)] - (*k* - 1)*p*(*t*)}\frac{{p}^{\xe2\u20ac\xb2}(t)}{p(t)}.

By inspection, a solution of this equation is given by

*Î¸*_{
k
}(*t*) = -\frac{{p}^{\xe2\u20ac\xb2}(t)}{p(t)}(*n*_{0} + *k* - 1), *k* â‰¥ 1.

As this solution gives the correct result (56) for *k* = 1 and a first-order linear difference equation needs only one boundary condition to uniquely determine the solution, we have proved that the birth rate is

{\mathrm{\xce\xb8}}_{k}(t)=\frac{\mathrm{\xce\xb3}(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})({n}_{0}+k\xe2\u02c6\u20191)}{\mathrm{\xce\xb3}+\mathrm{\xce\xbb}{e}^{\xe2\u02c6\u2019(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})t}},k\xe2\u2030\yen 1.

Puri's [[14], Theorem 8] condition for an order-statistic process on (0, *Ï„*) requires the existence of a positive, continuous and integrable function, *h*(*t*) and positive constants *L*(*k*) for *k* = 1, 2,..., with *L*(1) = 1 such that

{\mathrm{\xce\xb8}}_{k}(t)\mathrm{exp}\left\{{\displaystyle {\xe2\u02c6\xab}_{0}^{t}[{\mathrm{\xce\xb8}}_{k+1}(u)\xe2\u02c6\u2019{\mathrm{\xce\xb8}}_{k}(u)]du}\right\}=\frac{h(t)L(k+1)}{L(k)},

In the present case this is satisfied with

*h*(*t*) = *n*_{0}*Î³e*^{(Î»+Î³)t}

and *L*(*k*) = (*n*_{0} - 1)_{
k
}/{n}_{0}^{k}. Also from Puri's Theorem 8,

E({K}_{t})=1+{\displaystyle {\xe2\u02c6\xab}_{0}^{t}h(u)du}=1+\frac{{n}_{0}\mathrm{\xce\xb3}}{\mathrm{\xce\xbb}+\mathrm{\xce\xb3}}[{e}^{(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})t}\xe2\u02c6\u20191].

This agrees with the direct derivation of the expectation from the pgf of *K*_{
t
}(10) and enables the computation of the joint distribution of the times of establishment of derived genera as that of the order statistics of a random sample of size *k* from a distribution with cdf

\begin{array}{cc}F(t)=\frac{{e}^{(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})t}\xe2\u02c6\u20191}{{e}^{(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})\mathrm{\xcf\u201e}}\xe2\u02c6\u20191},& 0\xe2\u2030\xa4t\xe2\u2030\xa4\mathrm{\xcf\u201e}.\end{array}\left(57\right)

Thus it follows that at time *Ï„* the times since establishment of all derived genera are independent random variables with the truncated exponential distribution with pdf *f*_{
K
}(*t*) given in (11).

To establish the truncated exponential nature of the distributions (*f*_{
N
}(*t*) and *f*_{
L
}(*t*), given in second section) of the times since establishment of species in respectively the pioneering genus and the pioneering family, is much easier. From the facts (established in the second section) that {*N*_{
t
}} and {*L*_{
t
}} are pure birth processes with both *N*_{
t
}and *L*_{
t
}having negative binomial distributions with *E*(*N*_{
t
}) = *n*_{0}*e*^{Î»t}and E(*L*_{
t
}) = *n*_{0}*e*^{(Î»+Î³)t}, and the well-known fact that a pure birth process is an order statistic process (Feigin[13]), one can easily establish (using (50)) the cdfs of the times since establishment of non-pioneering species in respectively the pioneeing genus and family. The pdfs, *f*_{
N
}(*t*) and *f*_{
L
}(*t*) given in (12) and (13) follow.

To establish the relationship (43) between the generating functions of *G*_{
t
}(the number of genera which have existed by time *t*) and *L*_{
t
}(the number of species which have existed by *t*), first let

*Y*_{
t
}= *L*_{
t
}- *n*_{0} and *Z*_{
t
}= *G*_{
t
}- 1 Â Â Â (58)

denote the numbers of derived species and genera respectively. Since any speciation could have given rise to a new genus with probability *p* = *Î³*/(*Î»*+*Î³*), independently of other speciations, it follows that *Z*_{
t
}|*y* ~ Bin(*y*, *p*) and hence that

\begin{array}{lll}Pr({G}_{t}=g|{L}_{t}=l)\hfill & =\hfill & Pr({Z}_{t}=g\xe2\u02c6\u20191|{Y}_{t}=l\xe2\u02c6\u2019{n}_{0})\hfill \\ =\hfill & \left(\begin{array}{c}l\xe2\u02c6\u2019{n}_{0}\\ g\xe2\u02c6\u20191\end{array}\right){p}^{g\xe2\u02c6\u20191}{q}^{l\xe2\u02c6\u2019{n}_{0}\xe2\u02c6\u2019g+1}\hfill \\ =\hfill & \frac{{p}^{g\xe2\u02c6\u20191}}{(g\xe2\u02c6\u20191)!}{D}_{q}^{(g\xe2\u02c6\u20191)}{q}^{l\xe2\u02c6\u2019{n}_{0}}\hfill \end{array}\left(59\right)

where *q* = 1 - *p* and *D*_{
q
}is the differential operator \frac{d}{dq}. Multiplying by the pmf *f*_{
l
}= P(*L*_{
t
}= *l*) and summing from *l* = *n*_{0} to âˆž yields the marginal pmf of *G*_{
t
}, which can be written

\begin{array}{ccc}\mathrm{Pr}({G}_{t}=g)& =& \frac{{p}^{g\xe2\u02c6\u20191}}{(g\xe2\u02c6\u20191)!}{D}_{q}^{(g\xe2\u02c6\u20191)}{q}^{\xe2\u02c6\u2019{n}_{0}}{\displaystyle \underset{l={n}_{0}}{\overset{\mathrm{\xe2\u02c6\u017e}}{\xe2\u02c6\u2018}}{q}^{l}{f}_{l}}\\ =& \frac{{p}^{g\xe2\u02c6\u20191}}{(g\xe2\u02c6\u20191)!}{D}_{q}^{(g\xe2\u02c6\u20191)}\left[\frac{\stackrel{\xcb\u0153}{\mathrm{\xce\xa8}}(q)}{{q}^{{n}_{0}}}\right]\end{array}

where \stackrel{\xcb\u0153}{\mathrm{\xce\xa8}}(Â·) is the pgf of *L*_{
t
}which is the same as the pgf of *M*_{
t
}(see (27)), but with *Î»* replaced by *Î»* + *Î¼*. Thus

\mathrm{Pr}({G}_{t}=g)=\frac{{p}^{g\xe2\u02c6\u20191}}{(g\xe2\u02c6\u20191)!}{D}_{q}^{(g\xe2\u02c6\u20191)}{\left[\frac{\stackrel{\xcb\u0153}{\mathrm{\xce\xa8}}(q;t)}{q}\right]}^{{n}_{0}},\left(60\right)

where (using (27))

\stackrel{\xcb\u0153}{\mathrm{\xce\xa8}}(q;t)=\frac{{\mathrm{\xce\xb2}}^{\xe2\u20ac\xb2}(q\xe2\u02c6\u2019{\mathrm{\xce\pm}}^{\xe2\u20ac\xb2})+{\mathrm{\xce\pm}}^{\xe2\u20ac\xb2}({\mathrm{\xce\xb2}}^{\xe2\u20ac\xb2}\xe2\u02c6\u2019q){e}^{\xe2\u02c6\u2019(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})({\mathrm{\xce\pm}}^{\xe2\u20ac\xb2}\xe2\u02c6\u2019{\mathrm{\xce\xb2}}^{\xe2\u20ac\xb2})t}}{(q\xe2\u02c6\u2019{\mathrm{\xce\pm}}^{\xe2\u20ac\xb2})+({\mathrm{\xce\xb2}}^{\xe2\u20ac\xb2}\xe2\u02c6\u2019q){e}^{\xe2\u02c6\u2019(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})({\mathrm{\xce\pm}}^{\xe2\u20ac\xb2}\xe2\u02c6\u2019{\mathrm{\xce\xb2}}^{\xe2\u20ac\xb2})t}}.\left(61\right)

with *Î±'* and *Î²'* being the roots of (25) with *Î»* replaced by *Î»*+*Î¼*. The generating function of Î¨_{
G
}(*s*; *t*) can be obtained by multiplying (60) above by *s*^{g}and summing from *g* = 1 to âˆž:

\begin{array}{lll}{\mathrm{\xce\xa8}}_{G}(s;t)\hfill & =\hfill & {\displaystyle \underset{g=1}{\overset{\mathrm{\xe2\u02c6\u017e}}{\xe2\u02c6\u2018}}{s}^{g}}\mathrm{Pr}({G}_{t}=g)\hfill \\ =\hfill & {s{\displaystyle \underset{g=1}{\overset{\mathrm{\xe2\u02c6\u017e}}{\xe2\u02c6\u2018}}\frac{{[sp]}^{g\xe2\u02c6\u20191}}{(g\xe2\u02c6\u20191)!}}\frac{{d}^{g\xe2\u02c6\u20191}}{d{y}^{g\xe2\u02c6\u20191}}{\left[\frac{\stackrel{\xcb\u0153}{\mathrm{\xce\xa8}}(y;t)}{y}\right]}^{{n}_{0}}|}_{y=q}\hfill \\ =\hfill & s{\left[\frac{\stackrel{\xcb\u0153}{\mathrm{\xce\xa8}}(q+sp;t)}{q+sp}\right]}^{{n}_{0}}\hfill \end{array}\left(62\right)

since the penultimate line is simply a Taylor series expansion about *q* of the last line.

Thus we conclude that

{\mathrm{\xce\xa8}}_{G}(s;t)=s{\left[\frac{(\mathrm{\xce\xbb}+\mathrm{\xce\xb3})}{\mathrm{\xce\xbb}+\mathrm{\xce\xb3}s}\stackrel{\xcb\u0153}{\mathrm{\xce\xa8}}\left(\frac{\mathrm{\xce\xbb}+\mathrm{\xce\xb3}s}{\mathrm{\xce\xbb}+\mathrm{\xce\xb3}};t\right)\right]}^{{n}_{0}}.\left(63\right)

## References

Yule GU: A mathematical theory of evolution, based on the conclusions of Dr. J. C. Willis. F.R.S. R Soc Lond, Philos Trans (B). 1924, 213: 21-87.

Raup DM: Mathematical models of cladogenesis. Paleobiology. 1985, 11: 42-52.

Stoyan D: Estimation of transition rates of inhomogeneous birth-death processes with a paleontological application. Elektronische Informationsverabeitung u. Kybernetic. 1980, 16: 647-649.

Patzkowsky ME: A hierarchical branching model of evolutionary radiations. Paleobiology. 1995, 21: 440-460.

Przeworski M, Wall JD: An evaluation of a hierarchical branching process as a model for species diversification. Paleobiology. 1998, 24: 498-511.

Bailey NTJ: The Elements of Stochastic Processes. 1964, J Wiley and Sons: New York

Kendall DG: On the generalized "birth-and-death" process. Ann Math Stats. 1948, 19: 1-15.

Maple VII: 2001, Waterloo Maple Inc: Waterloo, Ontario

Hughes BD: Random Walks and Random Environments, Random Environments. 1966, Oxford University Press, 2:

Grimmett G: Percolation. 1999, Springer-Verlag: Berlin, 2

Raup DM, Sepkoski JJ: Periodicity of extinction in the geologic past. Proc Nat Acad Sci USA. 1984, 81: 801-805. 10.1073/pnas.81.3.801.

Hughes BD, Reed WJ: A problem in paleobiology. arXiv:physics/0211090. 2002

Feigin PD: On the characterization of point processes with the order statistic property. J Appl Prob. 1979, 16: 297-304. 10.2307/3212898.

Puri PS: On the characterization of point processes with the order statistic property without the moment condition. J Appl Prob. 1982, 19: 39-51. 10.2307/3213914.

## Author information

### Authors and Affiliations

### Corresponding author

## Authorsâ€™ original submitted files for images

Below are the links to the authorsâ€™ original submitted files for images.

## Rights and permissions

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## About this article

### Cite this article

Reed, W.J., Hughes, B.D. Theoretical size distribution of fossil taxa: analysis of a null model.
*Theor Biol Med Model* **4**, 12 (2007). https://doi.org/10.1186/1742-4682-4-12

Received:

Accepted:

Published:

DOI: https://doi.org/10.1186/1742-4682-4-12