# Point Processes and the Statistics of Cellular Neighborhoods
in Simple Multicellular Organisms

Anand Srinivasan<EMAIL_ADDRESS>Department of Applied Mathematics and
Theoretical Physics, Centre for Mathematical Sciences,
University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, United Kingdom
Steph S.M.H. Höhn<EMAIL_ADDRESS>Department of Applied Mathematics and
Theoretical Physics, Centre for Mathematical Sciences,
University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, United Kingdom
Raymond E. Goldstein<EMAIL_ADDRESS>Department of Applied
Mathematics and Theoretical Physics, Centre for Mathematical Sciences,
University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, United Kingdom

###### Abstract

Recent work on several distinct multicellular organisms has revealed a
hitherto unknown type of biological noise; rather than a highly regular
arrangement, cellular neighborhood volumes, obtained by Voronoi tessellations
of the cell locations, are broadly distributed and consistent with gamma
distributions. We propose an explanation for those observations in the case of
the alga Volvox, whose somatic cells are embedded in an extracellular matrix
(ECM) they export. We show that the simplest models for stochastic ECM
generation in one and two dimensions are point processes whose Voronoi
tesselations are demonstrably governed by gamma distributions. These results
highlight a universal consequence of intrinsic biological noise on the
architecture of certain tissues.

Some of the simplest multicellular organisms have an architecture consisting
of tens, hundreds, or thousands of cells arranged within an extracellular
matrix (ECM), a network of proteins and biopolymers secreted by the cells.
They often adopt regular geometric forms: the linear chains and rosettes of
choanoflagellates [1, 2], the sheets or spheres of green algae [3], and the
cylinders of sponges [4]. While the arrangement of cells within the ECM
appears regular, recent work [5] has revealed a hitherto undocumented form of
disorder found by assigning neighborhoods to the cells through a Voronoi
tessellation based on the cell centers. Strikingly, both the lab-evolved
organism “snowflake yeast” [6] (a ramified multicellular form found by
repeated rounds of selection for sedimentation speed) and the green alga
Volvox carteri display broad distributions of Voronoi volumes accurately fit
by gamma distributions. These observations are central to an understanding of
a general question in developmental biology: How do cells produce structures
external to themselves in an accurate and robust manner?

Volvox (Fig. 1(a)) is one of the simplest systems to understand the
statistical fluctuations in the generation of an ECM. The adult consists of
$\sim\\!\\!10^{3}$ somatic cells embedded at the surface of a transparent ECM,
a thin shell $\sim\\!\\!500\,\mu$m in diameter and $\sim\\!\\!30\,\mu$m thick;
the organism is $\sim 90\%$ ECM. Daughter colonies develop from germ cells at
the inner ECM surface; repeated rounds of binary division produce a raft of
cells held together by cytoplasmic bridges that remain after incomplete
cytokinesis. Following “embryonic inversion” that turns the raft inside-out
[7], daughter colonies enlarge (at fixed somatic cell number) by generation of
the ECM via export of proteins from those cells, expanding the colony radius
to the final size over the course of a day, during which the widely-
distributed neighborhood volumes appear.

Figure 1: The green alga Volvox and one-dimensional models for cell positions.
(a) Adult V. carteri with cell types labelled. Inset shows section of Voronoi
tessellation around somatic cells (visible as dots centered in each Voronoi
polygon). (b) Voronoi cell areas $v$ approximately follow the translated gamma
distribution (1) [5], computed by a maximum-likelihood fit. (c) The Poisson
process on $[0,\infty)$, where cell locations are indicated by vertical line
segments and intercalating ECM is colored. Conditioning on the cell count or
length leads to two different realizations of circular processes, IS.1 and
IS.2, respectively, whose numerically obtained empirical distribution are
shown in (d) and (e). Dashed and solid lines indicate approximate or exact
density functions, respectively. (f) Empirical distribution for model IS.4,
with no cell overlaps, compared to analytial result.

For Volvox, the general issue above becomes the question of how ECM is
generated such that the spheroidal form is maintained during the dramatic
enlargement of the colony and yet the cellular neighborhood volumes are
broadly distributed. A biological answer might invoke cell-cell signaling in
response to mechanical forces as a mechanism to coordinate growth and would
ascribe the distribution of neighborhood volumes to imperfections in that
process. Surprisingly, the novel problem of cellular neighborhood
distributions is so little-studied that we do not even understand
quantitatively the feedback-free case, which serves as a benchmark for any
analysis of the role of correlations. We note that in granular physics, gamma
distributions have been argued to arise from maximizing entropy when space is
divided up into a set of subvolumes constrained to add up to a fixed total
volume [8]. But this then begs the question of why biological systems should
follow a maximum-entropy principle.

Here, we examine a range of idealized models for distributions of cells within
a thin growing ECM. Regardless of the process that generates ECM, the
assumption that cell counts within any two disjoint regions are independent
random variables—along with three technical and physically reasonable
assumptions about the vanishing probability of realizing particular
configurations [9]—is the defining feature of a Poisson point process [11]. In
static and dynamic one-dimensional models, where cells sit in a circular ECM,
we show analytically that Voronoi segments obey gamma distributions. In two
dimensions we describe strong numerical evidence for consistency with gamma-
type distributions and conjecture an analytical form. The following is a non-
technical summary, with details provided in the Supplementary Material [9].

Our analysis is motivated by light-sheet measurements of somatic cell
positions in adult V. carteri [5]. Figure 1(a) shows a darkfield image of V.
carteri and a section of the Voronoi tessellation around the cells. The area
distribution of Voronoi partitions across $6$ organisms is shown in Fig. 1(b)
along with a maximum-likelihood fit to the affine-transformed $k$-gamma
distribution,

$p(v)=\frac{1}{\bar{v}-v_{c}}\frac{k^{k}x^{k-1}}{\Gamma(k)}\exp(-kx),\quad
x=\frac{v-v_{c}}{\bar{v}-v_{c}},$ (1)

where $\bar{v}$ is the mean and $v_{c}$ is the somatic cell area, the lower
bound on $v$ due to the non-overlapping of cells. We wish to explain the
origin of the gamma distribution itself and its validity in the form
translated by $v_{c}$.

Notation. In the following, random variables are denoted by capital letters
$X_{i;\alpha,\beta,\ldots}$ accompanied by indices $i$ and parameters
$\alpha,\beta,\cdots$. Variables $i\neq j$ are independent unless otherwise
noted. The letters $X_{i;\lambda}$, $Y_{i;k,\lambda},Z_{i;\alpha,\beta}$, and
$N_{i;\lambda}$ are reserved for exponential, gamma, beta, and Poisson random
variables respectively. $X\overset{\text{pdf}}{\ \sim\ }f_{X}(x)$ indicates
that $X$ has the probability density function $f_{X}(x)$, with
$X\overset{\text{cdf}}{\ \sim\ }F_{X}(x)$ indicating the same for cumulative
density. $X_{i}\overset{\text{d}}{\to}Y$ indicates $X_{i}$ converges to $Y$ in
distribution. Models are named by the ECM dimension $d$ (I or II) and the
property of being static (S) or dynamic (D), as summarized in Table 1.

Table 1: Summary of models Model | Description
---|---
IS.0 | fixed N on the half-line
IS.1 | fixed N, variable C on the circle
IS.2 | fixed N, fixed C on the circle
IS.3 | variable N, fixed C on the circle
IS.4 | fixed N and C with minimum spacing on the circle
ID.1 | Brownian motion on the circle
ID.2 | noncolliding Brownian motion on the unit circle
ID.3 | noncolliding Brownian motion on a growing circle
IIS.1 | variable N on the periodic unit square
IIS.2 | variable N with minimum spacing on the unit sphere
IID.1 | fixed N, Brownian motion on the unit sphere

One-dimensional equilibrium models. We first consider static models in which
cells sit within a one-dimensional ECM, and view the arclength of ECM
generated by each cell as a random variable. As in the Gibbs, microcanonical,
and grand canonical ensembles, we consider three types of random
configurations of cells in the ECM: (i) fixed cell counts $N=n$, with random
intercell spacings and circumferences; (ii) fixed circumferences $C=c$ and
counts $N=n$, with random positions $\\{R_{i}\\}_{i=1}^{N}$; and (iii) fixed
circumferences $C=c$, with random cell counts $N$ and positions. Like the
convergence of the ensembles in the thermodynamic limit, we show that the same
gamma distribution arises in the large-$n,c$ limits of these cases.

Model IS.0 (on the half-line) —Consider a semi-infinite Volvox modelled as a
homogeneous Poisson process of rate $\lambda$ on the half-line $[0,\infty)$,
as in Fig. 1(c). The sorted cell positions $R_{1},R_{2},R_{3},\cdots$ have
exponentially distributed spacings $X_{i;\lambda}$. The Voronoi lengths
$V_{i}=(X_{i;\lambda}+X_{i+1;\lambda})/2$ are proportional to the sum of
exponentially-distributed random variables with the same rate, which, by a
convolution [9], are gamma-distributed with $k=2$,

$V_{i}=\frac{1}{2}Y_{;2,\lambda}\overset{\text{pdf}}{\ \sim\
}4\lambda^{2}ve^{-2\lambda v}.$ (2)

This result already shows explicitly the deep link between the Voronoi
construction and gamma distributions.

Model IS.1 (fixed $N$, variable $C$ on the circle) —Now consider a circular
Volvox. Select a fixed number $N_{;\lambda}=n+1$ of successive points
$\\{R_{i}\\}_{i=k}^{k+n+1}$ from the half-line process in IS.0 and form a
circle by identifying the first and last points, as in Fig. 1(c). This yields
a variable circumference $C$ which is $n$-gamma distributed, i.e.
$C=Y_{;n,\lambda}$, and Voronoi lengths $V_{i}$ which by the same convolution
argument follow the identical $k=2$ gamma distribution as in IS.0. Fig. 1(d)
shows a numerical demonstration of this using $10^{4}$ realizations.

Model IS.2 (fixed $N$, fixed $C$)—Let us now constrain the circumference $C=c$
in IS.1 as follows. Take the points $\\{R_{i}\\}_{i=1}^{n+1}$, assume wlog
$R_{1}=0$, and form the normalized sequence
$cR_{1}/R_{n+1},cR_{2}/R_{n+1},\ldots,c$. These are the order statistics of
$n-1$ uniform random variables on the interval $(0,c)$ (see §4.1 of [12]), as
seen in Fig. 1(c). Identifying the first and last points as in IS.1, the
arclengths $\\{\ell_{i}\\}_{i=1}^{n}$ are the ratio of i.i.d. exponential
random variables to their sum,

$\displaystyle\ell_{i}$ $\displaystyle=$ $\displaystyle
c\,\frac{X_{i;\lambda}}{\sum_{j=1}^{n}X_{j;\lambda}}.$ (3)

Thus, the corresponding Voronoi lengths $\\{V_{i}\\}$ are

$\displaystyle V_{i}$ $\displaystyle=$
$\displaystyle\frac{c}{2}\,\frac{Y_{1;2,\lambda}}{Y_{1;2,\lambda}+Y_{2;n-2,\lambda}}=\frac{c}{2}Z_{;2,n-2}$
(4) $\displaystyle\overset{\text{pdf}}{\ \sim\ }$
$\displaystyle\frac{4v}{c(n-1)(n-2)}\left(1-\frac{2v}{c}\right)^{(n-3)},$ (5)

where $v\in\left[0,\frac{c}{2}\right]$, which is computed using the joint
distribution $(Y_{1}+Y_{2},Y_{1}/(Y_{1}+Y_{2})$ [9]. As can be seen, $V_{i}$
is $\alpha=2$ beta-distributed, taking values in $[0,\frac{c}{2}]$, which is
the second order statistic of $n-1$ uniform random variables. Taking the
limits $n,c\to\infty$ with density $\rho=m/c=(n-2)/c$ fixed—a kind of
thermodynamic limit—we obtain, by taking the pointwise limit of cumulative
density functions [9],

$\displaystyle V_{i}$ $\displaystyle=$
$\displaystyle\frac{m}{2\rho}Z_{2;m}\overset{n\to\infty}{\longrightarrow}\frac{1}{2}Y_{;2,\rho},$
(6)

converging in distribution to the $k=2$ gamma distribution (with cell
density/rate $\rho$) as in IS.0 & IS.1.

Model IS.3 (variable $N$, fixed $C$)—In contrast to IS.1, let us now constrain
the circumference $C=c$, leaving the count $N_{;c\lambda}$ a Poisson-
distributed random variable. Using IS.2, the distribution of $V_{i}$ can be
computed by marginalizing the joint distribution $(Z_{;2,n-1},N_{;\lambda})$
over $n$—the compound beta-Poisson distribution—which is a $k=2$ gamma
distribution (see SM [9]). A simpler approach is to place an observation
window $[x,x+c]$ at random on the half-line process and note that the inter-
point spacings $\ell_{i}$ follow a truncated exponential distribution:

$\displaystyle\ell_{i}$ $\displaystyle\overset{\text{pdf}}{\ \sim\ }$
$\displaystyle\frac{\lambda e^{-\lambda x}}{1-e^{-\lambda c}},\qquad
x\in[0,c]$ (7)

Taking $c\to\infty$, leaving the rate $\lambda$ fixed (analogous again to the
thermodynamic limit), we have $\ell_{i}\overset{\text{d}}{\to}X_{i;\lambda}$,
with $V_{i}$ converging to a $k=2$-gamma distribution as in IS.0.

Thus far, we have considered models in which cells in the Volvox ECM, viewed
as the final configuration of some stochastic growth process, occur as points
in space that are independently scattered. Gamma distributions govern their
Voronoi tessellations in various limits, independent of the constraints
applied to the scattering.

Model IS.4 (fixed $N$, fixed $C$, no overlaps)—Starting with fixed $C=c$ and
$N_{;c\lambda}=n$ per IS.2, now suppose that points have a finite “diameter”
$v_{c}$. In $d=1$, the shape of a segment $\ell_{i}$ is invariant of the
allocation $v_{c}$; hence, we may start with $v_{c}$-length segments and add
the random spacings of a subinterval of length $(c-nv_{c})$. Using IS.2,
$V_{i}$ therefore is a $v_{c}$-shifted beta random variable, converging in the
limits $n,c\to\infty$ with fixed density $\rho=(n-2)/(c-nv_{c})$ to a shifted
gamma random variable,

$V_{i}=\frac{n-2}{2\rho}Z_{;2,n-2}+v_{c}\overset{n\to\infty}{\longrightarrow}\frac{1}{2}Y_{;2,\rho}+v_{c}.$
(8)

Therefore, (8) has the distributional form

$f_{Y_{(k,\lambda)}+v_{c}}(v)=\frac{\lambda^{k}}{\Gamma(k)}(v-v_{c})^{k-1}e^{-\lambda(v-v_{c})},$
(9)

which is precisely the empirical law (1) found for Volvox with
$\lambda=k/(\bar{v}-v_{c})$. These hard-core point processes (including the
Matérn and Gibbs processes, 5.4-5.5 [14]) relax the independence property
characterizing Poisson processes by allowing energetic interactions—and
therefore correlations—between points.

Figure 2: Dynamic models of cellular positions during growth converge to gamma
distributions in various limits. (a1,a2) Standard Brownian motion in angular
coordinates converges to beta-distributed segments in large $t$, and gamma-
distributed segments in large $n$. (b1,b2) Dyson Brownian motion in angular
coordinates converges to gamma-distributed segments at large $t$ and $n$. Non-
conformance to gamma distributions is observed at low $n$ due to pair-
repulsion. (c1,c2) Dyson Brownian motion on a growing domain satisfying
particular growth constraints converges in large-$t$ and $n$ to gamma-
distributed segments.

One-dimensional dynamic models. In IS.1-IS.4, fluctuations in Voronoi segments
arise from random cell positions representing static configurations, with
different combinations of fixed and free random variables. Focusing now on a
single organism with fixed $N_{;c(0)\lambda}=n$ cells on a time-varying
circumference $c(t)$, we now consider fluctuations in Voronoi segments arising
from random motions of cells during ECM generation, represented by identically
distributed stochastic processes $(\\{\ell_{t}^{(i)}\\}_{i=1}^{n})_{t\geq 0}$.
Physically, the random motion of cell locations arises from stochastic ECM
production. Initially, let us consider a simplified view in which the ECM is
treated like a fluid, and cell motions obey diffusive dynamics without
constraints on inter-cell distances. Gamma distributions are shown to arise
once again in the large-$t$ and $n$ limits.

Model ID.1 (Brownian motion on the circle)—Let $\\{R_{t}^{(i)}\\}_{i=1}^{n}$
be $n$ Brownian motions (BMs) on a circle of radius $r$ given by

$\displaystyle R_{t}^{(j)}$ $\displaystyle=$ $\displaystyle
r\exp(i\theta_{t}^{(j)}),$ (10)

with $\theta_{t}^{(j)}$ standard BMs on $\mathbb{R}$ and initial conditions
for the probability density $p^{(j)}(0)=\delta_{R_{0}^{(j)}}$. The time-
depenent probability density $p^{(i)}(t)$ tends exponentially in $L^{2}$ to
the uniform distribution, hence $R_{t}^{(i)}$ converge in distribution to
i.i.d. uniform r.v.s. Therefore, in large $t$ the configuration approaches
Model IS.2, and $V_{i}$ are beta-distributed (4) in large $t$ and gamma-
distributed in large $t,n$ (6)—see Fig. 2(a1,a2). Yet, the sample paths of
(10) almost surely intersect, unphysically reordering the cells.

Model ID.2 (noncolliding Brownian motion)—Let $\\{R_{t}^{(i)}\\}_{i=1}^{n}$ be
samples of the conditional distributions of $n$ circular Brownian motions
whose angles are in ascending order modulo $2\pi$, thereby lying in the set

$\displaystyle D_{n}$ $\displaystyle=$ $\displaystyle\\{x\in\mathbb{R}^{n}\ |\
x_{1}<\cdots<x_{n}<x_{1}+2\pi\\},$ (11)

a construction known as noncolliding Brownian motion or Brownian motion within
the Weyl chamber $D_{n}$ [15]. In [16] eq. 4.1, it is shown that the
conditional fluctuations are Gaussian plus a singular $r^{-1}$ pair-repulsion,

$\displaystyle d\theta_{t}^{(i)}$ $\displaystyle=$ $\displaystyle\sigma
dB_{t}^{(i)}+\frac{\sigma^{2}}{2}\sum_{j\neq
i}\cot\left(\frac{\theta_{t}^{(i)}-\theta_{t}^{(j)}}{2}\right)dt.$ (12)

A physical interpretation of (12) is that of a gas confined to the unit circle
with the pair-potential [9]

$\displaystyle W$ $\displaystyle=$
$\displaystyle-\sum_{j<k}\log\left|\exp(i\theta_{k})-\exp(i\theta_{j})\right|,$
(13)

constituting a simple model of non-colliding cell motion. Eq. 12 is precisely
the eigenvalue dynamics $\lambda_{t}^{(j)}=\exp(i\theta_{t}^{(j)})$ of a
Brownian motion $U_{t}$ on the unitary group $\mathbb{U}(n)$, known as Dyson
Brownian motion [17, 9].

Being confined to $D_{n}$ (11), the positions $R_{t}^{(i)}$ do not converge to
uniform r.v.s on the circle as in I.5—compare Figs. 2(a1) and 2(b1). The
stationary distribution of (12) is the circular unitary ensemble (CUE) [17],

$\rho_{*}(\theta_{1},\cdots,\theta_{n})=\frac{1}{Z_{n}}\prod_{j<k}|\exp(i\theta_{j})-\exp(i\theta_{k})|^{2}.$
(14)

We use this result to derive the Voronoi length distribution in the
large-$t,n$ limits. Let $\\{\theta^{(i)}\\}_{i=1}^{n}$ be a sample of the
stationary distribution (14), the spectrum of a uniform sample of
$\mathbb{U}(n)$; the empirical distribution
$\mu_{\theta}(n)=n^{-1}\sum_{i=1}^{n}\delta_{\theta^{(i)}}$ converges almost
surely in large $n$ (Theorem 3, [18]) to the uniform distribution on the unit
circle. Thus, the spacings $\theta^{(i+1)}_{t}-\theta^{(i)}_{t}$ converge in
large-$t$ and $n$ (in either order [9]) to the spacings of the order
statistics of $n$ uniform random variables as in IS.2. Taking $n\to\infty$
with density $\rho=(n-2)/c$ fixed, $V_{i}$ converge to $k=2$ gamma random
variables as in (6).

Convergence to gamma laws depending on particle count $n$ is shown for BM and
DBM in Figs. 2(a2,b2) [9]; in contrast to BM, DBM spacings lose the long tail
at low $n$ due to repulsion (14).

Model ID.3 (noncolliding Brownian motion, growth)—To account for growth of the
ECM, let the radius $r(t)$ in ID.2 be time-dependent, scaling the
configuration as $R_{t}^{(i)}=r(t)\exp(i\theta_{t}^{(i)})$ with
$\theta_{t}^{(i)}$ given by (12) as shown in Fig. 2(c1). Applying Itô’s lemma,
$R_{t}^{(i)}$ satisfies the stochastic differential equation

$dR_{t}^{(i)}=\dot{r}\exp(i\theta_{t}^{(i)})dt+iR_{t}^{(i)}d\theta_{t}^{(i)}-\frac{r}{2}\theta_{t}^{(i)}dt.$
(15)

Substituting $d\theta_{t}^{(i)}$ from (12), the diffusion constant for (15) is
$D=r^{2}\sigma^{2}/2$. Requiring that the lateral diffusion is time invariant
implies $\dot{D}=0$, and thus that the standard deviation of the radius-
normalized dynamics (12) should decay as $\sigma(t)\propto r(t)^{-1}$.
Exponential convergence of (12) to the stationary solution (14) is ensured by
a growth condition on $r(t)$. By standard Fourier arguments (e.g. §2.2, [19],
and [9]), solutions to the time-dependent diffusion equation for the
probability density of $p(t)$ of the unitary Brownian motion $U_{t}$ satisfy,
for some constant $C_{n}>0$ depending only on $n$,

$d(p(t),\nu)\leq
d(p(0),\nu)\exp\left(-C_{n}\int_{0}^{t}\sigma(s)^{2}ds\right),$ (16)

with $d$ the $L^{2}$ metric and $\nu$ the uniform (Haar) measure on the
unitary group $\mathbb{U}(n)$. Therefore, if $r(t)$ satisfies:

$\displaystyle\lim_{t\to\infty}\int_{0}^{t}\frac{1}{r(s)^{2}}ds$
$\displaystyle=$ $\displaystyle\infty,\qquad r(0)>0$ (17)

the RHS in (16) vanishes in large $t$, and global exponential stability is
assured. Then, as in ID.2, the positions $R_{t}^{(i)}$ approach a random
uniform spacing of a circle of radius $r(t)$ in large $n$, resulting in gamma-
distributed Voronoi segments $V_{i}$. Rapid convergence in $t$ to a gamma law
(with $D=0.1,n=1000,\dot{r}=1$) is shown in Fig. 2(c2). Formally, (17) is not
satisfied by such a linear growth law, yet we observe empirical convergence to
a gamma law on the order of a unit of (scaled) time.

Condition (17) can be understood by a Péclet number relating drift and
diffusion timescales in growth. Considering the radial drift velocity $v_{r}$
in (15), let

$\text{Pe}=\frac{\tau_{d}}{\tau_{v_{r}}}=\frac{2L\dot{r}}{r^{2}\sigma^{2}}\propto
L\dot{r}$ (18)

with $L$ a test length section. In the limit $\text{Pe}\to 0$, condition (17)
is satisfied; when $\text{Pe}\to\infty$, the exponential multiplier in the
bound (16) approaches $1$, “freezing” the angles $\theta^{(i)}_{t}$ to initial
conditions. Finally, conditions $\dot{D}=0$, $\dot{\text{Pe}}=0$, and (17)
cannot simultaneously be satisfied.

Let us now consider the ECM as a solid, with cell motions arising from
stochastic growth of the inter-cell segments. As a simple starting point,
suppose that the time-varying intercell lengths $\ell_{t}^{(i)}$ are
stochastic processes with i.i.d. nonnegative increments $G_{t}$ for
$t=0,1,\ldots$, representing a strictly additive growth process in discrete
time. Naturally if $G_{t}$ have finite second moment, by the central limit
theorem,

$\displaystyle\frac{\ell_{t}^{(i)}}{\sqrt{t}}-\sqrt{t}\mathbb{E}[G_{t}]\ $
$\displaystyle\overset{t\to\infty}{\longrightarrow}\
W_{;0,\text{Var}(G_{t})},$ (19)

where $W$ is a Gaussian random variable. For example, exponentially
distributed fluctuations $G_{t}=X_{;\lambda}$, as in adding the spacings of
$T$ successive Poisson processes, results in $(k=T)$-gamma-distributed
segments $\ell_{i}$, converging in large $T$ to a translated Gaussian random
variable. This convergence to a Gaussian distribution restates the point made
earlier [5] regarding the surprising appearance of gamma distributions and, in
light of the results from static models above, indicates that not every random
growth process produces snapshots of the configurations that conform to gamma
distributions.

Figure 3: Poisson-Voronoi tessellation on the periodic unit square. Panels
(a1) - (d1) show the Voronoi cells colored by four geometric quantities
defined for closed plane curves (see [9]). Panels (a2) - (d2) show empirical
and maximum-likelihood estimations of gamma (denoted $Y$), lognormal (denoted
$\exp(W)$), and beta prime (denoted $Z^{\prime}$) distributions where
applicable. Poisson-Voronoi tessellations do not have gamma-distributed aspect
ratios, instead following an approximate beta-prime law.

Model ID.4 (maximum-entropy growth rates)—Let the segments grow linearly in
time as $\dot{\ell}_{i}=G_{i;\mu}$, with i.i.d. growth rates $G_{i}$ which are
positive continuous random variables with some common mean growth rate $\mu$.
If the distribution of $G_{i}$ maximizes entropy subject to the mean and
nonnegativity constraint—perhaps more interpretable as uncertainty about
cellular behavior than global maximum-entropy assumptions [8]—then
$G_{i;\mu}=X_{i;1/\mu}$ is an exponential random variable of rate
$\lambda=1/\mu$ (see [9]). Then, for any time $t$, the normalized
configuration

$\displaystyle\frac{\ell_{i}(t)}{\sum_{j=1}^{n}\ell_{j}(t)}$
$\displaystyle=\frac{\ell_{i}(0)+tG_{i}}{\sum_{j=1}^{n}\ell_{j}(0)+tG_{j}}\overset{\text{d}}{\to}\frac{X_{i;1/\mu}}{\sum_{j=1}^{n}X_{j;1/\mu}}$
(20)

converges in large $t$ to a uniform spacing as in IS.2, (3) and therefore has
beta-distributed Voronoi lengths (4) converging to gamma-distributed lengths
in large $n$ (6).

Two-dimensional equilibrium and dynamic models. In dimension $d\geq 2$, known
analytical results concerning Voronoi tessellations of point processes are
limited to lower-dimensional facets, such as edge (2D) and face (3D)
distributions (see [23] and §4.4, [24]). The difficulty in higher dimensions
arises primarily from the loss of uniqueness for shapes satisfying geometric
properties—such as fixed measure (length, area, volume)—for which one must
first consider distributions over shapes, then marginalize over the level-sets
satisfying a scalar geometric property, such as aspect ratio. For this reason,
in the following models we present primarily numerical results (with partial
analytical arguments where applicable), and present a validation of the
numerical method in the Supplementary Material [9].

Model IIS.1 (periodic unit square)—Consider a homogeneous Poisson point
process on the unit square $[0,1]^{2}$ with periodic boundary conditions,
denoted $\mathbb{T}^{2}$. On general $d$-dimensional domains, the point
process is specified by an intensity measure $\lambda(A)$ for subsets
$A\subseteq\mathbb{T}^{2}$ in which the count $N_{A;\lambda(A)}$ is a Poisson
random variable of rate $\lambda(A)$ (see [9] for sampling methods). A
homogeneous Poisson process —one whose intensity $\lambda$ is constant on sets
$A$ of constant measure—is realizable by sampling the total count
$N_{;\lambda(\mathbb{T}^{2})}$ and assigning the positions
$\\{R_{i}\\}_{i=1}^{n}$ conditional on $N=n$ as i.i.d. uniform random
variables. Figure 3 shows numerical simulations for $k=1000$ trials with
intensity $\lambda(\mathbb{T}^{2})=10^{3}$, which is on the order of the
number of somatic cells in Volvox carteri.

The periodic Voronoi tessellation on $\mathbb{T}^{2}$, shown in the small-$n$
example in Figure 3(a1)-(d1), is constructed by copying $\\{R_{i}\\}$ in four
quadrants around $\mathbb{T}^{2}$ and selecting the sub-tessellation
corresponding to the original points. Areas $a$, nondimensionalized as
$\lambda a$, conform to a gamma random variable with $k\approx 3.5$, as in
Fig. 3(a2). The isoperimetric deficit $D=L/\sqrt{4\pi a}-1$ (see [9] and [37])
with $L$ the perimeter, conforms after an appropriate rescaling to a log-
normal random variable with $\sigma\approx 0.6$, shown in Fig. 3(c2). The
aspect ratio AR of a Poisson-Voronoi tessellation does not conform to a gamma
distribution (Fig. 3(b2)), in contrast to prior studies [20] finding gamma-
distributed aspect ratio in a diverse range of densely-packed cell types and
inert matter. Instead, AR conforms approximately to a beta-prime distribution,
which naturally arises as the ratio of independent gamma random variables [9].
Figure 3(d2) shows that the major axis length is gamma-distributed; while we
do not in general expect the major and minor axis to be independent random
variables, we show in the SM [9] that both axes are gamma-distributed, with
estimated parameters within the same order of magnitude as the estimated beta-
prime distribution for aspect ratio.

Figure 4: Spherical point processes with a hard-core repulsion demonstrate a
transition to gamma laws in geometric features, such as aspect ratio and
isoperimetric deficit, not normally found in Poisson-Voronoi tessellations.
Panel (a2) shows a realization of a Matérn Type-II process with the minimum
distance $r_{d}$ set to twice the apothem of a regular spherical hexagon of
area $4\pi/\lambda$.

Model IIS.2 (sphere with hard-core repulsion)—Consider a homogeneous Poisson
process of rate $\lambda$, thinned (see [9]) to some rate
$\lambda^{\prime}<\lambda$ such that all points $\\{R_{i}\\}_{i=1}^{n}$ have a
minimum pairwise distance $d(R_{i},R_{j})\geq r_{d}$ where $d(\cdot,\cdot)$ is
the geodesic distance on the sphere. This generalizes IS.4 to the two-
dimensional case, in which the geometric realization of such a hard-core
repulsion is not unique, a fact intimately linked to the existence of an
ordering of points in 1 dimension. Realizing a minimum spacing requires
dependent thinning, as the inclusion of a particular point $R_{i}$ depends
upon the locations of nearest-neighbors $R_{j}$, and an independently thinned
Poisson process is again a Poisson process ([9]). A commonly used procedure is
the Matérn Type-II thinning rule (see [9] and 5.4, [14]):

$\displaystyle\\{R^{\prime}_{i}\\}$ $\displaystyle=\\{R_{i}\ |\ M_{i}<M_{j}\
\forall j\neq i,\ d(R_{i},R_{j})\leq r_{d}\\}$ (21)

where $M_{i}$ are i.i.d. uniform random variables called marks. Heuristically,
(21) imposes an upper bound on the aspect ratio, as the “minor axis” has a
lower-bounded arclength $r_{d}$. We expect this to result in vanishing tails
in both aspect ratio and isoperimetric deficit, confirmed by numerical
simulation with $k=100$ trials of $\lambda=2000$, with $r_{d}$ twice the
apothem of a regular spherical hexagon (see [9]) of area $4\pi/\lambda$. The
convergence to approximate gamma laws in area, aspect ratio, and anisotropy-
adjusted isoperimetric deficit is illustrated in Figure 4. Notably, the
isoperimetric deficit is no longer well-approximated by a log-normal law as in
the Poisson case (Figure 3).

Model IID.1 (Brownian motion on the sphere)—Analogous to ID.1, consider a
fluid ECM in which cells $R_{t}^{(i)}$ undergo diffusive motion without
constraints on the inter-cell distances. Identically to the 1-dimensional
case, $R_{t}^{(i)}$ convergence exponentially in distribution to iid uniform
random variables on the sphere, which is the conditional distribution of a
Poisson process. The resulting statistics on the sphere (Figure 4 (a1)-(d1))
closely follow those of the flat torus (Figure 3), with only area well-
approximated by a gamma law.

A conjectured distribution. In dimension $d\geq 2$, exact distributions for
the cell measures (areas, volumes, and so on) of Delaunay triangulations
(denoted $D_{i}$) of a Poisson point process are known [22], with the
distribution in $\mathbb{R}^{2}$ given by a modified Bessel function of the
second kind,

$\displaystyle f_{D_{i}}(v)$ $\displaystyle=$ $\displaystyle
c_{1}v^{\alpha_{1}}K_{n}(c_{2}v^{\alpha_{2}}),$ (22)

with parameters $c_{i},\alpha_{j},n$. On the other hand, exact distributions
for their vertex-cell duals—the Voronoi tessellations—are presently unknown.
We conjecture, however, that (22) should also govern the Voronoi areas $V_{i}$
based on the following heuristic argument. Conditional independence of the
volume of the fundamental region (a particular set containing the vertices of
a Voronoi cell and the origin) and its shape (see [23]) suggests that one may
assume a particular approximate shape for the typical Voronoi cell, say, a
disk or regular hexagon. If its radius $L_{i}$ is gamma-(inclusively,
exponentially-)distributed, then $V_{i}$ is proportional to $L_{i}^{2}$, whose
distribution is precisely of the form (22) (see [25]).

The universality of gamma distributions in random Voronoi
tessellations—particularly beyond the Poisson-Voronoi case, and in geometric
features beyond area—is remarkable. Concerning area, let us remark on this
observation using a well-known characterization of the gamma distribution
which seems to be analogous to the independent scattering property
characterizing Poisson processes [9]. Let $A$ be the typical partition size of
a random tessellation and $A_{c}$ be the size of the remainder. Then, the
total $T=A+A_{c}$ and fraction $F=A/(A+A_{c})$ are independent if and only if
$A$ and $A_{c}$ are independent gamma random variables of the same rate
parameter (Lukacs’ theorem, [13]). This recapitulates the story of gamma laws
describing entropy-maximizing configurations [8]; random tessellations which
maximize the joint entropy $H(T,F)$—where geometrically realizable, such as
one-dimensional settings—have gamma-distributed partition sizes.

Even more remarkably, the emergence of these distributions in biological
tissues, such as areas in Volvox and the volumes in “snowflake yeast” [5], as
well as aspect ratios in a wide range of confluent tissues and inert matter
[20] suggests that point processes are suitable idealized models for a
quantitative description of the structure of a multitude of biological
systems. In this sense, they play a role analogous to the idealized random
walk models of polymers and the hard sphere model of fluids.

###### Acknowledgements.

We are grateful to T. Day, P. Yunker, and W. Ratcliff for numerous
discussions. This work was supported in part by the Cambridge Trust (AS), The
John Templeton Foundation and Wellcome Trust Investigator Grant 207510/Z/17/Z
(REG).

## References

  * [1] M.J. Dayel, R.A. Alegado, S.R. Fairclough, T.C. Levin, S.A. Nichols, K. McDonald, and N.King, Cell differentiation and morphogenesis in the colony-forming choanoflagellate Salpingoeca rosetta, Dev. Biol. 357, 73–82 (2011).
  * [2] S.R. Fairclough, M.J. Dayel and N. King, Multicellular development in a choanoflagellate, Curr. Biol. 20, R875–R876 (2010).
  * [3] A. Hallmann, Morphogenesis in the Family Volvocaceae: Different Tactics for Turning an Embryo Right-side Out, Protist 157 445–461 (2006).
  * [4] E.E. Ruppert, R.S. Fox, and R.D. Barnes, Invertebrate Zoology, 7th ed. (Thomson-Brooks/Cole, Belmont, CA, 2004).
  * [5] T. C. Day, S. S. Höhn, S. A. Zamani-Dahaj, D. Yanni, A. Burnetti, J. Pentz, A. R. Honerkamp-Smith, H. Wioland, H. R. Sleath, W. C. Ratcliff, _et al._ , Cellular organization in lab-evolved and extant multicellular species obeys a maximum entropy law, eLife 11, e72707 (2022).
  * [6] W.C. Ratcliff, R.F. Denison, M. Borrello and M. Travisano, Experimental evolution of multicellularity, Proc. Natl. Acad. Sci. USA 109, 1595–1600 (2012).
  * [7] S. Höhn, A. Honerkamp-Smith, P.A. Haas, P. Khuc Trong, and R.E. Goldstein, Dynamics of a Volvox Embryo Turning Itself Inside Out, Phys. Rev. Lett. 114, 178101 (2015).
  * [8] T. Aste and T. Di Matteo, Emergence of Gamma distributions in granular materials and packing models, Phys. Rev. E 77, 021309 (2008).
  * [9] See Supplemental Material at http://link.aps.org/supplemental/xxx for analytical and numerical details.
  * [10] R. Durrett, _P_ robability: Theory and Examples, Vol. 49 (Cambridge University Press, 2019).
  * [11] D. J. Daley, D. Vere-Jones, _et al._ , An Introduction to the Theory of Point Processes: Volume I: Elementary Theory and Methods (Springer, 2003).
  * [12] R. Pyke, Spacings, J. R. Stat. Soc. B 27, 395–436 (1965).
  * [13] E. Lukacs, A Characterization of the Gamma Distribution, Ann. Math. Stat. 26, 319–324 (1955).
  * [14] D. Stoyan, W. S. Kendall, S. N. Chiu, and J. Mecke, Stochastic Geometry and its Applications (John Wiley & Sons, 2013).
  * [15] D. J. Grabiner, Brownian motion in a Weyl chamber, non-colliding particles, and random matrices, Ann. IHP Prob. Stat. 35, 177–204 (1999).
  * [16] D. G. Hobson and W. Werner, Non-Colliding Brownian Motions on the Circle, Bull. London Math. Soc. 28, 643–650 (1996).
  * [17] F. J. Dyson, A Brownian-Motion Model for the Eigenvalues of a Random Matrix, J. Math. Phys. 3, 1191–1198 (1962).
  * [18] P. Diaconis and M. Shahshahani, On the eigenvalues of random matrices, J. Appl. Prob. 31, 49–62 (1994).
  * [19] T. Kemp, Heat Kernel Empirical Laws on $\mathbb{U}(n)$ and $\mathbb{GL}(n)$, J. Theor. Prob. 30, 397–451 (2017).
  * [20] L. Atia, et al., Geometric constraints during epithelial jamming, Nat. Phys. 14, 613–620 (2018).
  * [21] R. Li, et al., E2 and gamma distributions in polygonal networks, Phys. Rev. Res. 3, L042001 (2021).
  * [22] P. Rathie, On the volume distribution of the typical Poisson–Delaunay cell, J. App. Prob. 29, 740–744 (1992).
  * [23] J. Møller and S. Zuyev, Gamma-type results and other related properties of Poisson processes, Adv. Appl. Prob. 28, 662–673 (1996).
  * [24] J. Møller, Lectures on random Voronoi tessellations, Vol. 87 (Springer Science & Business Media, 2012).
  * [25] H. J. Malik, Exact Distribution of the Product of Independent Generalized Gamma Variables with the Same Shape Parameter, Ann. Math. Stat. 39, 1751–1752 (1968).
  * [26] E. Gilbert, Random Subdivisions of Space into Crystals, Ann. Math. Stat. 33, 958–972 (1962).
  * [27] M. Tanemura, Statistical distributions of poisson voronoi cells in two and three dimensions, FORMA-TOKYO- 18, 221 (2003).
  * [28] J.-S. Ferenc and Z. Néda, On the size distribution of Poisson Voronoi cells, Physica A 385, 518–526 (2007).
  * [29] D. Weaire, J. Kermode, and J. Wejchert, On the distribution of cell areas in a Voronoi network, Phil. Mag. B 53, L101–L105 (1986).
  * [30] J. van den Berg, R. Meester, and D. G. White, Dynamic Boolean models, Stoch. Proc. Appl. 69, 247–257 (1997).
  * [31] R. Miles, Random points, sets and tessellations on the surface of a sphere, Ind. J. Stat. A 33, 145–174 (1971).
  * [32] M. Van Den Berg and J. Lewis, Brownian motion on a hypersurface, Bull. Lon. Math. Soc. 17, 144–150 (1985).
  * [33] E. Cépa and D. Lépingle, Brownian particles with electrostatic repulsion on the circle: Dyson’s model for unitary random matrices revisited, ESAIM: Prob. Stat. 5, 203–224 (2001).
  * [34] M. L. Mehta, Random Matrices (Elsevier, 2004).
  * [35] G. Menon, Lesser known miracles of Burgers equation, Acta Math. Sci. 32, 281–294 (2012).
  * [36] P. Findeisen, A simple proof of a classical theorem which characterizes the gamma distribution, The Annals of Statistics, 1165 (1978).
  * [37] N. Fusco, F. Maggi, and A. Pratelli, The sharp quantitative isoperimetric inequality, Ann. Math. 941, (2008).
  * [38] S. Diamond and S. Boyd, CVXPY: A Python-embedded modeling language for convex optimization, The Journal of Machine Learning Research 17, 2909 (2016).
  * [39] D. Kraft, A software package for sequential quadratic programming, Forschungsbericht- Deutsche Forschungs- und Versuchsanstalt fur Luft- und Raumfahrt (1988).

## Supplemental Material

This file contains analytical and computational details pertaining to the
models discussed in the main body.

##  I Random variables

###  I.1 Transforms and convergence

######  Definition S1 (Pushforward measure).

Let $g:\mathbb{R}^{n}\to\mathbb{R}^{n}$ be a diffeomorphism and $Y=g(X)$. The
pushforward probability measure $\mu_{Y}$ is, for all measurable
$U\subset\mathbb{R}^{n}$,

$\displaystyle\mu_{Y}(U)$ $\displaystyle=\mu_{X}(g^{-1}(U)).$ (S1)

##### Transforms of random variables.

When $\mu_{X}$ has a Radon-Nikodym derivative $f_{X}$ with respect to the
Lebesgue measure —i.e. is expressible by the probability density
$\mu_{X}(U)=\int_{U}f_{X}(x)dx$—then by the rule for integration under a
diffeomorphic change of coordinates $y=g(x)$,

$\displaystyle\mu_{X}(U)(g^{-1}(U))$
$\displaystyle=\int_{g^{-1}(U)}f_{X}(x)dx=\int_{U}f_{X}(g^{-1}(y))\left|\det\frac{\partial
g^{-1}}{\partial y}(y)\right|dy.$ (S2)

As this is true for all $U$ we conclude

$\displaystyle f_{Y}(y)$
$\displaystyle=f_{X}(g^{-1}(y))\left|\det\frac{\partial g^{-1}}{\partial
y}(y)\right|.$ (S3)

When preferable to work with $g$ rather than $g^{-1}$, we may apply the chain
rule to $g^{-1}\circ g=1$ to convert (S3) to

$\displaystyle f_{Y}(g(x))$
$\displaystyle=f_{X}(x)\left|\det\left(\frac{\partial g}{\partial
x}(x)\right)^{-1}\right|.$ (S4)

Due to this fact, we will abbreviate the scaling factor as $J^{-1}$, denoting
the inverse Jacobian determinant. For affine transforms $Y=cX+b$, (S3) becomes

$f_{Y}(y)=\frac{1}{c}f_{X}\left(\frac{y-b}{c}\right).$ (S5)

##### Sums of random variables.

Let $X,Y$ be independent random variables taking values in
$U_{X},U_{Y}\subseteq\mathbb{R}$. Then their sum is distributed as the
convolution

$\displaystyle X+Y=Z$ $\displaystyle\sim f_{Z}(z)=\int_{x\in U_{x},\ x\leq
z}f_{X}(x)f_{Y}(z-x)dx=f_{X}*f_{Y}.$ (S6)

This can alternatively be deduced by applying the transform rule (S3) to the
map $(X,Y)\mapsto(X,X+Y)$.

######  Definition S2 (Convergence in distribution).

A sequence of random variables $\\{X_{n}\\}$ taking values in an interval
$U\subseteq\mathbb{R}$ is said to converge in distribution if, for all $x$ at
which the cumulative distribution function $F_{X}$ is continuous, the c.d.f.s.
$F_{X_{n}}$ converge pointwise, i.e.

$\displaystyle\lim_{n\to\infty}F_{X_{n}}(x)$ $\displaystyle=F_{X}(x).$ (S7)

###  I.2 Basic properties of the exponential and gamma random variables

Let us recall with proof several basic properties of the gamma (and
exponential, which is a special case) random variable.

######  Definition S3 (Gamma random variable.).

The gamma random variable $Y_{k,\lambda}$ with shape parameter $k>0$ and rate
$\lambda>0$ is the continuous random variable with probability density
function

$\displaystyle f_{Y}(y)$ $\displaystyle=\frac{\lambda^{k}y^{k-1}\exp(-\lambda
y)}{\Gamma(k)}.$ (S8)

Its namesake is the gamma function

$\displaystyle\Gamma(k)$ $\displaystyle=\int_{0}^{\infty}y^{k-1}e^{-y}dy.$
(S9)

######  Lemma S1 (Gamma random variables are closed under addition).

The sum of two independent gamma random variables $Y_{1},Y_{2}$ with
$k_{1},k_{2}$ of common rate $\lambda$ is $(k_{1}+k_{2})$-gamma distributed.

###### Proof.

Since $Y_{1},Y_{2}\in[0,\infty)$, by (S6) we have

$\displaystyle Y_{1}+Y_{2}$ $\displaystyle\overset{\text{pdf}}{\ \sim\
}f_{Y_{1}}*f_{Y_{2}}=\frac{\lambda^{k_{1}+k_{2}}}{\Gamma(k_{1})\Gamma(k_{2})}\int_{0}^{z}y^{k_{1}-1}e^{-\lambda
x}(z-y)^{k_{2}-1}e^{-\lambda(z-y)}dy.$ (S10)

Taking the Beta function

$\displaystyle\mathrm{B}(k_{1},k_{2})$
$\displaystyle=\frac{\Gamma(k_{1})\Gamma(k_{2})}{\Gamma(k_{1}+k_{2})}=\int_{0}^{1}t^{k_{1}-1}(1-t)^{k_{2}-1}dt$
(S11)

with the change of variable $t=y/z$ yields

$\displaystyle=\frac{1}{z^{k_{1}+k_{2}-1}}\int_{0}^{z}y^{k_{1}-1}(z-y)^{k_{2}-1}dx.$
(S12)

Substituting into (S10), we obtain

$\displaystyle\eqref{eq:sum_gamma}$
$\displaystyle=\frac{\lambda^{k_{1}+k_{2}}z^{k_{1}+k_{2}-1}e^{-\lambda
z}}{\Gamma(k_{1}+k_{2})}\ \overset{\text{pdf}}{\ \sim\ }\
Y_{;k_{1}+k_{2},\lambda}.$ (S13)

∎

######  Corollary.

The sum of $n$ iid exponential random variables of rate $\lambda$ are gamma-
distributed with shape parameter $n$ and rate $\lambda$.

######  Lemma S2 (Beta-gamma convergence).

Let $Z_{2,m}$ be a beta random variable. Then $\frac{m}{\alpha}Z_{2,m}$
converges in distribution to the gamma random variable $Y_{2,\alpha}$ of rate
$\alpha$.

###### Proof.

By direct calculation,

$\displaystyle\frac{m}{\alpha}Z_{2,m}\overset{\text{cdf}}{\ \sim\
}\mathrm{P}\left[Z_{2,m}\leq\frac{z\alpha}{m}\right]$
$\displaystyle=\frac{1}{\mathrm{B}(2,m)}\int_{0}^{\frac{z\alpha}{m}}t(1-t)^{m-1}dt\qquad
z\in\left[0,\frac{m}{\alpha}\right],$ (S14)

and with the change of variables $s=1-t$,

$\displaystyle=\frac{1}{\mathrm{B}(2,m)}\int_{1-\frac{\alpha
z}{m}}^{1}(1-s)s^{m-1}ds$ (S15)
$\displaystyle=m(m+1)\left[\frac{1}{m}s^{m}-\frac{1}{m+1}s^{m+1}\right]_{1-\frac{z\alpha}{m}}^{1}$
(S16)
$\displaystyle=1-\left((m+1)\left(1-\frac{z\alpha}{m}\right)^{m}-m\left(1-\frac{z\alpha}{m}\right)^{m+1}\right)$
(S17) $\displaystyle=1-\left(1-\frac{z\alpha}{m}\right)^{m}(1+z\alpha)$ (S18)
$\displaystyle\overset{m\to\infty}{\longrightarrow}1-e^{-z\alpha}(1+z\alpha)\qquad
z\in[0,\infty).$ (S19)

Recalling the cumulative density function for $Y_{2,\alpha}$ (S8),

$\displaystyle Y_{2,\alpha}$ $\displaystyle\overset{\text{cdf}}{\ \sim\
}\int_{0}^{z}\alpha^{2}ye^{-\alpha y}dy=-\alpha ye^{-\alpha
y}\big{\rvert}_{0}^{z}+\int_{0}^{z}\alpha e^{-\alpha
y}dy=1-e^{-z\alpha}(1+z\alpha).$ (S20)

Hence, we have the distributional convergence
$\frac{m}{\alpha}Z_{2,m}\overset{\text{d}}{\to}Y_{2,\alpha}$. ∎

######  Lemma S3 (Memoryless characterization of the exponential).

The only continuous random variable $X$ which (i) possesses a cumulative
distribution function $F_{X}$ such that $F^{\prime}_{X}(0)$ exists and (ii)
satisfies the memoryless property

$\displaystyle\mathrm{P}[X>x+y\ |\ X>x]$ $\displaystyle=\mathrm{P}[X>y]$ (S21)

is the exponential.

###### Proof.

By Bayes’ theorem, (S21) is equivalent to

$\displaystyle\mathrm{P}[X>x+y]$
$\displaystyle=\mathrm{P}[X>x]\mathrm{P}[X>y]\ \Leftrightarrow\
1-F(x+y)=(1-F(x))(1-F(y)).$ (S22)

Let $G=1-F$; then (S22) is $G(x+y)=G(x)G(y)$. Then for all $z$,

$\displaystyle G^{\prime}(z)$ $\displaystyle=\lim_{h\to
0}\frac{G(z)G(h)-G(z)G(0)}{h}=G(z)G^{\prime}(0).$ (S23)

Since by hypothesis (i) $F^{\prime}(0)$ exists, this implies $F^{\prime}$
exists everywhere. Let $u$ be such that $F(u)<1$. Since $F$ is nondecreasing,
this implies $G>0$ in the interval $(-\infty,u]$. Then for $w\in(-\infty,u]$,
letting $G^{\prime}(0)=-F^{\prime}(0)=c$,

$\displaystyle c$ $\displaystyle=\frac{G^{\prime}(w)}{G(w)}=\frac{d}{dw}\log
G(w).$ (S24)

Integrating, we obtain $G(x)=b\exp(cx)$ for some $b$. By the conditions
$F(0)=0,\lim_{t\to\infty}F(x)=1$, we have $b=1,c<0$. Letting $c=-\lambda$ for
$\lambda>0$ we obtain $F_{X}(x)=1-\exp(-\lambda x)$. Hence,
$f_{X}(x)=F_{X}^{\prime}(x)=\lambda\exp(-\lambda x)$, and $X$ is an
exponential random variable of rate $\lambda$. ∎

######  Lemma S4 (Maximum-entropy characterization of the exponential).

The only nonnegative continuous random variable $X$ with density $f_{X}$ which
maximizes the entropy with fixed mean $\mu>0$ is the exponential.

###### Proof.

By hypothesis, $f_{X}$ is a critical point of the functional

$\displaystyle J[f]$
$\displaystyle=\int_{0}^{\infty}L(f(x),\lambda_{0},\lambda_{1})dx,$ (S25)

for Lagrange multipliers $\lambda_{0},\lambda_{1}$ constraining the 0th, 1st
moments, and Lagrangian density

$\displaystyle L(f(x),\lambda_{0},\lambda_{1})$ $\displaystyle=f(x)\log
f(x)+\lambda_{1}f(x)+\lambda_{1}xf(x).$ (S26)

Since for all test functions $\varphi$, the Fréchet derivative vanishes,

$\displaystyle 0$ $\displaystyle=\langle
DJ[f],\varphi\rangle=\int_{0}^{\infty}\frac{\partial L}{\partial
f(x)}\varphi(x)dx,$ (S27)

by the fundamental lemma of the calculus of variations,

$\displaystyle 0$ $\displaystyle=\frac{\partial L}{\partial f(x)}=1+\log
f(x)+\lambda_{0}+x\lambda_{1},$ (S28)

we have $f(x)=\exp(-1-\lambda_{0}-x\lambda_{1})$. Applying the total
probability constraint,

$\displaystyle 1$
$\displaystyle=\int_{0}^{\infty}\exp(-1-\lambda_{0}-x\lambda_{1})dx=\frac{1}{\lambda_{1}}\exp(-1-\lambda_{0}),$
(S29)

where $\lambda_{1}>0$ necessarily. Then,
$f_{X}(x)=\lambda_{1}\exp(-\lambda_{1}x)$ is the density of an exponential
random variable with rate $\lambda_{1}$. ∎

###  I.3 Characterization theorem for gamma random variables

######  Theorem S5 (Lukacs 1955, [13]).

Let $Y_{1},Y_{2}$ be independent random variables. Then

$\displaystyle A$ $\displaystyle=Y_{1}+Y_{2},\qquad
B=\frac{Y_{1}}{Y_{1}+Y_{2}}$ (S30)

are independent if and only if $Y_{1},Y_{2}$ are gamma random variables of the
same rate $\lambda$.

###### Proof ($\implies$).

Consider the map

$\displaystyle g(Y_{1},Y_{2})$
$\displaystyle=\left(Y_{1}+Y_{2},\frac{Y_{1}}{Y_{1}+Y_{2}}\right)=:(A,B).$
(S31)

Then for $a\neq 0$,

$\displaystyle g^{-1}(a,b)$ $\displaystyle=(ab,a-ab).$ (S32)

Therefore

$\displaystyle\left|\det\frac{\partial g^{-1}}{\partial(a,b)}\right|$
$\displaystyle=\left|\det\begin{bmatrix}b&1-b\\\ a&-a\end{bmatrix}\right|=a.$
(S33)

By (S4), the pushforward density is

$\displaystyle f_{g(Y_{1},Y_{2})}(a,b)$
$\displaystyle=af_{Y_{1}}(ab)f_{Y_{2}}(a-ab)$ (S34)
$\displaystyle=\frac{a(ab)^{k_{1}-1}(a-ab)^{k_{2}-1}}{\Gamma(k_{1})\Gamma(k_{2})}\lambda^{k_{1}+k_{2}}e^{-\lambda
a}$ (S35)
$\displaystyle=\frac{b^{k_{1}-1}(1-b)^{k_{2}-1}}{\Gamma(k_{1})\Gamma(k_{2})}a^{k_{1}+k_{2}-1}\lambda^{k_{1}+k_{2}}e^{-\lambda
a}.$ (S36)

By inspection, the total $A$ and fraction $B$ are independent. Substituting
the Beta function (S11),

$\displaystyle\eqref{eq:lukacs_joint_density}$
$\displaystyle=\frac{b^{k_{1}-1}(1-b)^{k_{2}-1}}{\mathrm{B}(k_{1},k_{2})}\frac{a^{k_{1}+k_{2}-1}\lambda^{k_{1}+k_{2}}e^{-\lambda
a}}{\Gamma(k_{1}+k_{2})}=:f_{B}(b)f_{A}(a).$ (S37)

Then $f_{A}(a)$ is a gamma distribution (as expected, S1) and $f_{B}(b)$ is a
beta distribution. Therefore,

$\displaystyle
B=\frac{Y_{1;k_{1},\lambda}}{Y_{1;k_{1},\lambda}+Y_{2;k_{2},\lambda}}$
$\displaystyle=Z_{k_{1},k_{2}}$ (S38)

is beta-distributed. ∎

Note that (S38) is independent of the rate $\lambda$ of $Y_{1},Y_{2}$. The
converse $(\impliedby)$, that gamma distributions are unique in possessing
this independence property, is not proven here, but we refer the reader to a
straightforward proof [36] relying on the fact that the gamma distribution is
uniquely determined by its moments (Exercise 3.3.25, [10]).

##  II Point processes

###  II.1 Definitions

The usual Poisson process of rate $\lambda$ on $[0,\infty)$ is realizable as
the cumulative sum of i.i.d. exponential random variables $X_{i;\lambda}$. The
construction of general point processes on a domain $K\subseteq\mathbb{R}^{d}$
is more technical, formalizing the notion of a “random almost-surely finite
subset,” for which we recall several standard definitions [11]. Throughout, we
will assume $K$ is a closed set.

######  Definition S1 (Finite point process —2.2-2.4, [11]).

A finite point process $N$ on a complete separable metric space $K$ is a
family of random variables $N(E)$ for each Borel set $E\in\mathcal{F}_{K}$,
such that, for every bounded $E$,

$\displaystyle\mathrm{P}[N(E)<\infty]$ $\displaystyle=1.$ (S39)

Formally, $N$ is a random measure (see e.g. 5.1 [11] or Chapter 9 [11]).
Without delving too deeply into this formalism, let us introduce the following
definition based upon samples of $N$.

######  Definition S2 (Simple point process).

A simple point process $N$ is one whose points are non-overlapping. In other
words, every sample of $N$ can be written as the counting measure

$\displaystyle\nu$ $\displaystyle=\sum_{i\in I}\delta_{x_{i}}$ (S40)

where $I$ is an index set and $\delta$ denotes the Dirac measure and
$\mathrm{P}[x_{i}=x_{j}]=0$ for all $i\neq j$.

######  Definition S3 (Non-atomic point process).

A nonatomic point process $N$ is one for whom the probability of realizing any
particular point $x\in K$ is zero. That is,

$\displaystyle\mathrm{P}[N(\\{x\\})>0]$ $\displaystyle=0.$ (S41)

These definitions are sufficient to define and analyze the Poisson point
process.

######  Definition S4 (Poisson point process).

A Poisson point process $N$ on a c.s.m.s. $K$ is defined by an intensity
measure $\Lambda(E)$ such that, for all Borel sets $E\in\mathcal{F}_{K}$,
$N(E)$ is a Poisson random variable, i.e.

$\displaystyle N(E)$ $\displaystyle\overset{\text{pdf}}{\ \sim\
}\frac{\Lambda(E)^{k}\exp(-\Lambda(E))}{k!}.$ (S42)

######  Example II.1 (Stationary Poisson process).

A stationary Poisson process of rate $\lambda$ is given by the intensity
measure $\Lambda(E)=\lambda\mu(E)$, where $\mu$ is the Lebesgue measure.

######  Definition S5 (Independent scattering / complete independence).

A point process $N$ satisfies the independent scattering or complete
independence property if, for all $n>1$ and disjoint Borel sets
$E_{1},\ldots,E_{n}\in\mathcal{F}_{K}$, the variables
$N(E_{1}),\ldots,N(E_{n})$ are mutually independent.

###  II.2 Characterization theorem for Poisson point processes

Immediately, we see that Poisson point processes are simple and finite if the
intensity measure $\Lambda(E)$ is given by

$\displaystyle\Lambda(E)$ $\displaystyle=\int_{E}\lambda(x)dx$ (S43)

for some function $\lambda:K\to\mathbb{R}_{+}$. Poisson point processes S4
satisfy the complete independence property S5, but remarkably, these
properties are not logically independent.

######  Theorem S1 (Prekopa 1957, Theorem 2.4.V [11]).

A point process $N$ is a non-atomic Poisson point process if and only if it is
finite S1, simple S2, non-atomic S3, and completely independent S5.

The characterization theorem S1 motivates and justifies the use of Poisson
processes (i.e. Poisson-distributed count variables) in any scenario where
points in any realization are non-interacting.

###  II.3 Sampling Poisson and Matérn point processes

######  Lemma S2.

The conditional distribution of a Poisson point process of intensity
$\lambda(x)$ on a domain $K$ is $\frac{1}{\Lambda(K)}\lambda(x)$.

Lemma S2 allows one to sample a homogeneous Poisson point process of rate
$\lambda$ by first sampling a Poisson random variable $N_{\lambda\mu(K)}=n$,
and scattering the $n$ as points as i.i.d. uniform random variables in $K$.
Point processes with a hard-core repulsion (such as the Matérn Type-II point
process [14]) can often be realized as a thinning of a Poisson point process.

######  Definition S6 (Matérn type-II hard-core point process, [14]).

Let $N$ be a homogeneous Poisson point process of rate $\lambda$. To each
point $X_{i}$ of $N$, assign a mark $M_{i}$ which are i.i.d. uniform random
variables on $(0,1)$. Then, construct the Matérn process $N^{\prime}$ with
hard-core repulsion distance $r$ as

$\displaystyle N^{\prime}$ $\displaystyle=\\{X_{i}\in N\ |\ M_{i}<M_{j}\
\forall M_{j}\in B(X_{i},r),\ i\neq j\\}.$ (S44)

The points of $N^{\prime}$ are situated at a minimum distance $r$ from one
another.

####  II.3.1 Poisson and Matérn processes on $\mathbb{S}^{d}$

Recall that a zero-mean multivariate Gaussian $W_{0,\Sigma}$ with $\Sigma$ a
$d\times d$ symmetric positive-definite covariance matrix has the ellipsoids
$\\{x\ |\ x^{\intercal}\Sigma^{-1}x=r^{2}\\}$ as level sets of constant
density. When $\Sigma=I$, this is the property of spherical symmetry. Combined
with the fact that the map $x\mapsto x\left\lVert x\right\rVert^{-1}$ is
surface-area-preserving up to a constant multiple, this yields a
computationally efficient and numerically stable method for generating iid
uniform random variables $U_{i}$ on $\mathbb{S}^{d}$:
$U_{i}=W_{i;0,I}\left\lVert W_{i;0,I}\right\rVert^{-1}$.

Lemma S2 then allows the realization of homogeneous Poisson processes on
$d$-spheres of radius $r$ as $N_{\lambda Ar^{d}}=n$, $\\{U_{i}\\}_{i=1}^{n}$,
with $A$ the surface area of the unit sphere $\mathbb{S}^{d}$ and $\lambda$
the intensity per unit area. Conveniently, the parametrization of $U_{i}$ in
Cartesian coordinates allows the realization of Matérn processes (S6) using
the geodesic (great-circle) distance $d(U_{i},U_{j})=r\arccos(U_{i}\cdot
U_{j})$.

##  III Geometry of closed curves

###  III.1 Moments of area

The area moments of simple closed curves on planar surfaces are useful for
quantifying for geometric properties of internally connected objects such as
cells. All tensorial quantities and identities below are for dimension 2.

######  Definition S1 (Planar $n$th area moment tensor).

Let $x\in D$ be a two-dimensional domain bounded by a simple closed curve
$\partial D=C$. The $n$th moment tensor of the area $D$ with respect to a
point $v$ is defined as

$\displaystyle I^{(n)}_{i_{1}i_{2}\cdots i_{n}}(v)$
$\displaystyle=\int_{D}(x_{i_{1}}-v_{i_{1}})(x_{i_{2}}-v_{i_{2}})\cdots(x_{i_{n}}-v_{i_{n}})$
(S45)

for $i_{1},\cdots,i_{n}\in\\{1,2\\}$. When the argument $v$ is dropped, it is
assumed $v=0$.

We will assume $D$ is a uniform mass-density domain in the following.

####  III.1.1 First moment, center of mass, and centrality

Applying the divergence theorem, $I^{(1)}$ becomes the surface integral

$\displaystyle I^{(1)}$
$\displaystyle=\int_{D}x=\frac{1}{3}\int_{D}\partial_{i}x_{i}x_{j}=\frac{1}{3}\int_{C}n_{i}x_{i}x_{j}$
(S46)

where $n_{i}$ is the unit boundary normal.

######  Definition S2 (Center of mass of a closed curve).

The center of mass of $D$ is

$\displaystyle\mu$ $\displaystyle=\frac{I^{(1)}}{I^{(0)}}.$ (S47)

The notation $\mu$ is used for the probabilistic interpretation as the
expected value of a uniform distribution supported on $D$.

######  Definition S3 (Variance-adjusted centrality of a test point within a
polygon).

For a test point $y$, we define the centrality metric

$\displaystyle\text{centrality}(y)$ $\displaystyle=\left\lVert
W(y-\mu)\right\rVert^{2}=(y-\mu)^{\intercal}\Sigma^{-1}(y-\mu)$ (S48)

with $W$ a matrix defined as $W^{\intercal}W=\Sigma^{-1}$, and $\Sigma$ the
second central moment defined in the next part (S50).

This whitening procedure (using terminology from the probabilistic
interpretation of $\Sigma$) enables comparison across curves $C$ of varying
moments—that is, centrality$(y)$ is dimensionless. In probability, (S48) is
the Mahalanobis distance of $y$ and $D$.

####  III.1.2 Second moment, area moment of inertia, and isotropy

Again by divergence theorem,

$\displaystyle I^{(2)}$
$\displaystyle=\int_{D}x_{i}x_{j}=\frac{1}{4}\int_{D}\partial_{i}x_{i}x_{j}x_{k}=\frac{1}{4}\int_{C}n_{i}x_{i}x_{j}x_{k}.$
(S49)

######  Definition S4 (Second moment of area).

The second central moment of area (i.e. the area moment of inertia about the
center of mass) is

$\displaystyle\Sigma$ $\displaystyle=I^{(2)}(\mu)=I^{(2)}-(\mu\otimes
I^{(1)}+I^{(1)}\otimes\mu)+I^{(0)}\mu\otimes\mu$ (S50)

where $x\otimes y=xy^{\intercal}$.

In connection to probability, $\Sigma$ is precisely the covariance matrix of a
mean-zero uniform distribution supported on the domain $D$.

######  Definition S5 (Principal axes of a polygon).

Using (S50), we may define the principal axes and stretches of a polygon
(represented by a closed curve $C$) via Hermitian eigendecomposition

$\displaystyle\Sigma$ $\displaystyle=P\Lambda P^{-1}$

with $v_{i}$ the columns of $P$ being the principal axes and and
$\lambda_{i}>0$ in ascending order (as $\Sigma$ is symmetric positive-definite
for non-degenerate curves).

We then define the “anisotropy” of the shape $D$ as the distance of the second
central moment $\Sigma$ to its $\ell^{2}$-orthogonal projection onto the
subspace $V\in\mathbb{R}^{d\times d}$ of isotropic matrices. There are two
linearly independent isotropic matrices in dimension $d=2$, the matrices
$\Sigma=\delta_{ij}$ and $\Sigma=\epsilon_{ij}$. The latter identity cannot
hold since $\Sigma_{ii}\neq 0$ for all $i$ as defined by (S49), hence we take
$V$ to be the one-dimensional subspace $\alpha I$ with $\alpha\in\mathbb{R}$.
Computing

$\displaystyle 0$
$\displaystyle=\frac{\partial}{\partial\alpha}\left\lVert\Sigma-\alpha
I\right\rVert^{2}=-2\langle\Sigma-\alpha
I,I\rangle=-2\operatorname{tr}(\Sigma)+2d\alpha,$ (S51)

then the orthogonal projection of $\Sigma$ onto $V$ is simply
$\frac{\operatorname{tr}(\Sigma)}{d}I$ as expected. Thus in dimension $d$
($d=2$ for planar curves), we define the anisotropy of $D$ (in dimensionless
units) as

$\displaystyle\text{anisotropy}(D)$
$\displaystyle=\left\lVert\frac{d}{(\operatorname{tr}\Sigma)}\Sigma-I\right\rVert^{2}.$
(S52)

Clearly anisotropy$(D)$ is nonnegative. Furthermore, anisotropy$(D)$ vanishes
if and only if $D$ is a regular polygon. Furthermore, the “orientation” of an
anisotropic domain $D$ is given by the principal axes of $\Sigma$ and
“elongation” by the principal stretches. For an ellipse, $\Sigma^{-1}=M$ is
precisely the ellipse matrix (S71).

####  III.1.3 Approximation by an ellipse.

######  Definition S6 (Aspect ratio of a polygon).

Let $\lambda_{1},\lambda_{2}$ and $v_{1},v_{2}$ be the principal stretches and
axes (in ascending order) of the second central moment $\Sigma$ (S50). Define
the aspect ratio:

$\displaystyle R$ $\displaystyle=\sqrt{\frac{\lambda_{2}}{\lambda_{1}}}$ (S53)

Then an ellipse with the same orientation, aspect ratio, and area $(\pi ab)$
as $D$ (as defined previously by its $n$th moments) has the minor and major
axes:

$\displaystyle a$ $\displaystyle=\sqrt{\frac{I^{(0)}}{\pi R}},\quad b=Ra$
(S54)

Equation (S54) may be considered a superior approximation of the polygon $D$
over e.g. a weighted $\ell^{2}$-minimization of vertex distance or a bounding
ellipse approach.

####  III.1.4 Whitening a polygon.

Let a polygon $C$ have second central moment $\Sigma$ (S50) and centroid $\mu$
(S47). With $C$ given in counterclockwise oriented coordinates as
$\\{x_{i}\\}_{i=1}^{n}$, we define the whitened polygon $C_{w}$,

$\displaystyle C_{w}$ $\displaystyle=\\{W(x_{i}-\mu)\ |\ \forall i\\},\quad
W=\Sigma^{-1/2},$ (S55)

with the matrix square root $W$ typically approximated by singular value
decomposition (SVD) as

$\displaystyle W=U(S+\varepsilon)^{-1/2}U^{*},\quad W=USV^{\intercal},$ (S56)

and $W$ then referred to as the ZCA whitening matrix with a regularization
constant $\varepsilon\ll 1$.

###  III.2 Isoperimetric inequality

We recall here several quantities which can be used as measures of the
deviation of a closed plane curve $C$ from a circle. The classical
isoperimetric inequality is

$\displaystyle 4\pi A$ $\displaystyle\leq L^{2},$ (S57)

where $A$ is the area bounded by $C$ (denoted by the set $E$) and $L$ its
total arclength, which is an equality only for circles. Accordingly, one may
define the isoperimetric quotient

$\displaystyle Q$ $\displaystyle=\frac{4\pi A}{L^{2}}\in[0,1],$ (S58)

which is maximized for circles. A natural question, however, is whether (S57)
can be used to define a set distance, in the sense of Hausdorff metric, of the
area bounded by $C$ to the “best” circle. This turns out [37] to be related to
the problem of making (S57) quantitative, in the sense of a nonnegative
quantity $\nu(E)$ such that

$\displaystyle 4\pi A+\nu(E)$ $\displaystyle\leq L^{2}$ (S59)

for all simple closed curves $C$, with $\nu(E)=0$ iff $C$ is a circle. Without
proof, we cite [37] the result that the isoperimetric deficit, defined as the
dimensionless quantity

$\displaystyle D(E)$ $\displaystyle=\frac{L}{\sqrt{4\pi A}}-1$ (S60)

upper-bounds any such quantity $\nu(E)$ via

$\displaystyle\nu(E)$ $\displaystyle\leq C_{2}\sqrt{D(E)}$ (S61)

for some dimension-dependent constant $C_{2}$. For convex $E$ (applicable e.g.
to the cells of a Voronoi tessellation of Euclidean space), (S60) upper-bounds
the Hausdorff distance to the best-fit equal-volume ball as

$\displaystyle\inf_{x\in\mathbb{R}^{2}}d_{H}(E,B+x)$ $\displaystyle\leq
C_{2}D(E)^{\alpha_{2}}$ (S62)

for dimension-dependent constants $C_{2},\alpha_{2}$ [37]. Naturally this is
sharp for $E$ which is a ball.

###  III.3 Curves on surfaces of constant Gaussian curvature

######  Identity S1 (Area-angle formula for a polygon on $\mathbb{S}^{2}$ or
$\mathbb{H}^{2}$).

Let $n$ denote the number of sides of a regular spherical or hyperbolic
polygon with interior angles $\theta_{i}$. Then by Gauss-Bonnet,

$\displaystyle KA+\sum_{i=1}^{n}(\pi-\theta_{i})$ $\displaystyle=2\pi$ (S63)

with $K$ the Gaussian curvature.

##  IV Models

Figure S1: Gamma-distributed stretches and approximately beta-prime
distributed aspect ratio.

Model IS.3 (variable $N$, fixed $C$)—As in Model IS.2, let us consider the
$k$th order statistic of $N$ uniform random variables (representing the fixed-
circumference constraint), with $N_{\lambda}$ now a Poisson-distributed random
variable conditioned to be minimum $k$. The marginal distribution of this
order statistic given $\lambda$ is (a particular) compound beta-Poisson
distribution, given by

$\displaystyle\mathrm{P}(Z_{k,N-k+1}=x|\lambda,N\geq k)$ $\displaystyle=$
$\displaystyle\sum_{n=k}^{\infty}\mathrm{P}(Z_{k,n-k+1}=x|N=n)\mathrm{P}(N=n|\lambda,N\geq
k).$ (S64)

To see the parametrization more clearly, let $m=n-k+1$; then,

$\displaystyle\eqref{eq:compound_order_statistic}$ $\displaystyle=$
$\displaystyle\frac{1}{1-\mathrm{P}(N<k)}\sum_{m=1}^{\infty}\frac{x^{k-1}(1-x)^{m-1}}{\mathrm{B}(k,m)}\frac{\lambda^{m-1+k}e^{-\lambda}}{\Gamma(m+k)}$
(S65) $\displaystyle=$
$\displaystyle\frac{1}{1-\mathrm{P}(N<k)}\frac{\lambda^{k}x^{k-1}e^{-\lambda}}{\Gamma(k)}\sum_{m=1}^{\infty}\frac{(\lambda(1-x))^{m-1}}{\Gamma(m)}$
(S66) $\displaystyle=$
$\displaystyle\frac{1}{1-\mathrm{P}(N<k)}\frac{\lambda^{k}x^{k-1}e^{-\lambda
x}}{\Gamma(k)}$ (S67) $\displaystyle=:$ $\displaystyle
c_{0}\frac{\lambda^{k}x^{k-1}e^{-\lambda x}}{\Gamma(k)},\qquad x\in[0,1].$
(S68)

Recognizing the second factor as the gamma distribution, the first factor
$c_{0}$ is simply a normalizing constant restricting the support to $[0,1]$.
Here we see the emergence of gamma distributions from the order statistics of
a Poisson-distributed number of uniform random variables, which is in fact a
roundabout way to construct the 1D Poisson process. As we take the support of
this distribution (the circumference $C$) to infinity, we have
$P(N_{C\lambda}<k)\to 0$ for any fixed $k$, hence $c_{0}\to 1$ and we recover
the true gamma distribution.

Model IIS.1 (periodic unit square)—Figure S1(b2-c2) shows that both major and
minor axes of the Poisson-Voronoi tessellation, nondimensionalized to
empirical mean 1, follow gamma distributions.

##  V Computational methods

###  V.1 Validation of the numerical method for Poisson-Voronoi tessellations

Let $V_{i}$ be the measure (area, volume, etc.) of the typical Poisson-Voronoi
cell. While the distribution of $V_{i}$ is presently unknown, exact second
moments of $V_{i}$ in $\mathbb{R}^{2}$ and $\mathbb{R}^{3}$ are known [26],
facilitating comparison with numerical study. Large simulations [27, 28] with
$n>10^{6}$ cells have found that gamma distributions, and in particular a
3-parameter generalization [27]

$f_{Y_{;k,\lambda,a}}(v)=\frac{a\lambda^{k/a}v^{k-1}e^{-\lambda
v^{a}}}{\Gamma(k/a)}$ (S69)

achieve good maximum-likelihood fit to data with $<1\%$ error relative to the
analytical second moment. In Fig. S2, the estimated second moment $\langle
V_{i}^{2}\rangle$, shape parameter $k$, and CDF root mean square error are
displayed for $500$ trials of $N\overset{\text{pdf}}{\ \sim\
}\text{Poisson}(10^{3})$ points. The average empirical, gamma, and beta second
moments show good agreement with Gilbert’s [26] numerically integrated value
of $1.280$ and are within $1\%$ relative error, validating the numerical
method. The estimated value of $k\approx 3.7$ for gamma on the torus
$\mathbb{T}$ is consistent with prior results finding $k\approx 3.6$ [29] in
the plane $\mathbb{R}^{2}$. On the other hand, the estimated value of
$k\approx 3.2$ for beta is lower than gamma and closer to Tanemura’s [27]
generalized-gamma (S69) fit finding $k=3.315$, suggesting that a beta
hypothesis is a good substitute for the generalized gamma distribution.
Lastly, we observe that the beta RMSE is slightly decreased compared to gamma.

Figure S2: Gamma- and Beta-maximum likelihood fits for Poisson-Voronoi
tessellations of the flat torus $\mathbb{T}^{2}$. The numerically integrated
value of the second moment $\mathbb{E}[V_{i}^{2}]=1.280$ found by [26] is
plotted as a dashed line in the left plot. The experimentally determined shape
parameter $k=3.315$ of the generalized-gamma fit found by [27] is plotted as a
dashed line in the middle plot.

###  V.2 Fitting ellipsoids in $d$ dimensions

Let $x,v$ be vectors in some $d$-dimensional basis. The equation of a
$(d-1)$-sphere $S$ centered at $v$ is

$\displaystyle\left\lVert x-v\right\rVert_{2}^{2}$
$\displaystyle=1\quad\forall x\in S.$ (S70)

Applying a rotation $P$ (an orthogonal matrix) of the sphere onto a set of
principal axes and stretching along those axes by $\Lambda$ (a positive
diagonal matrix), one generalizes (S70) to an ellipsoid $E$ via a symmetric
positive-definite matrix $M=P\Lambda P^{\intercal}$ such that $x$ satisfies

$\displaystyle(x-v)^{\intercal}M(x-v)$ $\displaystyle=1\quad\forall x\in E.$
(S71)

The elliptic radii $r_{i}$ are then $r_{i}=\Lambda_{ii}^{-\frac{1}{2}}$ and
the elliptic axes are the columns of $P$.

####  V.2.1 Minimum-volume bounding ellipsoid

Given the volume $V$

$V\propto\prod_{i=1}^{d}r_{i}=\sqrt{\det(M^{-1})},$ (S72)

it follows that maximizing $\log\det(M)$ minimizes $V$. Hence for a given
dataset $\\{x_{i}\\}_{i=1}^{n}$, the following convex program computes the
minimum-volume bounding ellipsoid:

$\displaystyle\sup_{M,v}\ $ $\displaystyle\log\det(M)$ (S73) subject to
$\displaystyle(x_{i}-v)^{\intercal}M(x-v_{i})\leq 1\quad\forall i$ (S74)
$\displaystyle M>0.$ (S75)

Here, $M>0$ is in the sense of linear matrix inequality (LMI), i.e. $M$ is
constrained to lie in the positive-definite cone. However, the offset
$x_{i}-v$ produces a variable-product constraint which is not disciplined-
convex. The following reparametrization uses the invertibility of $M$ to
eliminate the product from the problem:

$\displaystyle\sup_{A,b}\ $ $\displaystyle\log\det(A)$ (S76) subject to
$\displaystyle\left\lVert Ax_{i}-b\right\rVert_{2}^{2}\leq 1\quad\forall i$
(S77) $\displaystyle A$ $\displaystyle>0$ (S78) $\displaystyle M$
$\displaystyle=A^{2}=A^{\intercal}A$ (S79) $\displaystyle v$
$\displaystyle=A^{-1}b.$ (S80)

DCP solvers such as CVXOPT [38] solve this problem. We reduce the problem size
by considering the subset of $X$ lying on the convex hull, computed in
$O(n\log n)$.

####  V.2.2 Minimal $\ell^{2}$-distance projection to ellipsoid

Given a representation of an ellipsoid as $(M,v)$ in the same basis as $x$,
define the following convex program:

$\displaystyle\inf_{\hat{Y}}\ $
$\displaystyle\left\lVert\hat{Y}-\hat{X}\right\rVert_{2}^{2}$ (S81) subject to
$\displaystyle\left\lVert A\hat{x}_{i}\right\rVert_{2}^{2}=1\quad\forall i$
(S82) $\displaystyle A^{\intercal}A=A^{2}=M$ (S83) $\displaystyle\hat{X}=X-v$
(S84) $\displaystyle\hat{Y}=Y-v,$ (S85)

with $A$ computed either by (Hermitian) eigendecomposition or Cholesky
factorization of $M$. Then $Y$ is the minimum-$\ell^{2}$-distance projection
of $X$ onto $(M,v)$. This problem is not DCP; however, since there are no
matrix cone constraints, we can simply use non-DCP solvers compatible with
nonlinear constraints, such as SLSQP [39]. The constraint Jacobian for (S82)
is $2M\hat{x}_{i}$.

###  V.3 Fitting hyperplanes in $d$ dimensions

The equation of a hyperplane $H$ in $d$ dimensions is

$\displaystyle n\cdot(x-v)$ $\displaystyle=0\quad\forall x\in H$ (S86)

for $n,v\in\mathbb{R}^{d}$. Without loss of generality, we may assume that
$\left\lVert n\right\rVert=1$, so that $n$ is a unit normal to $H$; expressing
$H$ in the form

$\displaystyle n\cdot x-b$ $\displaystyle=0$ (S87)

we see that $b=n\cdot v$ is the distance from the origin to $H$. It further
follows that the distance from an arbitrary point $y\in\mathbb{R}^{d}$ to the
plane is

$\displaystyle d(y,H)$ $\displaystyle=\left|n\cdot(y-v)\right|=\left|n\cdot
y-b\right|.$ (S88)

Now, let $\\{x_{i}\\}_{i=1}^{N}=X\in\mathbb{R}^{N\times d}$ be a set of data
points with $N\geq 3$. The best-fit hyperplane (in the $\ell^{2}$ sense) is

$\displaystyle\inf_{n,v}\ $ $\displaystyle\
\left\lVert(X-v)n\right\rVert_{2}^{2}$ (S89) subject to $\displaystyle\
\left\lVert n\right\rVert_{2}=1.$ (S90)

The cost is bi-convex in the parameters $v$ and $n$. Calculating the gradient
of the cost with respect to $v$, we find:

$\displaystyle 0$ $\displaystyle=\frac{\partial}{\partial
v}\left\lVert(X-v)n\right\rVert_{2}^{2}=\sum_{i=1}^{N}(x_{i}-v)$ (S91)

Thus $v=\frac{1}{N}\sum_{i=1}^{N}x_{i}$ is the centroid. Letting
$\overline{X}=X-v$, the gradient of the cost with respect to $n$ is:

$\displaystyle 0$ $\displaystyle=\frac{\partial}{\partial
n}\left\lVert\overline{X}n\right\rVert_{2}^{2}=2\overline{X}^{\intercal}\overline{X}n$
(S92)

A common numerically stable method to estimate $n$ is the singular vector of
$\overline{X}$ corresponding to the smallest singular value, which is $0$ if
$\\{x_{i}\\}$ are coplanar. Moreover, by (S88) it follows that the
$\ell^{2}$-orthogonal projection of $X$ onto $H$ is

$\displaystyle Y$ $\displaystyle=X-(X-v)nn^{\intercal}.$ (S93)

Let $u\in\mathbb{R}^{d}$ be a random vector such that $u\times n\neq 0$ (e.g.
a Gaussian vector, for which this is almost surely true); then a planar
embedding $X_{H}$ of $X$ is defined by a random orthonormal basis:

$\displaystyle B$ $\displaystyle=\frac{1}{\left\lVert
v\right\rVert}\begin{bmatrix}v&v\times n\end{bmatrix},\quad v=u\times n,$
(S94) $\displaystyle X_{H}$ $\displaystyle=XB.$ (S95)

###  V.4 3D reconstruction & analysis of Volvox

Figure S3: Computational pipeline for extracting geometric features of V.
carteri. (a) The somatic cell positions, available from [5], displayed in
three-dimensional space. (b) The minimum-volume bounding ellipsoid containing
the somatic cells is computed as in Section V.2.1. This establishes the
approximate anterior-posterior axis of the spheroid. (c) The
$\ell^{2}$-orthogonal projection of the somatic cell positions $X$ to this
ellipsoid are computed as in Section V.2.2, and the Delaunay
tetrahedralization of this projected point cloud $\hat{X}$ is computed. As
$\hat{X}$ is now its own convex hull, the index triplets corresponding to the
triangular faces lying on the convex hull of the tetrahedral complex are taken
to be an approximate Delaunay triangulation of the original point cloud $X$.
(d) An approximate Voronoi tessellation of the surface is computed by sub-
triangulating the Delaunay triangulation, connecting the circumcenters of
adjacent triangular faces. The best-fit planes approximating these curved
faces are computed as in Section V.3, and the corresponding planar embeddings
are treated as closed plane curves for downstream analysis as in Section III.

First, 3D meshes and an ellipsoidal approximation of the organism’s surface
are constructed using the procedure detailed in Figure S3. Then, the anterior
section—typically free of reproductive cells (gonidia) which interrupt the
positions of somatic cells at the surface—is determined by separating cells
into two subsets lying on either side of the plane normal to the ellipsoid’s
major axis. The subset with smaller cell counts is used as a (generally good)
approximation of the anterior side. Finally, the areas of the Voronoi faces
are converted to solid-angles ($4\pi$ times the area fraction of total) and
are filtered by the cutoff $v_{c}=0.007$ specified in [5]. The empirical area
distribution for each organism is shifted by this cutoff and
nondimensionalized to empirical mean $1$, at which point they are combined
across all 6 organisms.