# How to Count Fish Using Mathematics

"How many fish are in the sea?" Fields Medalist Charlie Fefferman asked at the end of my talk about a coagulation-fragmentation model of animal group sizes during a workshop held in Princeton, N.J., in September 2015. Fefferman’s naïve-sounding question cut to the heart of the resource estimation issues that the model was designed to address. My work on this model, with collaborators Pierre Degond (Imperial College London) and Jian-Guo Liu (Duke University) [1], began with our interest in Japanese fisheries scientist Hiro-Sato Niwa’s investigation of the frequency distribution of school sizes for pelagic fish – fish that roam in the mid-ocean.

Niwa analyzed a large amount of observational data and found that, when scaled by an expected school size, it collapses well onto a single, highly non-Gaussian curve. Niwa worked under the hypothesis that simple rules for random merging and splitting of schools could explain the observed size distribution. He modeled this situation with a stochastic differential equation and performed merging-splitting simulations to estimate variance. Impressively, he then used Itô calculus to find an explicit solution of the model, in the form

\[f(s) = s^{-1} \exp(-s + \tfrac12 se^{-s}). \tag1 \]

**Figure 1.**Empirical school-size distribution of six types of pelagic fish. Data types and sources are listed in Table 1 of [3]. Data is scaled by empirical average. The solid line corresponds to (1). Figure reprinted from [3], with permission from Elsevier.

\[\int_0^\infty s f(s)\,ds = \int_0^\infty s^2 f(s)\,ds.\]

Niwa mentions that he checked this numerically to an accuracy of 128 digits! (Proving this with elementary calculus is quite a nice challenge.)

Our own work [1] involves solving and analyzing a more faithful model of Niwa’s simulated merging-splitting process. He was very aware of the type of coagulation-fragmentation model that we study, but its mathematical treatment turns out to require new developments in complex function theory related to Laplace transforms. One result of our analysis is that the coagulation-fragmentation equation has an exact solution that takes the form

\[f_\star(s) = g(s) \exp(-\tfrac89 s)\,,\]

where \(g(s)\) is a *completely monotone* function—a smooth positive function that is decreasing, convex, and has derivatives that alternate in sign indefinitely—that behaves as

\[g(s) \sim c_0 s^{-2/3} \ \ \textrm{as} \enspace s\to 0, \qquad g(s) \sim c_\infty s^{-3/2} \ \ \textrm{as} \enspace s \to \infty.\]

Moreover, the profile \(f_\star\) is computable in terms of a series in powers of \(s^{1/3}\), namely

\[f_\star(s) = (6s)^{-2/3} \sum_{n=0}^\infty \frac{2(-1)^n}{\Gamma(\frac43-\frac23n)}(6s)^{n/3}.\]

Niwa’s prefactor \(s^{-1}\) has exponent \(−1\) that interpolates nicely between the \(-\frac23\) and \(-\frac32\) that we find valid in the small and large limits, respectively. Furthermore, the distribution \(f_\star\) differs minimally from Niwa’s \(f\) (less than 20%, hardly noticeable on the scale of Figure 1) in the shoulder region of the log-log plot. The crossover in power-law behavior seems likely to complicate the fitting of exponentially-truncated power laws to empirical data on animal group size, as has been attempted by a number of researchers since the mid-1990s.

The coagulation-fragmentation model is essentially a simple kind of Boltzmann equation from the kinetic theory of gases. Instead of modeling collisions and the scattering of gas molecules, it models the merging and splitting of clusters. Prominent Polish physicist Marian Smoluchowski first derived coagulation equations a century ago. Scientists have since used them to model cluster-size distributions in a wide variety of physical systems, from soot and smog particles in aerosol physics to planetesimals and dark-matter galactic halos in astrophysics. In the case of pelagic fish, Niwa argued that the collision and splitting rates would not permit detailed balance, and that no analog of Boltzmann’s H-theorem could explain equilibrium. Consequently, studying the dynamics posed a novel problem.

Degond, Liu, and I made progress by rewriting the equations using what we now call *Bernstein transforms*. Though these are “just” antiderivatives of familiar Laplace transforms, they have a number of distinctive properties that make them a worthy subject of a recent book by René Schilling, Renming Song, and Zoran Vondracek [5]. For example, the composition of two Bernstein transforms is a Bernstein transform, and so is the inverse to a Bernstein transform’s antiderivative.

In the case of the Niwa-motivated model in the continuum limit of large populations, the Bernstein transform of the size distribution, given by

\[\varphi(q,t) = \int_0^\infty (1-e^{-qs}) f(s,t) ds,\]

turns out to satisfy an attractive integro-partial differential equation, namely

\[\frac{\partial \varphi}{\partial t} = -\varphi^2 -\varphi + \frac2q\int_0^q \varphi(r,t)\,dr. \tag2 \]

The steady states of this equation have a scale-invariant shape \(\varphi_\star\), determined implicitly by the simple algebraic relation

\[\varphi_\star = q(1-\varphi_\star)^3. \tag3 \]

To recover information about the equilibrium size distribution \(f_\star\) from its Bernstein transform, we use a curiously strong theorem in Bernstein function theory involving global analyticity properties. This theorem states that \(f_\star\) itself is a Laplace transform (and thus completely monotone) if and only if its own Bernstein transform \(\varphi_\star\) is a Pick function – a globally-complex analytic function on the upper half plane \(\mathcal{H}\), mapping \(\mathcal{H}\) into \(\mathcal{H}\). This turns out to be true for the \(\varphi_\star\) given by \((3)\).

We can also prove that \(f_\star\) is a global attractor for the model’s dynamics. This is derived from \((2)\) by improving the classical continuity theorem for Laplace transforms to demonstrate that pointwise convergence of Bernstein transforms corresponds to weak convergence of corresponding measures on the compactified half-line \([0,\infty]\).

A further twist in the story pertains to discrete size distributions, to which Niwa’s merging-splitting model leads for finite total population. We were unable to solve this model with Niwa’s original postulated splitting rates, namely equal rates for a group of size s splitting into two groups of size \(j\) and \(s-j\) for \(j=1\) to \(s-1\). However, we noticed that if we modify the model slightly to include the trivial cases \(j=0\) and \(j=s\), the Bernstein transform again solves the same integro-differential equation \((2)\) after a change of variables!

The curiously strong Bernstein theorem doesn’t work to describe equilibrium in the discrete-size case, but we discovered a discrete analog. This analog explains which sequences \((c_j)_{j=0}^\infty\) are sequences of moments of measures on positive intervals \([0,T]\). The characterization is in terms of simple global analyticity properties of the generating function

\[F(z) = \sum_{j=0}^\infty c_jz^j.\]

For example, one characterization says \(z F(z)\) is a Pick function analytic on \((-\infty,T^{-1})\). Through generating functions, discrete equilibria in Niwa’s model are related to the classical *Fuss-Catalan numbers* of combinatorics, given by

\[A_n(p,r) = \frac{r}{pn+r}\binom{pn+r}{n},\]

as well as sequences of moments of random matrix ensembles, infinitely divisible probability distributions on the natural numbers, and convolution semigroups of completely monotone sequences. The connections are laid out in [2]. An unusual detail, related to the fact that the map \(r\mapsto (A_n(p,r))_{n\ge0}\) is a convolution semigroup, is the formula

\[\sum_{n=0}^\infty \binom{pn}n z^n =\int_0^1\frac1{1-f_p(u)z}du,\]

where

\[f_p(u) = \frac{\sin^p\pi u}{\sin{(\frac1p \pi u)}\sin^{p-1}(\frac{p-1}p \pi u)}.\]

This formula was discovered due to Tewodros Amdeberhan’s use of the amazing *Zeilberger algorithm* [4], which derives and proves combinatorial formulas, and a striking integral expression for binomial coefficients, namely

\[\binom nk = \frac1\pi \int_0^\pi \frac{ \sin^n x} {\sin^k(\frac kn x) \sin^{n-k}( \frac{n-k}n x) }\,dx.\]

I would like to know a simple proof, but please don’t assign this as a calculus problem! It follows from results mentioned in [6].

Returning to Fefferman’s question about the number of fish in the sea, I didn’t have such a number handy. However, I could mention something about data collection methods, such as the use of sonar-measured thickness as a proxy for school size. Pelagic fish schools can extend over several kilometers in diameter and be tens of meters thick.

A central conclusion of Niwa’s modeling and data analysis is that the fish school size distribution is highly non-Gaussian. It is not characterized by a normal distribution about a mean, but instead has exponential-like tails that are much ‘fatter’ than Gaussian ones. Thus, observations of large schools are likely to be much more frequent than Gaussian statistics suggest. This indicates that the use of familiar Gaussian models may easily overestimate the total population.

Our work provides a mathematically consistent foundation for this line of thought. One hopes that the scientists involved in resource estimation are aware of the unsuitability of Gaussian models for this purpose, and that the management of ocean fisheries is not based on inaccurate models and estimates.

**References**

[1] Degond, P., Liu, J.-G., & Pego, R.L. (2017). Coagulation-fragmentation model for animal group-size statistics. *J. Nonlinear Sci., 27*(2), 379-424.

[2] Liu, J.-G., & Pego, R.L. (2016). On generating functions of Hausdorff moment sequences. *Trans. Amer. Math. Soc., 368*, 8499-8518.

[3] Niwa, H.S. (2003). Power-law versus exponential distributions of animal group sizes. *J. Theo. Biol., 224*, 451-457.

[4] Petkovšek, M., Wilf, H.S., & Zeilberger, D. (1996). *A = B*. Wellesley, MA: AK Peters, Ltd.

[5] Schilling, R.L., Song, R., & Vondracek, Z. (2010). *Bernstein Functions*. Vol. 37 of *De Gruyter Studies in Mathematics*. Berlin, Germany: Walter de Gruyter & Co.

[6] Simon, T. (2014). Comparing Fréchet and positive stable laws. *Electron. J. Probab., 19*(16), 1-25.