# Critical and Near-Critical Branching Processes

## I Introduction

Scale-free distributions, or power laws, have been observed in a variety of biological, chemical and physical systems. Such distributions can arise from different underlying mechanisms, but always involve a separation of scales, which forces the distribution to take a standard form. Scale-free distributions are most often observed in the distribution of sizes of events (such as the Gutenberg-Richter law ), the distribution of times between events (e.g., the inter-event interval distribution in neuronal spike trains ), and frequencies. An example of the latter is the well-known and ubiquitous $`1/f`$ noise. Some systems are even more interesting because they seem to exhibit self-organization or self-tuning, concomitant with scale-free behavior as an inherent and robust property of the system, not due to the tuning of a control parameter by the experimenter.

Two systems to which such spontaneous scale-free behavior has been attributed are sandpile models and taxon creation in biological systems. The former has served as the paradigm of “self-organized criticality”(SOC) , while the latter, manifested in the form of near power law shapes of rank-abundance curves, has been advanced as evidence of a fractal geometry of evolution .

A much simpler system where power laws are observed is the random walk . For example, the waiting times $`t`$ for first return to zero of the simple random walk in one dimension (starting at $`x=0`$, at each time step, $`x(t+1)=x(t)\pm 1`$ with equal probability) have a probability distribution $`t^{3/2}`$. Closely related to random walks, branching processes  can also create power law distributions. They have been used to model the dynamics of many systems in a wide variety of disciplines, including demography, genetics, ecology, physiology, chemistry, nuclear physics, and astrophysics. Here, we use a branching process to model the creation and growth of evolutionary taxa, and discuss its application to avalanches in SOC sandpile models.

In Section II, we examine the properties of the Galton-Watson process. We find that this process can generate power laws by appropriate tuning of a control parameter, and examine the dynamics of the system both at the critical point and away from it. In Section III, we apply this branching process model to the taxonomic rank-frequency abundance patterns of evolution, and discuss the universality of their underlying dynamics. Finally, in Section IV, we discuss the implications of our work, including a discussion of the order and control parameters for the branching process and its applications, and suggest further questions.

## II The Branching Process

The Galton-Watson branching process was first introduced in 1874 to explain the disappearance of family names among the British peerage . It is the first branching process in the literature, and also one of the simplest. Consider an organism that replicates. The number of replicants (daughters) it spawns is determined probabilistically, with $`p_i`$ ($`i=0,1,2`$…) being the probability of having $`i`$ daughters. Each daughter replicates (with the same $`p_i`$ as the original organism) and the daughter’s daughters replicate and so on. We are interested in the rank-frequency probability distribution $`P(n)`$ of the total number of organisms descended from this organism plus $`1`$ (for the original organism), i.e., the historical size of the “colony” the ancestral replicant has given rise to. Note that this is equivalent to asking for the probability distribution of the length of a random walk starting from 1 and returning to 0 with step sizes given by $`P(\mathrm{\Delta }n)=p_{i1}`$ ($`i=0,1,2`$…) .

The abundance distribution $`P(n)`$ can be found by defining a generating function

$`F(s)={\displaystyle \underset{i=1}{\overset{\mathrm{}}{}}}P(i)s^i.`$ (1)   
This function satisfies the relationship

$`F(s)=s{\displaystyle \underset{i=0}{\overset{\mathrm{}}{}}}p_i[F(s)]^i,`$ (2)   
from which each $`P(n)`$ can be determined by equating coefficients of the same order in $`s`$ . This result can also be written as

$`P(n)={\displaystyle \frac{1}{n}}Q(n,n1)(k1),`$ (3)   
where $`Q(i,j)`$ is defined as the probability that $`j`$ organisms will give birth to a total of $`i`$ true daughters . However, these approaches are not numerically efficient, as the calculation of $`P(n)`$ for each new value of $`n`$ requires re-calculation of each term in the result.

For our present purposes, we approach the problem in a different manner. Let $`P_{k|j}`$ be the probability that given $`j`$ original organisms, we end up with a total of $`k`$ organisms after all organisms have finished replicating. Obviously,

$`P_{k|j}=0(k<j),`$ (4)   
since it is impossible to have less total organisms than one starts out with, and

$`P_{1|1}=p_0,`$ (5)   
i.e., the probability for one organism to have no daughters. A little less obviously,

$`P_{k|1}`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{k1}{}}}p_jP_{(k1)|j},`$ (6)   
$`P_{k|j}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{k1}{}}}P_{i|1}P_{i|(j1)}(jk>1).`$ (7)   
These equations allow us to use dynamic programming techniques to calculate $`P(n)`$ ($`=P_{n|1}`$), significantly reducing the computational time required. Also, from Eq. (6), we can write

$`{\displaystyle \frac{P_{n|1}}{P_{(n1)|1}}}=p_1+p_2{\displaystyle \frac{P_{(n1)|2}}{P_{(n1)|1}}}+p_3{\displaystyle \frac{P_{(n1)|3}}{P_{(n1)|1}}}+\text{… }.`$ (8)   
Since, for $`n\mathrm{}`$, $`P_{n|j}`$ is uniformly decreasing, we see

$`{\displaystyle \frac{P(n)}{P(n1)}}={\displaystyle \frac{P_{n|1}}{P_{(n1)|1}}}C\text{ }\text{ as }n\mathrm{},(C1)`$ (9)   
where $`C`$ is a constant. $`C`$ indicates the asymptotic behavior of $`P(n)`$ as $`n\mathrm{}`$. If $`C<1`$, the probability distribution is asymptotically exponential, while for $`C=1`$, the probability distribution is a power law with exponent $`3/2`$.

Let us now examine the behavior of $`P(n)`$ when $`n10^4`$, the more relevant case in the examples to follow. Using Eqs. (4)-(7), we can numerically calculate $`P(n)`$ for different sets of $`p_i`$. We define $`m`$ as the expected number of daughters per organism, given a set of probabilities $`p_i`$;

$`m={\displaystyle \underset{i}{}}ip_i.`$ (10)   
We see that the branching rate $`m`$ (the control parameter) is a good indicator of the shape of the probability curve (Fig. 1). When $`m`$ is close to $`1`$, the distribution is nearly a power law, and the further $`m`$ diverges from $`1`$, the further the curve diverges from a power law towards an exponential. When $`m=1/2`$, the curve is completely exponential. For a population of organisms, $`m`$ is a measure of the tendency for new generations to grow, or shrink, in number. A value of $`m>1`$ indicates a growing generation size, which implies that there will, on average, be no generation with no daughters, and that the expected number of total organisms is infinite. Conversely, $`m<1`$ indicates a shrinking population size: There will be a final generation with no daughters, and the expected number of organisms is finite. When $`m=1`$, the system is in between the two regimes (the system is said to be “critical”), and only then is a power law distribution found. In general, the branching rate is determined by the ratio of the rate of introduction of competitors $`R_c`$ to the intrinsic rate of growth of existing assemblages $`R_p`$ via

$`m=(1+{\displaystyle \frac{R_c}{R_p}})^1,`$ (11)   
as can be shown by assuming stationarity. As this ratio goes to $`0`$, $`m1`$ and the system becomes critical.

In the following section, we explore systems where the “organisms” are individual members of species or taxa in a taxonomic tree, and $`m`$ is the average number of exact copies an individual makes of itself, or the average number of new taxa of the same supertaxon a taxon spawns, respectively. The same thinking can be applied to tumbling sites in a sandpile model, where $`m`$ would stand for the average number of new tumbles directly caused by a tumbling site.

## III Applications

### A Neutral Model

We first present a simple simulation to test our analysis and lay the groundwork for the exploration of more complicated systems. Consider a population of organisms on a finite two-dimensional Euclidean lattice, one organism to a grid square. Each organism can be viable or sterile. All viable organisms replicate approximately every $`\tau `$ time steps (there is a small random component to each individual’s replication time to avoid synchronization effects), while sterile organisms do not replicate. When an organism replicates, its daughter replaces the oldest organism in the parent’s $`9`$-site neighborhood (Fig. 2). We define the fidelity $`F`$ as the probability that the organism will create a daughter of the same type as itself and the corresponding genomic mutation rate $`R`$ ($`=1F`$) at which it creates copies different from itself. The genomic mutation rate is actually the sum of two rates, a probability $`R_n`$ for the daughter to be viable but to be of a new genotype, different from that of the parent (neutrality rate), and a probability $`R_s`$ of the daughter being sterile. Of course, $`R_n+R_s=R`$. Note that all viable mutant daughters still share the same replication time $`\tau `$—all mutations are neutral (See Fig. 3). Such a system gives rise to abundance distributions of power law and near-power law type, which can be predicted as follows.

The total number of organisms is determined by the size of the grid. We write equilibrium conditions for the total number of organisms $`\rho _A`$, and for the total number of viable organisms $`\rho _V`$,

$`\mathrm{\Delta }\rho _A`$ $``$ $`a\rho _V\rho _A=0,`$ (12)   
$`\mathrm{\Delta }\rho _V`$ $``$ $`v\rho _V\rho _V=0,`$ (13)   
where $`a`$ is the average number of daughters (viable and sterile) a viable organism has, and $`v`$ is the average number of viable daughters a viable organism has. Introducing $`m`$—the average number of true daughters (daughters which share the parent’s genotype) for a viable organism—we see that

$`v={\displaystyle \frac{F+R_n}{F}}m=(F+R_n)a.`$ (14)   
From Eqs. (12)-(14), we obtain steady state solutions for $`a`$ and $`m`$,

$`a`$ $`=`$ $`{\displaystyle \frac{F^1}{1+\frac{R_n}{F}}},`$ (15)   
$`m`$ $`=`$ $`{\displaystyle \frac{1}{1+\frac{R_n}{F}}}.`$ (16)   
Using the branching process model, we can predict the abundance curve from the values of $`a`$ and $`m`$ (or conversely, $`F`$ and $`R_n`$). Fig. 4 shows abundance data for two neutral model runs with differing values of $`R_n`$ (and consequently $`m`$), along with predicted distributions (which use only $`R_n`$ and $`F`$ as parameters) based on the branching model. Although the distribution patterns are very different, both are fit extremely well by the branching process model’s predicted curves. In Eq. (16), note that $`R_n`$ is the rate of influx of new genotypes (and therefore new competitors for space), while $`F`$ is the rate of growth of existing genotypes. The value of $`m`$ is determined by the ratio of these two rates. Unless the total number of creatures is increasing, $`m1`$ ($`m=1`$ if and only if $`R_n0`$ and new competing genotypes are introduced at a vanishing rate).

### B Artificial Life

Our next system is the artificial life system sanda , an example of environments which host digital organisms . In this system, while the organisms occupy a two-dimensional grid as in the neutral model detailed above, the organisms are no longer simple, and instead each has a complex genotype consisting of a string of assembly language-like instructions (Fig. 5). Each organism independently executes the instructions of its genotype, and this genotype determines the organism’s replication time $`\tau `$. Unlike the neutral model, the system allows non-neutral mutations which lead to new genotypes with both lower and higher replication times than the parent.

The system and the instructions are designed so that the organisms can self-replicate by executing certain sequences of instructions. The replication time of an organism is not a predetermined constant, rather it is determined by the genotype of the organism: Organisms can replicate faster or slower than other competing organisms with different genotypes. For an organism to successfully replicate, its genotype must contain information which allows the organism to allocate temporary space (memory) for its daughter, replicate its genotype (one instruction at a time) into this temporary space, and then to divide, placing its daughter in a grid site of its own (Fig. 5). As in the neutral model, on division, the daughter replaces the oldest organism in its parent’s $`9`$-site neighborhood.

Organisms, depending on their genotype, may not be able to replicate (may be sterile) or may only be able to replicate imperfectly (resulting in no true daughters). Also, the copy instruction, which the organisms must use to copy instructions from their own code into that of their nascent daughters, has a probability of failure (copy mutation rate), which can be set by the experimenter. When the copy instruction fails, an instruction is randomly chosen from all the instructions available to the organisms (the instruction set) and written in the string location copied to. Copy mutations also lead to non-true daughters. The instruction set is robust; copy errors (mutations) induced during the replication of viable organisms have a non-vanishing probability of creating viable new organisms and genotypes. Indeed, by selecting for certain traits (such as the ability to perform binary logical operations) by increasing the relative speed at which instructions are executed in organisms which carry these traits, the system can be forced to evolve and find novel genotypes which contain more information (and less entropy) than their ancestors. Even without this external selection, the system evolves organisms (and genotypes) which replicate more efficiently in less executed instructions.

As a result of this evolution, the fidelity and neutral mutation rate are not fixed, but can vary with the length of an organism’s genome and the instructions contained therein. Also, new genotypes formed by beneficial mutations that allow faster replication than previously existing genotypes will have (on average) an increasing number of organisms—$`m>1`$—until the new, faster replicating genotypes fill up a sizable portion of the grid. All these factors combine to make predicting the abundance distributions for sanda much harder than for the neutral model.

Indeed, rather than being constant during the course of a sanda experiment, $`R_n`$ and $`F`$ will vary unpredictably as the population of organisms occupies different areas in genotypic phase space. Certain genotypes may be brittle, allowing very few mutations that result in new viable genotypes. The length of the organisms may change, changing both the genomic mutation rate and the neutrality rate. Genotypes exist which make systematic errors when copying, which decreases the fidelity. In short, the dynamics of these digital organisms are complex and messy, much like those of their biochemical brethren. These variations are observed at the same time across different organisms in the population, and are also observed with the progression of time. Still, we attempt to predict the abundance distributions by approximating the ratio of neutral mutations to true copies by the observed ratio of viable genotypes to total number of viable organisms ever created:

$`{\displaystyle \frac{R_n}{F}}{\displaystyle \frac{N_g}{N_v}},`$ (17)   
where $`N_g`$ is the total number of viable genotypes observed during a sanda run and $`N_v`$ is the total number of viable organisms. This relation should hold approximately under equilibrium conditions. Then, Eq. (16) becomes

$`m(1+{\displaystyle \frac{N_g}{N_v}})^1,`$ (18)   
and from Eq. (15)

$`a={\displaystyle \frac{m}{F}}.`$ (19)   
The fidelity $`F`$ is inferred from the average length $`l`$ of genotypes during a run and the (externally enforced) per-instruction copy mutation rate $`\gamma `$, $`F=(1\gamma )^l`$. Because we estimate $`m`$ and $`a`$ from macroscopic observables averaged over the length of a run, we expect some error in our results due to the shifting dynamics of the evolution of genotypes as the system moves in genotypic phase space.

The abundance data from two different sanda runs are shown in Fig. 6 with the predicted abundance curves. The two runs shared the same grid size and per-instruction copy mutation rate, and were started with the same initial genotypes, but the runs evolved into different regions of genotypic phase space and consequently had significantly different statistics. Considering the many gross approximations we have made, the agreement between our prediction and the experimental data is surprisingly good (especially as no fitting is involved). Sanda is most closely related to an asexually replicating biological population, such as colonies of certain types of bacteria occupying a single niche. The genotype abundance distributions measured in sanda are analogous to the species or subspecies abundance distributions of its biological counterparts. In general, species abundance distributions are complicated by the effects of sexual reproduction, and of the localized and variable influences of other species and the environment on species abundances. However, we believe the branching model—used judiciously—can be helpful in the study of such distributions as well.

### C Evolution

Rank-abundance distributions at taxonomic levels higher than species (e.g., the distribution of the number of families per order) are simpler to model than species abundance distributions, as the effects of the complications noted above are weak or nonexistent. We find that the available data is well fit by assuming no direct interaction or fitness difference between taxa . The shapes of rank-frequency distributions of taxonomic and evolutionary assemblages found in nature are surprisingly uniform. Indeed, Burlando has speculated that all higher-order taxonomic rank-frequency distributions follow power laws stemming from underlying fractal dynamics . We believe this conclusion is hasty: The divergence of the distributions from power law can be observed by applying appropriate binning methods to the data. (See Appendix.) Yule  attempted a branching process model explanation of these distributions, and claimed that divergence from power law of rank-abundance patterns was transient and indicated a finite time since the creation of the evolutionary assemblage. Our model indicates that this is not generally the case. We find that the divergence from power law is not a result of disequilibration, but is an inherent property of the evolutionary assemblage under consideration and that this divergence provides insight into microscopic properties of the assemblage (e.g., the rate of innovation).

Say, for example, that we are interested in the rank-frequency distribution of the number of families in each order for fossil marine animal orders. We assume that all new families and orders in this assemblage originate from mutations in extant families. Then, we can define rates of successful mutation $`R_f`$ for mutations which create new families in the same order as the original family, and $`R_o`$ for mutations which create an entirely new order. In this case, unlike the cases treated above, we approximate $`a\mathrm{}`$; many individual births and mutations occur, but the proportion that are family- or order-forming is minuscule. Finally, assuming a quasi-steady state (the total numbers of orders and families vary slowly ), we rewrite Eq. (16),

$`m`$ $``$ $`(1+{\displaystyle \frac{R_o}{R_f}})^1`$ (20)   
$``$ $`(1+{\displaystyle \frac{N_o}{N_f}})^1,`$ (21)   
in terms of $`N_o`$ and $`N_f`$, the total numbers of orders and families respectively. As in the previous systems studied, $`R_o`$ is the rate of creation of new—competing—orders, while $`R_f`$ is the rate of growth of existing orders, and $`m`$ is determined by their ratio.

Data for the abundance distribution of number of families in fossil marine animal orders  are shown in Fig. 7. We obtained values for $`N_o`$ and $`N_f`$ directly from the fossil data to generate the predicted curve with no free parameters. The agreement is very good, much better than that for the sanda runs where evolutionary parameters such as the fidelity $`F`$ and the neutrality $`R_n`$ were constantly changing. Comparing $`m`$ and the resultant abundance curves with those obtained above for the rank-abundance distribution of sanda genotypes leads us to the expected conclusion that the probability of creation of a new genotype in sanda per birth is much higher than the probability of a new family creating an order in natural evolution. Indeed, a wide variety of higher-order taxonomic assemblages have abundance distributions consistent with $`m`$ near 1 . We believe this is a robust result of the evolutionary process. Low values of $`m`$ may not be observed for large taxon assemblages for several reasons. A small value of $`m`$ implies either a small number of individuals in the assemblage, or a very specialized niche with a very low rate of taxon formation. A low number of individuals would lead to a low probability of the taxon being discovered and cataloged by biologists. A small number of individuals and taxa would result in an assemblage with too few taxa to give us a clear statistical picture. Also, since such an assemblage would have a small population, be incapable of further adaptation, or both, we expect it would be more susceptible to competition and environmental effects leading to early extinction.

### D Sandpile Models

It was originally suggested that the self-organization observed in the sandpile model of Bak, Tang and Wiesenfeld (BTW)  (and the power laws it displayed) was an inherent property of the system, while it now seems established that the system is actually tuned by waiting until avalanches are over before dropping new grains—this is equivalent to allowing non-local interactions . The same conclusion is reached when using a branching process to describe the avalanche dynamics. Branching processes have been applied to sandpile models as early as 1988  (see also, ). Using a mean-field approach in higher dimensions ($`d4`$), power law distributions for the size of avalanches $`s(n)`$ can be obtained analytically, and critical exponents can be calculated exactly to reveal $`s(n)n^{3/2}`$  in the limit of infinitesimally small driving. This is supported by numerical simulations. However, for lower dimensions, sandpiles will “interfere” with themselves, and a smaller exponent is found. Attempts to calculate the effects of this “final-state” interaction through renormalization have as yet not been completely successful . Still, the phenomenon of “violations” of power-law behavior due to $`m<1`$ (non-critical branching process) can be seen there as well.

## IV Discussion

The Galton-Watson branching process generates power law distributions when its control parameter $`m=1`$. In all the systems we have examined above,

$`m=(1+{\displaystyle \frac{R_c}{R_p}})^1`$ (22)   
is determined by the ratio of the rate of introduction of competitors $`R_c`$ to the intrinsic rate of growth of existing assemblages $`R_p`$. As this ratio goes to $`0`$, $`m1`$ and the system becomes critical.

This relation can be translated into the standard relation between an order parameter

$`\alpha ={\displaystyle \frac{R_c}{R_p}}`$ (23)   
and a new form for the control parameter

$`\mu =m^1.`$ (24)   
Writing $`\alpha `$ in terms of $`\mu `$,

$`\alpha =\{\begin{array}{cc}(\mu \mu _c)^\beta & (\mu >\mu _c),\hfill \\ 0& (\mu \mu _c),\hfill \end{array}`$

where $`\mu _c=1`$ and $`\beta =1`$ (Fig. 8). The order parameter represents the rate at which competition is introduced in the system (the strength of selection). A value of the control parameter $`\mu <\mu _c`$ implies a system with no competition and no selection—an exponentially growing population. Values of $`\mu `$ higher than $`\mu _c`$ indicate that new competition is always being introduced and that all existing species or avalanches must eventually die out. When $`\mu =\mu _c`$, competition is introduced at a vanishingly small rate, and we find the critical situation where separation of scales occurs.

For sandpile models, this $`\alpha `$ is arbitrarily set close to $`0`$ by using large lattice sizes (reducing dissipation) and waiting for avalanches to finish before introducing new perturbations (resulting in an infinitesimal driving rate and a diverging diffusion coefficient). For the biological and biologically-inspired systems we have considered, the control parameter is not set arbitrarily to a critical value. However, the dynamics of the evolutionary process, in which it is much harder to effect large jumps in fitness and function than it is to effect small ones, lead to naturally observed values of $`\alpha `$ being small, especially for higher taxonomic orders. The dynamics of evolution act, robustly, to keep $`\mu `$ near $`\mu _c`$. This in turn leads to a near power law pattern for rank-frequency distributions.

We have shown that the apparent power laws of avalanches in species-abundance distributions in artificial life systems, as well as rank-abundance distributions in taxonomy can be explained by modeling the dynamics of the underlying system with a simple branching process. This branching process model successfully predicts, with no free parameters, the observed abundance distributions—including their divergence from power law.

A branching process approach may allow the deduction of the microscopic parameters of the system directly from the macroscopic abundance distribution. We find that we can identify a control parameter—the average number of new events an event directly spawns, and an order parameter— the rate of introduction of competing events into the system, and that these are related in a form familiar from second order phase transitions in statistical physics.

###### Acknowledgements.

We are grateful to the late Prof. J. J. Sepkoski for kindly sending us his amended data set for fossil marine animal families, as well as an anonymous referee for insightful comments. J. C. thanks M. C. Cross for continued support and discussions. Access to the Intel Paragon XP/S was provided by the Center of Advanced Computing Research at the California Institute of Technology. This research was supported by the NSF under contract Nos. PHY-9723972 and DEB-9981397. Part of this work was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration.    

## A Binning Methods

When dealing with event distributions best plotted on single log or double log scales (such as exponential and power law distributions), care must be taken in the proper binning of the experimental data. Say we are interested in the probability distribution $`P(n)`$ of an event distribution over positive integer values of $`n`$. We conduct $`N`$ trials, resulting in a data set $`Q(n)`$ of number of events observed for every $`n`$ value. For ranges of $`n`$ where the expected or observed number of events for each $`n`$ is much higher than 1, normally no binning is required. However, for ranges of $`n`$ where $`Q(n)`$ or $`P(n)`$ is small, binning is necessary to produce both statistically significant data points, and intuitively correct graphical representations. A constant bin size has several drawbacks: One must guess and choose an intermediate bin size to serve across a broad range of parameter space, and the shape and slopes of the curve (even in a double log plot) are distorted . These disadvantages can be overcome by using a variable bin size. However, choosing bin sizes for variable binning is time-consuming and arbitrary—different choices will lead to different conclusions. We propose two related methods of systematically determining appropriate variable bin sizes. Both methods lead to binned data which help in visualizing the underlying distribution (slopes and shapes are conserved).

For the first method (the Data Threshold Method), we start by selecting a threshold value $`T`$. Starting from $`n=1`$ and proceeding to higher values, no binning is done until a value of $`n`$ is found for which $`Q(n)<T`$. When such a value $`n_s`$ is found, subsequent $`Q(n)`$ values are added to this amount until the sum of these values is greater than the threshold value,

$`{\displaystyle \underset{n=n_s}{\overset{n_l}{}}}Q(n)>T.`$ (A1)   
We then have a bin size $`(n_ln_s+1)`$, with value $`_{n=n_s}^{n_l}Q(n)`$. When plotting, it is convenient to plot this as a single point at the midpoint of $`[n_s,n_l]`$, with an averaged value,

$`({\displaystyle \frac{n_s+n_l}{2}},{\displaystyle \frac{_{n=n_s}^{n_l}Q(n)}{n_ln_s+1}}).`$ (A2)   
This yields a graphical representation with little distortion and good predictive power (Fig. 9). This binning procedure is continued until no more data remains to be binned.

The second binning method (the Template Threshold Method), uses a predicted probability distribution $`P(n)`$, or a reasonable surrogate. Again, we define a threshold value for fitting $`T`$. However, in this case, the bin sizes are determined by comparing values of the expected distribution

$`E(n)=P(n)\times N`$ (A3)   
to $`T`$. Starting from $`n=1`$ and proceeding to higher values, no binning is done until a value of $`n`$ is found for which $`E(n)<T`$. When such a value $`n_s`$ is found, subsequent $`E(n)`$ values are added to this amount until the sum of these values is greater than the threshold value,

$`{\displaystyle \underset{n=n_s}{\overset{n_l}{}}}E(n)>T.`$ (A4)   
We then have a bin of $`[n_s,n_l]`$ with corresponding size $`(n_ln_s+1)`$. The average value associated with this bin is

$`{\displaystyle \frac{_{n=n_s}^{n_l}Q(n)}{n_ln_s+1}}.`$ (A5)   
This procedure is repeated until the data is exhausted. For this method, the data may be graphically represented either as a single point per bin (as in the data threshold method above), or as a point (showing the associated average value) for each measured (non-zero) data point $`Q(n)`$.

The data threshold method requires no a priori knowledge, and is a good predictor of the underlying distribution. However, when there are few data points, the template threshold method is more reliable. For both methods, a range of $`T`$ should be tried and the best $`T`$ (neither over- or under-binning) chosen.