# Phenomenology and Dynamics of Competitive Ecosystems Beyond the Niche-
Neutral Regimes

Nava Leibovich  Jeremy Rothschild  Department of Physics, University of
Toronto  Sidhartha Goyal<EMAIL_ADDRESS>Anton Zilman
<EMAIL_ADDRESS>Department of Physics, University of Toronto
Institute for Biomedical Engineering, University of Toronto

###### Abstract

Structure, composition and stability of ecological populations are shaped by
the inter- and intra-species interactions within these communities. It remains
to be fully understood how the interplay of these interactions with other
factors, such as immigration, control the structure, diversity and the long
term stability of ecological systems in the presence of noise and
fluctuations. We address this problem using a minimal model of interacting
multi-species ecological communities that incorporates competition,
immigration and demographic noise. We find that the complete phase diagram
exhibits rich behavior with multiple regimes that go beyond the classical
‘niche’ and ‘neutral’ regimes, extending and modifying the ‘neutral-like’ or
‘niche-like’ dichotomy. In particular, we observe novel regimes that cannot be
characterized as either ‘niche’ or ‘neutral’ where a multimodal species
abundance distribution is observed. We characterize the transitions between
the different regimes and show how they arise from the underlying kinetics of
the species turnover, extinction and invasion. Our model serves as a minimal
null model of noisy competitive ecological systems, against which more complex
models that include factors such as mutations and environmental noise can be
compared.

ecological drift $|$ neutral-niche theory $|$ phase diagram $|$

##  I Introduction

Composition and behavior of ecological communities are shaped by direct and
indirect interactions between the occupant species, such as the competition
for the physical space and the intrinsic and the extrinsic resources. Examples
of such competitive ecosystems are microbial communities in various biomes
such as the soil [1], the ocean [2, 3] and the human body [4] \- in particular
the human gut which hosts a diverse microbiome whose dynamics are important
for human health [5, 6]. In the context of cellular populations within
organisms, the evolution of neoplasms and tumor cells [7, 8, 9], interactions
within the immune system [10, 11], as well as the appearance of dominant
clones during cell reprogramming [12], exhibit phenomenology akin to
ecological competition. Beyond biology [13, 14, 15], competitive interactions
shape behaviors in a vast array of systems such as competition economics [16]
and social networks [17].

A classical example of the effects of inter-species competition - which
inspired important ecological competition paradigms - is the differentiation
in beak forms of finches in the Galápagos islands [18, 19]. On these islands,
dissimilar finch species possess beaks of varying shapes and sizes allowing
them to consume different food sources and thus occupy distinct niches; this
type of ecosystem structure is commonly referred to as an ecological niche
model [20, 21]. Various niche models have been used to describe the community
structures observed in diverse ecosystems such as plant grassland communities
[22, 23], marine plankton [24] and conservation ecology [25, 26]. Commonly,
niche specialization results in weaker competition for resources between
individuals occupying separate niches (inter-species competition) compared to
the competition between individuals of the same kind residing in the same
niche (intra-species competition)[27, 28, 29].

Another paradigmatic class of ecological models that are used to describe
noisy ecosystems comprises neutral models. In contrast to niche models, in
neutral models individuals of all species are considered functionally
equivalent, and interactions between them are identical regardless of their
species [30, 31, 32, 33]. One classical example of a neutral model is the
Hubbell model, that showed that a neutral process underlying the population
dynamics of an ecosystem recovers experimental observations of species
abundances in tree communities [31]. Subsequently, neutral models have
commonly served as the paradigmatic null hypotheses for the exploration of
ecological processes in which the differences between inter-specific and
intra-specific interaction are functionally negligible [30, 34, 35, 36].
Neutral theories may be viewed as a limit of niche theories where inter-
specific and intra-specific interactions are equal: in other words, all
species reside in completely overlapping niches [37, 38, 21].

In multi-species communities, the intra- and inter-species interactions as
well as interactions with the environment, can lead to complex community
composition and population dynamics; some species survive in the long term,
while others are driven to extinction. However, in large communities with high
numbers of competing species, it is often impractical or impossible to
characterize the entire system composition by the assemblage of abundances for
each species. Hence, coarse-grained paradigmatic descriptions are often used
to provide general insights into the common behavior of these ecological
communities.

Two variables commonly used to characterize complex ecological communities are
1) the richness, reflecting the number of co-occurring species [39, 40], and
2) the species abundance distributions (SAD) - the number of species present
at a given abundance. The latter is closely related to the species
proportional abundance distribution (SPAD) as well as to the species rank
abundance (SRA) - the species ranked in terms of their abundance [33, 41, 42,
43, 44, 45]. These aggregate variables are observable experimentally and serve
as the reporters on the underlying community structure, dynamics and the
interaction network [46, 47, 48, 49]. Richness, for example, is commonly
considered to be an indicator of the competition strength and stability of the
ecosystem [50, 51, 52, 53, 54].

The shape of the SAD is also used as a proxy for the structure of the
underlying interactions’ network. For high immigration or weak inter-species
competition, the SAD commonly has a peak at high species abundance, away from
extinction. This community structure is closely related to the niche models
whereby different species co-exist: most species inhabit their own niches with
their species abundance fluctuating around the peak of the SAD. Conversely,
other ecosystems, such as many microbial communities and T-cell repertoires,
commonly comprise few high-abundance species alongside highly diverse
populations of low-abundance species [55, 43]. This unimodal, monotonically
decreasing SAD - sometimes called a ‘hollow-curved distribution’ - we refer to
as the ‘rare biosphere’ SAD. Interestingly, this behaviour is empirically
observed in many different ecosystems and is often considered universal (see
[56] and references therein). Neutral models have been championed to describe
the emergence of this universality, although other theoretical explanations
for the ‘rare biosphere’ SAD in competitive ecosystems have been suggested
[44, 57].

Theoretical studies commonly employ a small number of paradigmatic models to
quantify the competitive dynamics, the richness and the abundance
distributions in ecological populations. One common model of ecological
competition is the deterministic, competitive Lotka-Voltera (LV) model, which
has been especially useful in characterizing the niche regime by describing
stable species coexistence as stable fixed points of the model. Depending on
the ratios of inter- and intra-species competition strengths, deterministic LV
models provide examples of both the ‘niche-like’ regimes of multiple species
coexistence, and the competitive exclusion where species with weaker intra-
species interactions drive others to extinction [58, 59, 60, 61]. In complex
scenarios, such as when the strengths of inter-specific interactions are
randomly distributed among different species pairs, multi-species
deterministic LV models can exhibit not only deterministic fixed point
coexistence but also chaotic behavior reflected in the SAD shapes and richness
[62, 63, 64, 65, 66]. Beyond disorder in the interaction network, dynamical
noise from various sources - both extrinsic and intrinsic - has important
effects on the system composition and dynamics, especially in the neutral
regime. In order to capture experimentally observed stochastic fluctuations of
population abundances, environmental noise is often introduced into the
mathematical models [67, 55, 68, 69, 70]. In particular, by tuning the
strength of environmental noise the shape of the SAD can change from unimodal
to bimodal [67], indicating a transition between ‘niche-like’ and ‘rare
biosphere’ regimes. Incorporating both asymmetric interactions and
environmental noise, ‘patch models’ of communities have also been utilized to
study coexistence and abundances [71, 72] in island metacommunities.

Regardless of the presence of the external environmental noise or randomness
in the interaction network, the demographic noise - the inherent randomness of
birth and death events - is ever-present and has fundamental impact on the
community structure and stochastic population dynamics [33, 31, 73, 74]. In
particular, demographic noise in neutral systems has been shown to result in
an SAD shape characterized by a monotonically decreasing distribution often
referred to as a ‘rare biosphere’ distribution. Consequently, it has been
suggested that the ‘rare biosphere’ SAD observed in many experimental systems
is the outcome of neutral dynamics of ecological communities[33, 31, 75, 36].
On the other hand, neutral birth-death-immigration processes with demographic
noise have also been shown to exhibit bimodal SADs at very low immigration
rates [76] breaking from the paradigm wherein neutrality synonymously refers
to an SAD of the ‘rare biosphere’ type. Although demographic noise models have
been shown to reproduce the observed features of a number of ecological
systems [74, 28, 54, 77], a complete picture of the different regimes of
community structures, is still missing. In particular, it remains to be fully
understood how the interplay of the competition strength, the immigration
rate, demographic noise and the resulting dynamics of species turnover shape
transitions between these different community structure regimes.

A

B

Figure 1: Island model. Panel A: Conventionally, weak competition is
associated with ‘niche-like’ bimodal SAD, while strong competition is linked
to ‘rare biosphere’ monotonically decreasing SAD. However, this paradigm is
not complete, since the dependence on other parameters, such as immigration
rate $\mu$ or diversity $S$, is not fully investigated. Thus, the entire phase
space, e.g. $(\mu,\rho)$ or $(S,\rho)$, remains unexplored. Panel B: The model
illustration. An island with $J$ individuals from $S^{*}$ species. Each
individual may proliferate and die with some rate corresponding to inter- and
intraspecific interactions within the island. Here we consider deterministic,
symmetric, fully-connected interspecific interactions network, governed by
single parameter; the competition strength $\rho$. Additionally, individuals
may migrate from a cloud/mainland, contains $S$ species, into the island with
a constant rate $\mu$.

In this paper, we systematically investigate the full parameter space of the
community composition and structure using a competitive LV model with the
demographic noise and an interaction network of minimal complexity structure;
more complex scenarios may be examined by building on this paradigmatic null
model. We show that, beyond the perception of dichotomous neutral-niche
regimes, many different regimes of richness and SAD shape emerge from the
interplay between the competition strength and immigration in the presence of
stochasticity as illustrated in Fig. 1. These regimes exhibit contrasting
dynamics that underpin the differences in the community structures in
different regimes, and the transitions between them.

The paper is structured as follows. In Section II we introduce the minimal
model. In Section III we present our main results, including the regimes
boundaries, their richness and the abundance distributions, as well as their
associated underlying dynamics. Lastly, in Section IV, we discuss our results
in the context of experimental observations.

##  II Mathematical models and methods

The minimal model studied in this paper incorporates three essential features
of the ecological processes: competitive interactions, immigration and
intrinsic demographic noise [78, 74]. In the model, illustrated in Fig. 1B,
the community composition is characterized by the species abundances,
$\vec{n}=(n_{1},\dots n_{i}\dots n_{S})$ where the discrete random variable
$n_{i}$ represents the number of individuals of the $i$-th species, and $S$ is
the total number of species. The dynamics of the system are described by a
birth-death process with interactions, whereby the abundance (number of
individuals) of any species can increase by one with the birth rate $q^{+}$ or
decrease by one with the death rate $q^{-}$ defined as

$\displaystyle q_{i}^{+}(\vec{n})$ $\displaystyle=r^{+}n_{i}+\mu,$ (1)
$\displaystyle q_{i}^{-}(\vec{n})$
$\displaystyle=r^{-}n_{i}+\frac{r}{K}n_{i}\left(n_{i}+\sum_{j\neq
i}\rho_{j,i}n_{j}\right)$

for each species $i\in\\{1,2,\dots,S\\}$. These rates recover LV models which
have been extensively used to study deterministic multispecies coexistence,
and in particular to explore various stabilizing and equalizing mechanisms
[79, 70].

The birth rate incorporates two factors: the per-capita birth rate $r^{+}$
corresponding to procreation, and the constant and positive immigration rate
$\mu$ from an external basin which ensures that the system possesses no global
absorbing extinction state [28]. The death rates include the ‘bare’ per-capita
death rate of the organisms $r^{-}$ and the competitive interactions effects
that increase the mortality at high population numbers, incorporated through a
quadratic term in the death rates; Parameter $\rho_{j,i}$ quantifies the
competition strength between species $i$ and $j$; the competition strength
(analogous to the niche overlap [27, 28]) is defined as the ratio between the
inter-specific and intra-specific competition strengths. The carrying capacity
for each species is represented by $K$. The per-capita turnover rate is
$r=r^{+}-r^{-}$.

These aggregate coarse-grained parameters are determined by a variety of
system factors such as the efficiency of resource consumption, interactions
with the environment and external forces. Although it is possible to derive
these rates from explicit resource competition models in several special
cases, the expressions are highly model-dependent and are not explicitly
modeled here [80, 81, 82]. For biological reasons, $K$, $r^{+}$, $r^{-}>0$ are
all positive, which results in strictly positive transition rates for all
$n_{i}\geq 0$. In this paper, we focus on the homogeneous case where the
parameters ($\mu$, $K$, $\rho$, $r^{+}$, and $r^{-}$) are identical for all
species and the competitive interactions $\forall i,j:\rho_{j,i}=\rho$ for all
species pairs. This symmetric and homogeneous interaction network has been
used in [27, 54, 77, 74] in contrast to the models wherein the competition
strengths are inhomogeneous and drawn from a distribution [67, 83]. This
minimal complexity model allows us to investigate the full phase space of the
system to examine the underlying principle without extensive and impractical
multi-parameter sweeps.

The stochastic evolution of the system is described by the master equation

$\partial_{t}{\rm\mathcal{P}}(\vec{n};t)=\sum_{i}\left\\{-\left[q^{+}_{i}(\vec{n})+q^{-}_{i}(\vec{n})\right]{\rm\mathcal{P}}(\vec{n};t)\vphantom{\left[\sum
q^{+}_{i}\right]}\right.\\\ \left.\vphantom{\left[\sum
q^{+}_{i}\right]}+q^{+}_{i}(\vec{n}-\vec{e}_{i}){\rm\mathcal{P}}(\vec{n}-\vec{e}_{i};t)+q^{-}_{i}(\vec{n}+\vec{e}_{i}){\rm\mathcal{P}}(\vec{n}+\vec{e}_{i};t)\right\\},$
(2)

where $\vec{e}_{i}$ is the standard basis vector and
${\rm\mathcal{P}}(\vec{n},t)$ is the joint probability density function for
the system to exhibit the species composition $\vec{n}$ at time $t$ [84]. In
the long time limit, the system reaches a stationary state where
$\partial_{t}\mathcal{P}=0$, see SI Section 2 [85, 86, 87, 88].

The species abundance distribution (SAD) describing the mean fractions of
species with $n$ individuals, can be related to the marginal single species
probability distribution $P(n)$:

$\displaystyle{\rm SAD}$
$\displaystyle(n)=\frac{1}{S}\left\langle\sum_{i=1}^{S}\delta(n_{i}-n)\right\rangle$
(3)
$\displaystyle=\frac{1}{S}\sum_{i=1}^{S}\left[\sum_{n_{1}=0}^{\infty}\cdots\sum_{n_{i-1}=0}^{\infty}\sum_{n_{i+1}=0}^{\infty}\cdots\sum_{n_{S}=0}^{\infty}{\rm\mathcal{P}}(\vec{n})|_{n_{i}=n}\right]$
$\displaystyle=P_{i}(n)\equiv P(n),$

where $\delta$ is the Kronecker delta function, and using the fact that in
this homogeneous system the marginal distributions $P_{i}(n)=P(n)$ of
population abundance are identical for all species.

Dynamics of ecological populations can also be described using continuous
approximations (see the Supplementary Information; SI Section 1) [67, 55, 68,
69, 70, 33]. In particular, in the Fokker-Planck approximation, the continuous
deterministic limit of the master equation (Eq. 2) recovers the well-known
competitive Lotka-Volterra (LV) equations

$\displaystyle\frac{\partial x_{i}}{\partial t}$
$\displaystyle=q_{i}^{+}(\vec{x})-q_{i}^{-}(\vec{x})$
$\displaystyle=rx_{i}\left(1-\frac{x_{i}}{K}-\sum_{j\neq
i}\rho\frac{x_{j}}{K}\right)+\mu$ (4)

for the variable $x_{i}$, which corresponds to the continuous deterministic
limit of the discrete variable $n_{i}$ [84]; see SI Section 1 for further
details.

The deterministic steady state is given by

$\tilde{x}(S)=\frac{K}{2[1+\rho(S-1)]}\left\\{1+\sqrt{1+\frac{4\mu[1+\rho(S-1)]}{rK}}\right\\}.$
(5)

Note that in the deterministic LV process all species survive with abundance
$\tilde{x}$ as long as $\rho\leq 1$ and $\mu>0$ [28]. Conversely, in the
stochastic competitive environment the numbers of individuals of each species
fluctuate, occasionally reaching extinction. Thus, the number of co-existing
species $S^{*}$ is a stochastic variable as well, and may be smaller than the
overall number species $S$ in the immigration flux from the larger basin, with
$S^{*}\leq S$. The number of co-existing species has a corresponding
probability distribution whose evolution is governed by a master equation
derived from Eq. 2 (see SI Section 2) [89]. The richness, denoted as $\langle
S^{*}\rangle$, is defined as the average number of the (co-)existing species,
and is related to the SAD via

$\langle S^{*}\rangle=S(1-P(0)).$ (6)

Intuitively, this is the sum of the expectation of $S$ random indicator
variables; the richness is determined by $S$ times the probability that a
species is present in the system, $1-P(0)$ (see SI Section 3.A).

No exact analytical solution for the high-dimensional master equation Eq. 2 is
known for a general competition strength $\rho$. To understand the principles
of the community organization and the impact of competition, immigration and
demographic noise, we developed approximate analytical solutions to the master
equation verified by Gillespie simulations (see SI Section 2 for details).

##  III Results

###  III.1 Mean-Field Approximation

The full master equation Eq. 2 can be reduced to a one dimensional
approximation for the marginal distribution $P(n)$ with effective birth-death
rates (see SI Section 2.A). The SAD, $P(n)$, is obtained as a self-consistent
stationary solution of this equation as

$P(n)\equiv P_{i}(n_{i}=n)\\\
=P(0)\frac{(r^{+})^{n}(\mu/r^{+})_{n}}{n!\prod_{n_{i}=1}^{n}\left(r^{-}+rn_{i}/K+r\rho\sum_{j\neq
i}^{S}\langle n_{j}|n_{i}\rangle/K\right)}.\\\ $ (7)

To obtain an analytical approximation to $P(n)$ we use a mean field closure
for the unknown conditional averages $\langle n_{j}|n_{i}\rangle$ as
$\left\langle\sum_{j\neq i}n_{j}|n_{i}\right\rangle\approx(S-1)\langle
n\rangle$ (see SI Section 2.A for discussion and alternative approximations).
Thus, Eq. LABEL:eq:mean-field becomes a closed-form implicit equation for the
probability distribution $P(n)$ which can be solved numerically. We have found
a good agreement between exact stochastic simulation results and this mean-
field approximation for most of the parameter space examined.

Following Eq. 6, the average richness in the mean-field approximation is

$\displaystyle\langle
S^{*}\rangle=S\left(1-\frac{1}{{{}_{1}}F_{1}[a,b+1;c]}\right),$ (8)

where $P(0)=1/{{}_{1}F_{1}}[a,b+1+1;c]$ is the normalization constant of
$P(n)$ where ${{}_{1}F_{1}}[a,b;c]$ is the hypergeometric Kummer confluent
function, with $a=\mu/r^{+}$, ${b}=[r^{-}K+r\rho(S-1)\langle n\rangle]/r$, and
${c}={r^{+}K}/{r}$. The exact solution for the distribution of the number of
co-existing species, $S^{*}$, can be obtained for $\rho=0$ (a binomial
distribution) and $\rho=1$ (a sum of hypergeometric functions; see SI Section
3)[74]. For intermediate competition strengths, $0<\rho<1$, a mean-field
approximation results in a binomial distribution for the species richness
distribution as in [89]; however, we discuss how this mean-field solution
differs from the richness distribution from simulations in the SI Section 3.

###  III.2 The system exhibits rich behavior with distinct regimes of
population structures controlled by competition strength, immigration rate and
the species number

Depending on the values of the competition strength and the immigration rate,
the number of species and the system size, the population can exhibit a number
of different regimes of behavior which can be categorized by their richness
and the shape of their SAD, as visualized in Fig. 2 and described below.

####  III.2.1 Richness regimes

In the classical deterministic LV model, the systems exhibits either an
interior fixed-point with full coexistence of all species at abundances given
by Eq. 5, or mass extinction with a single surviving species, in agreement
with the well-known Gause’s law of deterministic competitive exclusion [28].
By contrast, the stochastic model may exhibit partial coexistence due to the
temporary extinctions of some species driven by the abundance fluctuations
arising from the demographic noise. Overall, the number of co-existing species
and their abundances are determined by the balance between the immigration and
the competition induced stochastic extinction events. Three distinct richness
regimes can be discerned as shown in Fig. 2, based on the variations of the
richness of the system $\langle S^{*}\rangle$ in different regions of the
($\rho$,$\mu$,$S$) parameter space.

At low competition strength - region (a) in Fig. 2A - all species co-exist so
that the richness of the system is equal to the total number of species
$\langle S^{*}\rangle\approx S$, similar to the deterministic regime. In this
regime, each species effectively inhabits its own niche because the inter-
species competition is not sufficiently strong to drive any of the species to
extinction in the presence of abundance fluctuations arising from the
demographic noise. The probability for a species to be present is determined
by the balance of its immigration rate and the extinction rate. At higher
immigration rates this regime extends into regions with higher competition
strength $\rho$: high immigration rates stabilize full richness populations
even with a relatively high competition strength.

In the second regime - region (b) in Fig. 2A - only a fraction of the species
are simultaneously present on average, which we denote as the partial
coexistence regime. In this regime, the immigration influx is not high enough
to prevent temporary stochastic extinctions of some species resulting from the
competition.

At very high competition strengths a complete exclusion regime - region (c) in
Fig. 2A - is found. High competition along with the very low immigration rates
act in unison to reduce the richness to below two species on average. Although
regime (c) may appear similar to regime (b) since both present partial
coexistence, they are distinguished by key behaviors as explained below.

Note that the stochasticity is central to the effect of the competition on the
observed richness. Stochastic fluctuations increase the risk of extinction
with increasing competition strength, unlike in the deterministic case where
the richness is independent of the competition strength for $\rho<1$ [28].

####  III.2.2 SAD shape and modality regimes

A

B

C

D

Figure 2: Phenomenology of the population structures. Panel A: The system
possesses three distinct richness phases. (a): full coexistence of all the
species $\langle S^{*}\rangle\approx S$; (b): partial coexistence with
$\langle S^{*}\rangle<S$; (c): a single species exists on average. Panel B:
Different population regimes are distinguished by different SAD modalities.
(I): immigration dominated regime with unimodal SAD at a typical abundance
given by the positive root of $\tilde{n}$; (II): bimodal regime with species
at non-zero abundance $\tilde{n}$ and a rapid species turnover peak a zero
abundance; (III): ‘rare biosphere’ regime of a unimodal SAD with peak at zero
abundance resulting from the rapid turnover of the temporarily extinct
species; (IV) multimodal regime. Panels C and D: Intersection of the modality
and richness regimes in the $(\mu,\rho)$ plane ( C) and $(S,\rho)$ plane (D);
see text for discussion. In panels  A, B and C the number of species $S=30$.
In panel D the immigration rate is $\mu=10^{-1}$. For all panels; Colored
regions represent data from simulation (see Methods), whereas boundaries from
the mean-field approximation are represented by solid black lines. The
solution for the master equation Eq. 2 is simulated using the Gillespie
algorithm with $6\cdot 10^{8}$ time steps, $r^{+}=2$, $r^{-}=1$, and $K=100$.

Besides determining the richness, the balance between immigration and
stochastic competitive extinctions also dictates the mean abundances of the
individual species and the species abundance distribution (SAD). When the
immigration influx of individuals into the system is higher than the average
out-flux due to the transient extinctions, shown in Fig. 2B as region (I),
most species are forced away from extinction. In this regime, the SAD is
unimodal with a peak at relatively high species abundances $\tilde{n}$
approximately located at

$\tilde{n}=\frac{K-\rho(S-1)\langle
n\rangle}{2}\left\\{1\pm\sqrt{1+4\frac{(\mu-r^{+})K}{r(K-\rho(S-1)\langle
n\rangle)^{2}}}\right\\},$ (9)

which agrees with the simulation results, as shown in Fig. 3; see also SI
Section 2.B.

At lower immigration rates - regime (II) in Fig. 2B - the immigration rate is
insufficiently strong to overcome the competition-driven temporary extinctions
of some species, and the SAD develops an additional peak around $n=0$
corresponding to the temporarily extinct species. The subset of the ‘quasi-
stable’ co-existing species dominate the population number with abundances
that fluctuate around the ‘niche-like’ abundance peak, $\tilde{n}$. Their
persistence at dominant abundances is punctuated by rare fluctuation-driven
extinctions and the occasional invasion of a temporarily extinct species into
the dominant population. By contrast, the dynamics of species in the $n=0$
zero peak is characterized by the rapid turnover of the remaining species
close to extinction. This balance between the immigration and the stochastic
competitive extinctions may be related to the trade-offs in competitive
ability and dispersal (immigration) in meta-community population dynamics [72,
90].

At low immigration rates, the peak at Eq. 9 coincides with the deterministic
stable solution in Eq. 5 (see SI Section 4.B)

$\lim_{\mu\rightarrow 0}\tilde{n}=\lim_{\mu\rightarrow
0}\tilde{x}\left(\langle S^{*}\rangle\right)=\frac{K}{1+\rho(\langle
S^{*}\rangle-1)}.$ (10)

Namely, in the bimodal regime the coexisting dominant species are fluctuating
around $\tilde{n}$ which, at low immigration, is the deterministic fixed point
with $\langle S^{*}\rangle$ species. In this regime, the dynamics of the
fluctuations of the abundant species around $\tilde{n}$ can be heuristically
understood as a spatially dependent diffusion in an effective potential well
of the Fokker-Plank Equation (See Section II and SI Section 1).

Somewhat unexpectedly, at low immigration rate $\mu\lesssim.05$, the bimodal
regime extends onto the neutral line at $\rho=1$ where the SAD has been
commonly believed to have the monotonically decreasing ‘rare biosphere’ shape
[31, 75]. Surprisingly, in this regime the competition is so strong that most
of the time either no species are present at high abundance, or only one
species survives in a kinetically ‘frozen’ and long lived quasi-stable state
with an abundance $\tilde{n}\simeq K$, as observed previously [76]; this is
region (IIc) in Fig. 1C.

Furthermore, at the intermediate immigration rates and relatively high
competition strengths we observe a unimodal behaviour with a peak at zero
rather than at a finite $\tilde{n}$ \- region (III) in Fig. 2B. In this
regime, the competition is strong enough so that the fluctuations
competitively drive species to temporary extinction before any of them is able
to establish a ‘quasi-stable’ state at a high abundance. All species undergo
rapid turnover around zero resulting from the balance between random
immigration and extinction events. This regime corresponds to what was
previously described as the ‘rare biosphere’: fewer number of species are
found at higher abundances resulting in a monotonically decreasing SAD. This
SAD shape is classically recognized as a hallmark of a ‘rare biosphere’
regime. However, as shown in Fig. 2 the unimodal regime (III) unexpectedly
extends substantially beyond the neutral manifold $\rho=1$, into the non-
neutral regions with $\rho<1$, and the monotonic-decreasing SAD persists even
for competition strengths as low as $\rho\approx 0.1$ \- an order of magnitude
weaker than the classical neutral regime. This challenges the common
perception that the ‘rare biosphere’ SAD is an indicator of neutrality. On the
neutral line $\rho=1$, for large $n$ the SAD asymptotically tends to a power
law with an exponential cutoff in line with similar functional forms found in
previous works (see SI Section 2.C)[31, 36, 75, 91].

Finally, we found an entirely novel multimodal regime with more than two peaks
- regime (IV) in Fig. 2 \- which possesses one rapid turnover peak around
extinction and multiple peaks at non-zero abundances. Similar to region (IIc),
the peak at $n=0$ comprises species which rapidly turnover around extinction.
However, in addition to the peak at positive abundances $K$ formed by one
surviving species ($S^{*}=1$) in a meta-stable frozen state, this regime
possesses a second peak at $\sim K/(1+\rho)$ with two simultaneously surviving
quasi-stable species ($S^{*}=2$. The abundance at these peaks are solutions to
Eq. 10 wherein $\langle S^{*}\rangle$ is replaced by the momentary $S^{*}$.
The slow fluctuations between the states with $S^{*}=1$ and $S^{*}=2$ result
in the appearance of the SAD with two non-zero modes at quasi-stable dominance
abundance, $\tilde{n}\sim K$ and $\tilde{n}\sim K/(1+\rho)$ observed in the
region (IV). These two peaks are only visibly separated when the richness is
low and carrying capacity is high, since solutions of Eq. 10 for different
$S^{*}$ are more distant in this regime.

The transitions between the different modality regimes and the corresponding
changes in the SAD shapes are illustrated in Fig. 3. Generally, at low
competition strength $\rho$ the species are practically independent of each
other, residing in largely non-overlapping niches and with their typical
abundance $\tilde{n}$ close to the carrying capacity $K$. Increasing
competition strength $\rho$ makes it harder to sustain the co-existing species
at high abundances, and accordingly $\tilde{n}$ decreases, as illustrated in
the top panels of Fig. 3A and Fig. 3B. With further increase in $\rho$ the
system behavior bifurcates depending on the immigration rate $\mu$. At high
immigration rates, $\mu\gtrsim 0.05$, the competition-driven decrease in
$\tilde{n}$ continues up to the critical competition strength (calculated in
the next section) where the peak around $\tilde{n}$ disappears (top right
panel of Fig. 3A) and Fig. 3B), as the system is not able to sustain ‘quasi-
stable’ niche-like species co-existence. This corresponds to the transition
from the bimodal region (II) to the ‘rare biosphere’ region (III) in Fig. 2).
At lower immigration rates (top left panel of Fig. 3A and Fig. 3B), further
increases in the competition strength eventually cause mass species
extinctions which allow the remaining few dominant species to maintain higher
abundances (region (III) Fig. 2). As $\rho\rightarrow 1$, the system
transitions to the region (IIc) of the Fig. 2: only one dominant species
remains, as described in [76], with abundance fluctuating around $K$.

####  III.2.3 Global Phase Diagram and Regime Boundaries

A

B

Figure 3: SAD changes between different regimes. Panel A: (upper left)
Simulation results for species abundance distributions (SADs) for fixed
$\mu=10^{-3}$ as a function of $\rho$. (upper right) same for $\mu=1$.
Different values of the competition strength $\rho$ are emphasized with
different colors indicated in the color-bar. (lower left) Simulation results
for SADs as a function of $\mu$ for fixed $\rho=0.5$ (lower right) same for
$\rho=1$. Different immigration rates $\mu$ are emphasised with different
color shown in the color-bar. Panel B: The non-zero mode of the SAD given by
the positive solution of $\tilde{n}$ representing the dominant species
abundance as a function of $\rho$ for different values of $\mu$. Markers and
dotted lines represent simulation results, while solid lines are given from
analytic analysis, Eq. 9.

In this section we describe the complete phase diagram of the system defined
by the intersection of the different richness and the SAD shape/modality
regimes, derive the regime boundaries and discuss the transitions between
them, as shown in the ($\mu,\rho$) space in Fig. 2C, and in ($S,\rho$) space
in Fig. 2D. We show that the boundaries between different regimes observed in
simulations can be understood within simple mean field theories, and discuss
the underlying physical factors responsible for the transitions between
different regimes.

We define the boundary between the full coexistence (a) and partial
coexistence (b) regimes to be at $\langle S^{*}\rangle=S-1/2$: the midpoint
between full richness $S^{*}=S$ and the loss of 1 species on average.
Similarly, the boundary between the partial coexistence (b) and exclusion (c)
regimes is located at $\langle S^{*}\rangle=3/2$, that is to say where the
richness is between one and two species such that on average only 1 species is
present in regime (c).

To derive the boundaries corresponding to the transitions of the SAD modality
regimes, we use discrete derivatives of the approximated SAD to determine the
existence of peaks and their location (see SI Section 4.C). The immigration
dominated regime (I) is characterized by a unimodal SAD with a peak at the
positive root of $\tilde{n}$ given in Eq. 9. Compared to this immigration
dominated regime, the neighboring bimodal and monotonically-decreasing
unimodal regimes - regions (II) and (III) respectively - differ by the
emergence of a new mode at zero abundance.

Thus, the boundary that defines transitions to either regime (II) or (III)
from the immigration dominated regime (I) is described by a flattening of SAD
at $n=0$: $\partial P(n)/\partial n|_{n=0}=0$. In the discrete case, this
heuristically corresponds to $P(0)=P(1)$. Combining this condition for the
boundary with the global-balance of the master equation Eq. 2 results in the
rate balance equation, $\langle q_{i}^{+}(\vec{n})|n_{i}=0\rangle=\langle
q_{i}^{-}(\vec{n})|n_{i}=1\rangle$.

In the mean-field approximation, this boundary is found at

$\mu=r^{-}+\frac{r}{K}[1+\rho(S-1)\langle n\rangle].$ (11)

This equation recovers the similar transition for $\rho=1$ derived
independently in [76].

The boundary between the bimodal regime (II) and the ‘rare biosphere’ regime
(III) is characterized by the disappearance of the peak at high abundance
$\tilde{n}$ in Eq. 9. In the bimodal regime at least one solution to
$\tilde{n}$ is real and positive; as such a maximal, real peak exists.
Conversely, in the ‘rare biosphere’ regime, both solutions of $\tilde{n}$ are
negative or imaginary. We find that the boundary between the real and
imaginary $\tilde{n}$ is

$r(K-\rho(S-1)\langle n\rangle)^{2}=4(r^{+}-\mu){K}$ (12)

and the transition line between positive and negative solutions,
$\tilde{n}=0$, is

$\frac{(K-\rho(S-1)\langle n\rangle)^{4}}{16}=1+\frac{K(\mu-r^{+})}{r}.$ (13)

The intersection of these two conditions defines the ‘rare biosphere’ regime
and is shown as the blue line in Fig. 2B,C.

The modality and the richness of the system are also affected by the number of
species $S$ as shown in Fig. 2D. In brief, the frequency of the immigration
events rises as more species are present in the immigration flux. Increased
immigration causes the total population to rise without providing more room
for each species in the system; this increases the stochastic competition,
driving more species to extinction. Hence, as $S$ increases, the transition
from the bimodal regime (II) to the unimodal regime (III) occurs at lower
values of competition strength $\rho$, and the fraction of the concurrently
surviving species decreases. This effect has been qualitatively observed
experimentally [92], and we return to it in the Discussion.

These analytical expressions for the regime boundaries - confirmed by
stochastic simulations - provide insights into the effects of different
control parameters on the regime boundaries. In particular, using the low
$\mu$ deterministic approximation for $\langle n\rangle\approx
K/\left[1+\rho(S-1)\right]$, shows that the location of the boundary of the
‘rare biosphere’ regime grows proportionally to the carrying capacity and is a
decreasing function of the number of species $S$. Thus, the size of the ‘rare
biopshere’ regime increases with the number of species $S$ as shown in Fig.
2D, whereas increasing the carrying capacity shrinks this regime (see SI
Section 5).

###  III.3 Kinetics of the species turnover, extinction and recovery underlie
the transitions between different regimes

To better understand the driving forces for the transitions between the
different regimes, we focused on the underlying kinetics of species turnover
and fluctuations. There is a stark contrast in the kinetics of an individual
species in the unimodal ‘rare biosphere’ regime (III) and the ‘niche-like’
regimes with a peak in SAD at non-zero abundance, as shown in Fig. 4A. In
regime (III), all species undergo rapid turnover in the relatively broad range
of abundances around extinction. By contrast, in regimes (I, II, and IV) the
‘quasi-stable’ dominant species undergo fast fluctuations around the co-
existence peak at $\tilde{n}$ in addition to fast turnover of the remaining
species near extinction. These fluctuations around the ‘quasi-stable’
abundance are punctuated by the temporary extinctions and the reverse
invasions of temporarily extinct species into the dominant ‘niche-like’ peak.

To characterize the kinetics in different regimes, we calculate the mean
first-passage times $T(a\rightarrow b)$ (MFPT) of the transitions between
different abundance levels ($a$ and $b$), using the one-dimensional backward
Master equation (see SI Section 6) [93, 94].

A

B

C

Figure 4: Kinetics of species extinction, invasion and turnover. Panel A:
Sample trajectories of the species abundances. (Upper panel): stable ‘niche-
like’ dynamics, where the dominant species fluctuate about $\tilde{n}$. The
red curve represents the corresponding bimodal SAD. (Lower panel): the erratic
dynamics in the ‘rare biosphere’ regime, where species fluctuate close to
extinction. The SAD is a monotonically decreasing function. Panel B: The MFPT
ratio $T(\tilde{x}\rightarrow 0)/T(\tilde{x}\rightarrow\tilde{x})$ as a
function of $\mu$ and $\rho$. Select contour lines are highlighted as dashed
lines. This ratio qualitatively captures the transition from ‘rare biosphere’
to ‘’niche-like’ regimes Fig. 2A. For weak immigration rates $\mu\approx
10^{-3}$ the ratio is non-monotonic in the competition strength, revealing
regime (c). Panel C: The MFPT ratio $T(0\rightarrow\tilde{x})/T(0\rightarrow
0)$. This ratio qualitatively captures the richness behaviour.

We first focus on the ratio of the MFPT of the transition from dominance to
exclusion to the MFPT of return to the dominant abundance level (starting from
the dominant abundance level), $T(\tilde{x}(\langle S^{*}\rangle)\rightarrow
0)/T(\tilde{x}(\langle S^{*}\rangle)\rightarrow\tilde{x}(\langle
S^{*}\rangle))$, shown in (Fig. 4B). Here, $\tilde{x}$, given in Eq. 5 is the
deterministic extension of the peak abundance $\tilde{n}$ in regimes without a
non-zero abundance peak. Large values of this ratio signify that the
extinction rate from $\tilde{x}(\langle S^{*}\rangle)$ is much slower than the
rate of local fluctuations in the effective potential well around
$\tilde{x}(\langle S^{*}\rangle)$. Accordingly, Fig. 4B shows that this ratio
is high in the bimodal and immigration-dominated regimes. Conversely, this
ratio is lower within the ‘rare biosphere’ regime which does not possess a
high abundance peak with ‘quasi-stable’ co-existing species. As shown in Fig.
4B, this ratio approximately delineates the ‘rare biosphere’ regime from the
‘niche-like’ regimes and its contour lines qualitatively recover the
boundaries of region IIIb in Fig. 2C; see SI Section 6 for further discussion.

The second ratio, which underlies the richness transitions in the system,
$T(0\rightarrow\tilde{x}(\langle S^{*}\rangle))/T(0\rightarrow 0)$ (Fig. 4
panel C) relates MFPT from extinction at zero abundance to dominance at
$\tilde{x}$ to the the mean return time to extinction. This ratio gives a
rough measure of the number of species present in the system: $T(0\rightarrow
0)$ is inversely proportional to $P(0)$ and $T(0\rightarrow\tilde{x}(\langle
S^{*}\rangle))$ is heuristically inversely proportional to the number of co-
existing species (see SI Sections 3A and 6). As such, this MFPT ratio
approximates the ratio of the average number of temporarily extinct species,
$S-\langle S^{*}\rangle$ to the average number of existing species, $\langle
S^{*}\rangle$, see Fig. 4C. As shown in Fig. 4C ratio quantitatively recovers
the boundaries of richness regimes in Fig. 2 in most regions of the parameter
space.

These MFPTs not only serve to interpret the underlying dynamics that give rise
to different regimes, they may also be more easily experimentally measured
than steady state distributions. Further discussion on the dynamical features
are presented in the SI Section 6.

##  IV Discussion

Ecological systems display a wide variety of different behavior regimes that
have been commonly analysed through a limited number of paradigmatic models
such as the ‘niche’ and ‘neutral’ theories. However, it remains incompletely
understood what features of ecological population structure and dynamics are
universal and which are system specific, how different models relate to each
other, and what behavior is expected in the full range of the parameter space.
Using a minimal model of the competitive population dynamics with demographic
noise, we have investigated the different regimes of the population structures
and dynamics as a function of the immigration rate $\mu$, the competition
strength $\rho$, as well as the number of species $S$. Although this minimal
model may not fully capture the more complex interaction structures of many
ecological communities, it exhibits rich and unexpected behaviours paralleling
many experimental observations (see Table 1), and illuminates the underlying
mechanisms that shape population structures in different ecosystems.

We have focused on the system richness reflecting the number of the co-
existing species, and the SAD shape as the characteristics of the different
population regimes, using a combination of simulations and analytical mean-
field approaches. Our analysis shows that the ecosystem behaviors can be
partitioned into different regimes of richness and SAD shape/modality,
parameterized by the immigration rate and the competition strength - which is
analogous to the competitive overlap discussed in other studies and references
therein [77, 27].

Our model recovers the limits of the well known ‘rare biosphere’ and the
‘niche-like’ regimes. In particular, at $\rho=1$ and intermediate values of
$\mu$, the SAD has the monotonically decreasing shape characteristic of the
classical neutral regime. On the other hand, at low competition strength, the
system SAD exhibits a peak at high species abundance where all species co-
exist. We recover the expected regime where different species effectively
occupy distinct ecological niches. Notably, even independent species with no
inter-species competition with $\rho=0$ may present either a unimodal or
bimodal SAD depending on the immigration rate, as shown in Fig. 2B, C. Unlike
the immigration dominated high abundance peak at high immigration rates, at
the very low immigration rates the SAD is peaked around zero due to high
extinction probability solely from the intra-species competition.

We found that, contrary to the common expectation that different species
inhabit separate niches away from neutrality, the system can maintain the
monotonically decaying ‘rare biosphere’ SAD even at low competition strength
(up to $\rho\approx 0.1$) as shown in the regime (III) in Fig. 2. Similarly,
unexpectedly, at the very low immigration rates, the system SAD maintains the
peak at non-zero abundance characteristic of ‘niche-like’ regimes even for the
high values of the competition strength $\rho$ usually considered to be in the
‘rare biosphere’ domain (regime (IIc)) in Fig. 2 and Section IIIB.

We have also uncovered an unusual - and to the best of our knowledge hitherto
not described - regime characterized by the multi-modal SAD with more than one
positive, ‘quasi-stable’ abundance peak (Regime (IV) in Fig 2). This multi-
modality arises from the richness fluctuations in this regime: the number of
co-existing species is switching randomly between two relatively long-lasting
states with $S^{*}=1$ and $S^{*}=2$. Thus, one peak of the SAD is found around
$\sim K$ and the other one in the vicinity of $\sim K/2$, as explained in
Section IIIB. We observe that for low $K$, the multimodal regime is non-
existent and appears as $K$ increases; see the corresponding phase diagrams in
SI Section 5.

Transitions of the SAD between different regimes occur through different
routes. In particular, as the immigration rate increases, the SAD peak of the
bimodal regime (II) at $\rho=1$, gradually decreases in height without
significantly shifting its location until it disappears at the boundary of the
‘rare biosphere’ regime (III). By contrast, at lower competition strengths
$\rho<1$, the transition from the bi-modality to the ‘rare biosphere’ regime
occurs via simultaneous changes in the peak’s height and location. This is
discussed in Section IIIB.

We show that the population structures in different regimes and the
transitions between them, stem from the underlying dynamics of species
fluctuations, extinctions and invasions. In the ‘rare biosphere’ regimes, all
species undergo relatively fast turnover around extinction. This is reflected
in the low ratio of the turnover to the extinction mean first-passage times.
Conversely, in the ‘niche-like’ regimes the system develops two additional
time scales: relatively fast fluctuations about the high abundance peak, and
the long waiting times for the transitions from the ‘quasi-stable’ co-
existence at high abundance to extinction. This is reflected in the fact that
the ratio of the mean extinction time to the mean time of return to dominance
is higher in the ‘niche-like’ regime, as discussed in Section IIIC.

Interestingly, ecological regimes akin to those predicted by our demographic
noise model (except for the multimodal SAD regime) have been also found using
deterministic, noiseless LV models with a random matrix of inter-species
competitive competition strengths [95, 96, 83, 64, 92]. However, the
underlying mechanisms that give rise to the apparently similar regimes in the
two model types are very different. In the demographic noise model, the
partial richness ‘niche-like’ regime (IIb) (Fig. 2C) comprises the quasi-
stable coexistence of a subset of species at a positive abundance in parallel
with the temporary stochastic extinctions of other species. By contrast, in
the deterministic LV models with random asymmetric interactions, the partial
richnes ‘niche-like’ regime comprises large number of saddle fixed points
where different sets of species are competitively excluded deterministically.
At higher competition strengths, the deterministic system transitions to the
chaotic behavior that resembles the ‘rare biosphere’ regime (III) (Fig. 2C),
however the nature of the species turnover and the shape of the SAD are
different from the results we presented in Section III [97, 64, 92].

The existence of the predicted regimes and the transitions between them can be
tested experimentally by measuring the SAD and the dynamics of the species
abundances in ecosystems with varying immigration and competition strengths,
numbers of species and effective carrying capacities. Measurements of the SADs
and the community compositions have become more attainable due to the advances
in single cell gene sequencing techniques [1, 92, 12], overcoming the
difficulties of SAD estimation due to data limitations. Long-term observations
may provide measurements of the stationary species abundance distributions
[98]. Although it may be difficult to experimentally determine and control the
immigration rate, the competition strength, and the carrying capacity,
practical proxies for these parameters exist. By way of an example, the flow
rate carrying bacteria into a chamber of a microfluidic device is a well
controlled quantity that approximates well the immigration rate for
populations encased in the chamber [99]. Another commonly used and robustly
estimated experimental observable is the species rank abundance (SRA), which
can be used to infer the SAD to which it is closely mathematically related
(see SI Section 7), although in practice the conversion might be constrained
by limitations of noise and quantity of the experimental data.

The asymptotic behaviour of the SADs may show qualitative dissimilarities
between distributions allowing one to discern different regimes of behavior
among the experimental observations. In the mean-field approximation, the
asymptotic behaviour of the model’s SAD on the neutral line $\rho=1$ is well
approximated by a power law with an exponential cutoff (see SI Section 2.C).
This asymptotic is similar in functional form to the SADs commonly found by
Hubbell models of a neutral birth-death process with a fixed total population
size [75, 36]. Notably, the Yule process that is often used to model neutral
processes also results in the SAD of a similar form. However, the Yule process
is substantially different from the model of this paper because it does not
include inter-species interactions and reaches the steady state SAD only if
the effective death rate is higher than the effective birth/immigration rate
[100].

In Table 1, we qualitatively compare the family of the regimes predicted by
our model to the various behaviors inferred from experimental findings based
on the SAD measurements and population abundance time series. The apparent
abundance of the neutral ecosystems observed experimentally - such as
gastrointestinal microbiomes - may pertain to our finding (Section IIIB) that
the ‘rare biosphere’ regime extends substantially beyond the neutral line
$\rho=1$: non-neutral communities may appear neutral as they exhibit SAD’s
characteristic of neutral communities [101]. Furthermore, multimodal SAD’s
predicted by our model that are related to the richness fluctuations may
provide an explanation for the multimodal SADs observed in some ecological
data, complementary to the existing explanations such as spatial heterogeneity
or emergent neutrality [102, 63]

System (Ref.) |  Regimes |  Observations
---|---|---
microbial competition [92] |  stable full coexistence (IIa), stable partial coexistence (IIb), persistent fluctuation (IIIb) |  Community composition/ richness/ fluctuating communities
global birds species [103] |  unimodal - log skew (I) |  SAD
plankton [104] |  power-law decay (III) |  SAD and SRA
coral [102] |  multimodal (IV) |  SAD
arthropods [105] |  multimodal (IV) |  SAD
T-cell receptors [106] |  bimodal (II) and unimodal (III) |  SAD
microbial competition [107] |  ‘rare biosphere’ (III) and ‘niche-like’ (I & II) |  SRA and time series
gastrointestinal microbiomes [101] |  ‘rare biosphere’ (III) |  SRA and operational taxonomic units (OTUs)
Table 1: Qualitative classification of observed population regimes in various
ecological systems.

One quantity that is relatively easy to control experimentally is the total
number of species $S$. The regimes predicted by the model and the transitions
between them are shown in Fig. 2D: our model yields ‘rare biosphere’ regimes
for high $S$ and $\rho$, which are characterized by high-turnover dynamics,
and ‘niche-like’ regines with more stable behavior at lower $S$ or $\rho$.
These predictions qualitatively agree with the experimentally observed phase-
space in [92], which were previously explained within the deterministic LV
models with a random interaction matrix [92]. The fact that both the
deterministic LV model with a random interaction matrix and the homogeneous LV
model with demographic noise are in qualitative agreement with the
experimental data raises interesting and important questions concerning the
interplay of stochastic and deterministic dynamics in determining the
community composition.

Another quantity that may enable qualitative and quantitative testing of
different models is the carrying capacity $K$ which may be controlled
experimentally in some systems. As shown in SM, the ‘rare biosphere’ regime
shrinks in size with increasing $K$ because a higher carrying capacity can
sustain higher average abundance, and larger (less likely) fluctuations are
needed for the extinction events to occur. Higher average abundance together
with insufficiently strong fluctuations result in longer MFPTs from dominance
to extinction abundances and vice-versa. These effects will be investigated in
future work

In the context of other ecological theories, the competition strength, as
defined in this work, can be viewed as a quantification of the heuristic
notion of the niche overlap, and we observe that decreasing niche overlap
results in richness increases as suggested previously [108, 109]. Our model
also serves as a quantitative example of some of the coexistence promoting
mechanisms of the contemporary ecological theory; we explore stabilizing
mechanisms that increase richness via decreases in niche overlap, such as
fluctuation-dependent processes and fitness-density covariance [29, 110]. In
particular, the demographic noise model studied here exhibits fluctuation-
dependent mechanisms that promote richness as species are able to coexist at
high abundance in our model.

We expect that the minimal model of this paper can be used for more
complicated scenarios, including more complex distributions of the interaction
network $\rho_{i,j}$, speciation to probe the interaction of the natural
selection, and inter-species interactions and population diversity and
structure.

Finally, our model of a local island community in the mainland-island
ecosystem (see Fig. 1), can be expanded to many-island models or many-patch
dynamics [71, 72]. These many-island and many-patch models examine the
interplay between competition and dispersal rate [72] and its effects on the
diversity of the metacommunity, a prominent topic in conservation ecology and
the study of the human microbiome. The patch models may address how
coexistence and persistence are influenced by spatial heterogeneity and
environmental noise in a demographic noise formulation. Future work will
explore integration of our model into other scenarios to predict species
fitness, non-equilibrium coexistence and their connection to broader
qualitative ideas in ecology.

## Methods

The solution for the master equation Eq. 2 is simulated using the Gillespie
algorithm with $10^{8}$ time steps. We use $r^{+}=2$, $r^{-}=1$, $K=100$.
Modalities’ classification is numerically executed after smoothing the
simulated SAD. The MFPT is evaluated via the simulated SAD ($\tilde{x}(S^{*})$
is rounded), where a uni-dimensional approximation of the process is
considered, see details in SI Section 6.

###### Acknowledgements.

The authors acknowledge helpful discussions and comments from all the members
of the Goyal and Zilman Groups. AZ acknowledges the support from the National
Science and Engineering Research Council of Canada (NSERC) through the
Discovery Grant Program. SG acknowledges the support from the National Science
and Engineering Research Council of Canada (NSERC) through the Discovery Grant
Program and from the Medicine by DEsign Program at the University of Toronto.

## References

  * Ratzke  _et al._ [2020] C. Ratzke, J. Barrere, and J. Gore, Strength of species interactions determines biodiversity and stability in microbial communities, Nature ecology & evolution 4, 376 (2020).
  * Tilman [1977] D. Tilman, Resource competition between plankton algae: an experimental and theoretical approach, Ecology 58, 338 (1977).
  * Strom [2008] S. L. Strom, Microbial ecology of ocean biogeochemistry: a community perspective, Science 320, 1043 (2008).
  * Foster  _et al._ [2017] K. R. Foster, J. Schluter, K. Z. Coyte, and S. Rakoff-Nahoum, The evolution of the host microbiome as an ecosystem on a leash, Nature 548, 43 (2017).
  * Coyte  _et al._ [2015] K. Z. Coyte, J. Schluter, and K. R. Foster, The ecology of the microbiome: networks, competition, and stability, Science 350, 663 (2015).
  * Gorter  _et al._ [2020] F. A. Gorter, M. Manhart, and M. Ackermann, Understanding the evolution of interspecies interactions in microbial communities, Philosophical Transactions of the Royal Society B 375, 20190256 (2020).
  * Merlo  _et al._ [2006] L. M. Merlo, J. W. Pepper, B. J. Reid, and C. C. Maley, Cancer as an evolutionary and ecological process, Nature reviews cancer 6, 924 (2006).
  * Kareva [2015] I. Kareva, Cancer ecology: Niche construction, keystone species, ecological succession, and ergodic theory, Biological Theory 10, 283 (2015).
  * Smart  _et al._ [2021] M. Smart, S. Goyal, and A. Zilman, Roles of phenotypic heterogeneity and microenvironment feedback in early tumor development, Physical Review E 103, 032407 (2021).
  * Tauber [2008] A. I. Tauber, The immune system and its ecology, Philosophy of Science 75, 224 (2008).
  * Schmid-Hempel and Ebert [2003] P. Schmid-Hempel and D. Ebert, On the evolutionary ecology of specific immune defence, Trends in Ecology & Evolution 18, 27 (2003).
  * Shakiba  _et al._ [2019] N. Shakiba, A. Fahmy, G. Jayakumaran, S. McGibbon, L. David, D. Trcka, J. Elbaz, M. C. Puri, A. Nagy, D. van der Kooy, _et al._ , Cell competition during reprogramming gives rise to dominant clones, Science 364 (2019).
  * Tilman [1982] D. Tilman,  _Resource competition and community structure_ (Princeton university press, 1982).
  * Morin [2009] P. J. Morin,  _Community ecology_ (John Wiley & Sons, 2009).
  * Tuljapurkar [2013] S. Tuljapurkar,  _Population dynamics in variable environments_ , Vol. 85 (Springer Science & Business Media, 2013).
  * Budzinski [2007] O. Budzinski, Monoculture versus diversity in competition economics, Cambridge Journal of Economics 32, 295 (2007).
  * Koura  _et al._ [2017] Y. H. Koura, Y. Zhang, and H. Liu, Competitive interaction model for online social networks’ users’ data forwarding at a subnet, Mathematical Problems in Engineering 2017 (2017).
  * Lewin [1983] R. Lewin, Finches show competition in ecology, Science 219, 1411 (1983).
  * Lack [1983] D. Lack,  _Darwin’s finches_ (Cambridge University Press, 1983).
  * Grant and Grant [1979] B. Grant and P. Grant, Darwin’s finches: population variation and sympatric speciation, Proceedings of the National Academy of Sciences 76, 2359 (1979).
  * Pocheville [2015] A. Pocheville, The ecological niche: history and recent controversies, in  _Handbook of evolutionary thinking in the Sciences_ (Springer, 2015) pp. 547–586.
  * Zuppinger-Dingley  _et al._ [2014] D. Zuppinger-Dingley, B. Schmid, J. S. Petermann, V. Yadav, G. B. De Deyn, and D. F. Flynn, Selection for niche differentiation in plant communities increases biodiversity effects, Nature 515, 108 (2014).
  * Silvertown [2004] J. Silvertown, Plant coexistence and the niche, Trends in Ecology & evolution 19, 605 (2004).
  * Cullen and MacIntyre [1998] J. J. Cullen and J. G. MacIntyre, Behavior, physiology and the niche of depth-regulating phytoplankton, Nato Asi Series G Ecological Sciences 41, 559 (1998).
  * Melo-Merino  _et al._ [2020] S. M. Melo-Merino, H. Reyes-Bonilla, and A. Lira-Noriega, Ecological niche models and species distribution models in marine environments: A literature review and spatial analysis of evidence, Ecological Modelling 415, 108837 (2020).
  * Aguirre-Gutiérrez  _et al._ [2015] J. Aguirre-Gutiérrez, H. M. Serna-Chavez, A. R. Villalobos-Arambula, J. A. Perez de la Rosa, and N. Raes, Similar but not equivalent: ecological niche comparison across closely–related m exican white pines, Diversity and distributions 21, 245 (2015).
  * Badali and Zilman [2020] M. Badali and A. Zilman, Effects of niche overlap on coexistence, fixation and invasion in a population of two interacting species, Royal Society open Science 7, 192181 (2020).
  * Capitán  _et al._ [2015] J. A. Capitán, S. Cuenda, and D. Alonso, How similar can co-occurring species be in the presence of competition and ecological drift?, Journal of the Royal Society Interface 12, 20150604 (2015).
  * Chesson [2000] P. Chesson, Mechanisms of maintenance of species diversity, Annual review of Ecology and Systematics 31, 343 (2000).
  * Bell [2001] G. Bell, Neutral macroecology, Science 293, 2413 (2001).
  * Hubbell [2001] S. P. Hubbell,  _The unified neutral theory of biodiversity and biogeography (MPB-32)_ , Vol. 32 (Princeton University Press, 2001).
  * Chave [2004] J. Chave, Neutral theory and community ecology, Ecology letters 7, 241 (2004).
  * Marquet  _et al._ [2017] P. A. Marquet, G. Espinoza, S. R. Abades, A. Ganz, and R. Rebolledo, On the proportional abundance of species: Integrating population genetics and community ecology, Scientific reports 7, 1 (2017).
  * Gotelli and McGill [2006] N. J. Gotelli and B. J. McGill, Null versus neutral models: what’s the difference?, Ecography 29, 793 (2006).
  * Blythe [2012] R. A. Blythe, Neutral evolution: a null model for language dynamics, Advances in complex systems 15, 1150015 (2012).
  * McKane  _et al._ [2004] A. J. McKane, D. Alonso, and R. V. Solé, Analytic solution of hubbell’s model of local community dynamics, Theoretical Population Biology 65, 67 (2004).
  * Grover  _et al._ [1997] J. P. Grover, J. Hudziak, and J. D. Grover,  _Resource competition_ , Vol. 19 (Springer Science & Business Media, 1997).
  * Begon  _et al._ [2006] M. Begon, C. R. Townsend, and J. L. Harper,  _Ecology: from individuals to ecosystems_ , Sirsi) i9781405111171 (2006).
  * Adams [2009] J. Adams,  _Species richness: patterns in the diversity of life_ (Springer, 2009).
  * Kéry and Royle [2020] M. Kéry and J. A. Royle,  _Applied Hierarchical Modeling in Ecology: Analysis of distribution, abundance and species richness in R and BUGS: Volume 2: Dynamic and Advanced Models_ (Academic Press, 2020).
  * Nias [1968] A. Nias, Clone size analysis: a parameter in the study of cell population kinetics, Cell Proliferation 1, 153 (1968).
  * Rulands  _et al._ [2018] S. Rulands, F. Lescroart, S. Chabab, C. J. Hindley, N. Prior, M. K. Sznurkowska, M. Huch, A. Philpott, C. Blanpain, and B. D. Simons, Universality of clone dynamics during tissue development, Nature physics 14, 469 (2018).
  * de Greef  _et al._ [2020] P. C. de Greef, T. Oakes, B. Gerritsen, M. Ismail, J. M. Heather, R. Hermsen, B. Chain, and R. J. de Boer, The naive t-cell receptor repertoire has an extremely broad distribution of clone sizes, Elife 9, e49900 (2020).
  * McGill  _et al._ [2007] B. J. McGill, R. S. Etienne, J. S. Gray, D. Alonso, M. J. Anderson, H. K. Benecha, M. Dornelas, B. J. Enquist, J. L. Green, F. He, _et al._ , Species abundance distributions: moving beyond single prediction theories to integration within an ecological framework, Ecology letters 10, 995 (2007).
  * Matthews and Whittaker [2015] T. J. Matthews and R. J. Whittaker, On the species abundance distribution in applied ecology and biodiversity management, Journal of Applied Ecology 52, 443 (2015).
  * Rahbek and Graves [2001] C. Rahbek and G. R. Graves, Multiscale assessment of patterns of avian species richness, Proceedings of the National Academy of Sciences 98, 4534 (2001).
  * Hong  _et al._ [2006] S.-H. Hong, J. Bunge, S.-O. Jeon, and S. S. Epstein, Predicting microbial species richness, Proceedings of the National Academy of Sciences 103, 117 (2006).
  * Adler  _et al._ [2011] P. B. Adler, E. W. Seabloom, E. T. Borer, H. Hillebrand, Y. Hautier, A. Hector, W. S. Harpole, L. R. O’Halloran, J. B. Grace, T. M. Anderson, _et al._ , Productivity is a poor predictor of plant species richness, Science 333, 1750 (2011).
  * Valencia  _et al._ [2020] E. Valencia, F. de Bello, T. Galland, P. B. Adler, J. Lepš, E. Anna, R. van Klink, C. P. Carmona, J. Danihelka, J. Dengler, _et al._ , Synchrony matters more than species richness in plant community stability at a global scale, Proceedings of the National Academy of Sciences 117, 24345 (2020).
  * Pimm [1984] S. L. Pimm, The complexity and stability of ecosystems, Nature 307, 321 (1984).
  * Ives  _et al._ [2000] A. R. Ives, J. L. Klug, and K. Gross, Stability and species richness in complex communities, Ecology Letters 3, 399 (2000).
  * Jousset  _et al._ [2011] A. Jousset, W. Schulz, S. Scheu, and N. Eisenhauer, Intraspecific genotypic richness and relatedness predict the invasibility of microbial communities, The ISME journal 5, 1108 (2011).
  * Mallon  _et al._ [2015] C. A. Mallon, J. D. Van Elsas, and J. F. Salles, Microbial invasions: the process, patterns, and mechanisms, Trends in microbiology 23, 719 (2015).
  * Capitán  _et al._ [2017] J. A. Capitán, S. Cuenda, and D. Alonso, Stochastic competitive exclusion leads to a cascade of species extinctions, Journal of Theoretical Biology 419, 137 (2017).
  * Lynch and Neufeld [2015] M. D. Lynch and J. D. Neufeld, Ecology and exploration of the rare biosphere, Nature Reviews Microbiology 13, 217 (2015).
  * Leidinger and Cabral [2017] L. Leidinger and J. S. Cabral, Biodiversity dynamics on islands: Explicitly accounting for causality in mechanistic models, Diversity 9, 30 (2017).
  * Magurran [2013] A. E. Magurran,  _Measuring biological diversity_ (John Wiley & Sons, 2013).
  * Hardin [1960] G. Hardin, The competitive exclusion principle, Science 131, 1292 (1960).
  * MacArthur and Levins [1967] R. MacArthur and R. Levins, The limiting similarity, convergence, and divergence of coexisting species, The american naturalist 101, 377 (1967).
  * Mac Arthur [1969] R. Mac Arthur, Species packing, and what competition minimizes, Proceedings of the National Academy of Sciences 64, 1369 (1969).
  * Gause [2019] G. F. Gause,  _The Struggle for Existence: A Classic of Mathematical Biology and Ecology_ (Courier Dover Publications, 2019).
  * Scheffer and van Nes [2006] M. Scheffer and E. H. van Nes, Self-organized similarity, the evolutionary emergence of groups of similar species, Proceedings of the National Academy of Sciences 103, 6230 (2006).
  * Vergnon  _et al._ [2012] R. Vergnon, E. H. Van Nes, and M. Scheffer, Emergent neutrality leads to multimodal species abundance distributions, Nature communications 3, 1 (2012).
  * Kessler and Shnerb [2015] D. A. Kessler and N. M. Shnerb, Generalized model of island biodiversity, Physical Review E 91, 042705 (2015).
  * Bunin [2016] G. Bunin, Interaction patterns and diversity in assembled ecological communities (2016), arXiv [Preprint] https://arxiv.org/abs/1607.04734.
  * Roy  _et al._ [2020] F. Roy, M. Barbier, G. Biroli, and G. Bunin, Complex interactions can create persistent fluctuations in high-diversity ecosystems, PLoS computational biology 16, e1007827 (2020).
  * Fisher and Mehta [2014] C. K. Fisher and P. Mehta, The transition between the niche and neutral regimes in ecology, Proceedings of the National Academy of Sciences 111, 13111 (2014).
  * Verberk [2011] W. Verberk, Explaining general patterns in species abundance and distributions, Nature Education Knowledge 3, 38 (2011).
  * Fowler and Ruokolainen [2013] M. S. Fowler and L. Ruokolainen, Colonization, covariance and colour: Environmental and ecological drivers of diversity–stability relationships, Journal of theoretical biology 324, 32 (2013).
  * Barabás  _et al._ [2016] G. Barabás, M. J. Michalska-Smith, and S. Allesina, The effect of intra-and interspecific competition on coexistence in multispecies communities, The American Naturalist 188, E1 (2016).
  * Evans  _et al._ [2013] S. N. Evans, P. L. Ralph, S. J. Schreiber, and A. Sen, Stochastic population growth in spatially heterogeneous environments, Journal of mathematical biology 66, 423 (2013).
  * Tejo  _et al._ [2021] M. Tejo, C. Quiñinao, R. Rebolledo, and P. A. Marquet, Coexistence, dispersal and spatial structure in metacommunities: a stochastic model approach, Theoretical Ecology 14, 279 (2021).
  * Alonso  _et al._ [2006] D. Alonso, R. S. Etienne, and A. J. McKane, The merits of neutral theory, Trends in ecology & evolution 21, 451 (2006).
  * Haegeman and Loreau [2011] B. Haegeman and M. Loreau, A mathematical synthesis of niche and neutral theories in community ecology, Journal of theoretical biology 269, 150 (2011).
  * Baxter  _et al._ [2007] G. J. Baxter, R. A. Blythe, and A. J. McKane, Exact solution of the multi-allelic diffusion model, Mathematical bioSciences 209, 124 (2007).
  * Xu and Chou [2018] S. Xu and T. Chou, Immigration-induced phase transition in a regulated multispecies birth-death process, Journal of Physics A: Mathematical and Theoretical 51, 425602 (2018).
  * Capitán  _et al._ [2020] J. A. Capitán, S. Cuenda, and D. Alonso, Competitive dominance in plant communities: Modeling approaches and theoretical predictions, Journal of Theoretical Biology 502, 110349 (2020).
  * Black and McKane [2012] A. J. Black and A. J. McKane, Stochastic formulation of ecological models and their applications, Trends in ecology & evolution 27, 337 (2012).
  * Chesson [2018] P. Chesson, Updates on mechanisms of maintenance of species diversity, Journal of ecology 106, 1773 (2018).
  * MacArthur [1970] R. MacArthur, Species packing and competitive equilibrium for many species, Theoretical population biology 1, 1 (1970).
  * Chesson [1990] P. Chesson, Macarthur’s consumer-resource model, Theoretical Population Biology 37, 26 (1990).
  * O’Dwyer [2018] J. P. O’Dwyer, Whence lotka-volterra?, Theoretical Ecology 11, 441 (2018).
  * Allesina and Tang [2012] S. Allesina and S. Tang, Stability criteria for complex ecosystems, Nature 483, 205 (2012).
  * Gardiner  _et al._ [1985] C. W. Gardiner _et al._ ,  _Handbook of stochastic methods_ , Vol. 3 (springer Berlin, 1985).
  * Grimmett and Stirzaker [2001] G. Grimmett and D. Stirzaker,  _Probability and Random Processes_ (Oxford University Press, 2001).
  * Schnakenberg [1976] J. Schnakenberg, Network theory of microscopic and macroscopic behavior of master equation systems, Reviews of Modern physics 48, 571 (1976).
  * Meyn and Tweedie [1993] S. P. Meyn and R. L. Tweedie, Stability of markovian processes iii: Foster–lyapunov criteria for continuous-time processes, Advances in Applied Probability 25, 518 (1993).
  * Gupta  _et al._ [2014] A. Gupta, C. Briat, and M. Khammash, A scalable computational framework for establishing long-term behavior of stochastic reaction networks, PLoS computational biology 10, e1003669 (2014).
  * Marquet  _et al._ [2020] P. A. Marquet, M. Tejo, and R. Rebolledo, What is the species richness distribution?, in  _Unsolved problems in ecology_ , edited by A. Dobson, R. D. Holt, and D. Tilman (Princeton University Press, USA, 2020).
  * MacArthur and Wilson [1967] R. H. MacArthur and E. O. Wilson,  _The Theory of Island Biogeography: By Robert H. MacArthur and Edward O. Wilson_ (Princeton University press, 1967).
  * Goyal  _et al._ [2015] S. Goyal, S. Kim, I. S. Chen, and T. Chou, Mechanisms of blood homeostasis: lineage tracking and a neutral model of cell populations in rhesus macaques, BMC biology 13, 1 (2015).
  * Hu  _et al._ [2021] J. Hu, D. R. Amor, M. Barbier, G. Bunin, and J. Gore, Emergent phases of ecological diversity and dynamics mapped in microcosms, https://www.biorxiv.org/content/early/2021/10/29/2021.10.28.466339.full.pdf (2021), biorXiv [Preprint] https://doi.org/10.1101/2021.10.28.466339.
  * Iyer-Biswas and Zilman [2016] S. Iyer-Biswas and A. Zilman, First-passage processes in cellular biology, Advances in chemical physics 160, 261 (2016).
  * Redner [2001] S. Redner,  _A guide to first-passage processes_ (Cambridge university press, 2001).
  * May [1972] R. M. May, Will a large complex system be stable?, Nature 238, 413 (1972).
  * Allesina and Pascual [2008] S. Allesina and M. Pascual, Network structure, predator–prey modules, and stability in large food webs, Theoretical Ecology 1, 55 (2008).
  * Bunin [2017] G. Bunin, Ecological communities with lotka-volterra dynamics, Physical Review E 95, 042414 (2017).
  * Weigelt  _et al._ [2010] A. Weigelt, E. Marquard, V. M. Temperton, C. Roscher, C. Scherber, P. N. Mwangi, S. Von Felten, N. Buchmann, B. Schmid, E.-D. Schulze, _et al._ , The jena experiment: six years of data from a grassland biodiversity experiment, Ecology 91, 930 (2010).
  * Durán  _et al._ [2021] D. C. Durán, C. A. Hernández, E. Suesca, R. Acevedo, I. M. Acosta, D. A. Forero, F. E. Rozo, and J. M. Pedraza, Slipstreaming mother machine: A microfluidic device for single-cell dynamic imaging of yeast, Micromachines 12, 4 (2021).
  * Bacaër [2011] N. Bacaër, Yule and evolution (1924), in  _A Short History of Mathematical Population Dynamics_ (Springer, 2011) pp. 81–88.
  * Jeraldo  _et al._ [2012] P. Jeraldo, M. Sipos, N. Chia, J. M. Brulc, A. S. Dhillon, M. E. Konkel, C. L. Larson, K. E. Nelson, A. Qu, L. B. Schook, _et al._ , Quantification of the relative roles of niche and neutral processes in structuring gastrointestinal microbiomes, Proceedings of the National Academy of Sciences 109, 9692 (2012).
  * Dornelas and Connolly [2008] M. Dornelas and S. R. Connolly, Multiple modes in a coral species abundance distribution, Ecology Letters 11, 1008 (2008).
  * Callaghan  _et al._ [2021] C. T. Callaghan, S. Nakagawa, and W. K. Cornwell, Global abundance estimates for 9,700 bird species, Proceedings of the National Academy of Sciences 118 (2021).
  * Ser-Giacomi  _et al._ [2018] E. Ser-Giacomi, L. Zinger, S. Malviya, C. De Vargas, E. Karsenti, C. Bowler, and S. De Monte, Ubiquitous abundance distribution of non-dominant plankton across the global ocean, Nature ecology & evolution 2, 1243 (2018).
  * Matthews  _et al._ [2014] T. J. Matthews, P. A. Borges, and R. J. Whittaker, Multimodal species abundance distributions: a deconstruction approach reveals the processes behind the pattern, Oikos 123, 533 (2014).
  * Oakes  _et al._ [2017] T. Oakes, J. M. Heather, K. Best, R. Byng-Maddick, C. Husovsky, M. Ismail, K. Joshi, G. Maxwell, M. Noursadeghi, N. Riddell, _et al._ , Quantitative characterization of the t cell receptor repertoire of naïve and memory subsets using an integrated experimental and computational pipeline which is robust, economical, and versatile, Frontiers in immunology 8, 1267 (2017).
  * Descheemaeker and de Buyl [2020] L. Descheemaeker and S. de Buyl, Stochastic logistic models reproduce experimental time series of microbial communities, Elife 9, e55650 (2020).
  * Li and Chesson [2018] Y. M. Li and P. Chesson, Seed demographic comparisons reveal spatial and temporal niche differentiation between native and invasive species in a community of desert winter annual plants, Evolutionary Ecology Research 19, 71 (2018).
  * Chesson and Kuang [2008] P. Chesson and J. J. Kuang, The interaction between predation and competition, Nature 456, 235 (2008).
  * HilleRisLambers  _et al._ [2012] J. HilleRisLambers, P. B. Adler, W. S. Harpole, J. M. Levine, M. M. Mayfield, _et al._ , Rethinking community assembly through the lens of coexistence theory, Annual Review of Ecology, Evolution and Systematics 43, 2012 (2012).


Fluctuations and Selection Bias in 5 and 13 TeV
p-p Collisions: Where are the jets?

Thomas A. Trainor1*

1 University of Washington, Seattle, USA

* CorrespondingAuthor<EMAIL_ADDRESS>
## Abstract

The ALICE collaboration recently reported high-statistics $\bf p_{t}$ spectra
from 5 TeV and 13 TeV p-p collisions with intent to determine the role of jets
in high-multiplicity collisions. In the present study a two-component (soft +
hard) model (TCM) of hadron production in p-p collisions is applied to ALICE
$\bf p_{t}$ spectra. As in previous TCM studies of A-B collision systems jet
and nonjet contributions to $\bf p_{t}$ spectra are accurately separated over
the entire $\bf p_{t}$ acceptance. The statistical significance of data-model
differences is established leading to insights concerning selection bias and
spectrum model validity.

##  1 Introduction

A trend has emerged to interpret certain phenomena commonly associated with
“collectivity” in A-A collisions as evidence for QGP formation in smaller
collision systems if the same phenomena appear there. To address claims of
collectivity in $p$-$p$ collisions the present study applies a two-component
(soft + hard) model (TCM) of hadron production near midrapidity [1, 2] to
evaluate jet production in 13 TeV p-p $p_{t}$ spectra [3]. The study reveals
how jets contribute to spectra, especially at lower $p_{t}$, and what
selection biases result from event sorting according to two different $\eta$
acceptances. Z-scores are employed to determine the statistical significance
of data-model differences. Based on Z-scores four models are evaluated: fixed
TCM, variable TCM, Tsallis model and blast-wave model. Spectrum data appear to
be inconsistent with collectivity (flows) or jet modification.

##  2 Published 13 TeV p-p $\bf p_{t}$ spectra vs TCM reference

The ALICE collaboration recently published high-statistics $p_{t}$ spectra for
5 and 13 TeV $p$-$p$ collisions [3]. The stated purpose was “…to investigate
the importance of jets in high-multiplicity pp collisions and their
contribution to charged-particle production at low $p_{T}$.”

Figure 1: Published $p_{t}$ spectra (points) and TCM (solid) for 13 TeV
$p$-$p$ collisions and for V0M (left) and SPD (right) event selection together
with data/TCM data-model ratios.

Figure 1 shows published 13 TeV spectra (points) and fixed two-component model
(TCM) for two methods (V0M and SPD $\eta$ acceptances) of sorting the same
event (and jet) ensemble. The TCM provides excellent descriptions of data at
lower $p_{t}$, but systematic deviations depending on the selection method
(V0M or SPD $\eta$ acceptance) arise at higher $p_{t}$ and for lower $n_{ch}$.
The TCM is then used as a fixed reference to study selection bias.

##  3 Two-component Model – Jets vs Nonjet Production

The two-component (soft + hard) model (TCM) of hadron production near
midrapidity is a predictive model that can be applied to any A-B collision
system [4, 5, 7, 6]. Its simple algebraic structure is described by (with
$\bar{\rho}_{x}\equiv n_{x}/\Delta\eta$)

$\displaystyle\bar{\rho}_{0}(y_{t};n_{ch})$ $\displaystyle\approx$
$\displaystyle\frac{d^{2}n_{ch}}{y_{t}dy_{t}d\eta}~{}\approx~{}\bar{\rho}_{s}\hat{S}_{0}(y_{t})+\bar{\rho}_{h}\hat{H}_{0}(y_{t}),$
(1)

where transverse rapidity $y_{ti}\equiv\ln[(m_{ti}+p_{t})/m_{i}]$ is defined
for hadron species $i$ and $\hat{S}_{0}(y_{t})$ and $\hat{H}_{0}(y_{t})$ are
unit-normal fixed model functions. The model is based on (a) factorization of
$n_{ch}$ and $y_{t}$ dependence for each component and (b) the critical
relation $\bar{\rho}_{h}\approx\alpha(\sqrt{s})\bar{\rho}_{s}^{2}$ with
$\alpha(\sqrt{s})\approx O(0.01)$, both features originally inferred from
spectrum data [4]. The soft component is interpreted to represent projectile-
nucleon dissociation along the $p$-$p$ collision axis. The hard component is
interpreted to represent large-angle scattered low-$x$ gluons fragmenting to
dijets. The direct relation between spectrum hard components and measured jet
properties has previously been demonstrated quantitatively [8, 9].

Figure 2: Normalized $p_{t}$ spectra $X(y_{t})$ and inferred jet-related hard
components $Y(y_{t})$ for V0M (left) and SPD (right) event selection.
$\hat{S}_{0}(y_{t})$ and $\hat{H}_{0}(y_{t})$ are fixed model functions.

Figure 2 illustrates decomposition of measured $p_{t}$ spectra into soft and
hard components. Normalized spectra $X(y_{t})$ as defined by the first
relation in Eqs. (2) are shown in first and third panels. The normalized
spectra all coincide with TCM model $\hat{S}_{0}(y_{t})$ (bold dashed) within
data uncertainties below 0.5 GeV/c ($y_{t}\approx 2$). The hard/soft density
ratio is defined by
$x(n_{s})\equiv\bar{\rho}_{h}/\bar{\rho}_{s}\approx\alpha\bar{\rho}_{s}$.

$\displaystyle
X(y_{t})\equiv\frac{\bar{\rho}_{0}(y_{t})}{\bar{\rho}_{s}}\approx\hat{S}_{0}(y_{t})+x(n_{s})\hat{H}_{0}(y_{t});~{}~{}~{}Y(y_{t})\equiv\frac{1}{x(n_{s})}\left[X(y_{t})-\hat{S}_{0}(y_{t})\right]\approx\hat{H}_{0}(y_{t}).$
(2)

Spectrum hard components $Y(y_{t})$ as defined by the second relation in Eqs.
(2) are shown in second and fourth panels compared to model function
$\hat{H}_{0}(y_{t})$ (bold dashed). The inferred spectrum hard components
accurately represent the full jet contribution to spectra over the entire
$p_{t}$ acceptance subject to certain biases depending on event selection
criteria [1, 2].

##  4 Measures of Spectrum Structure and Bias Trends

Reference [3] introduces several methods to assess jet contributions to
spectra. One is plots of data spectra in ratio to a common reference data
spectrum. Another is fitting a power-law model to data spectra to infer
exponent $n$ possibly related to a jet contribution.

Figure 3 (left) shows spectrum ratios $X_{i}(y_{t})/X_{ref}(y_{t})$ [see Eq.
(2)] where $X_{ref}(y_{t})$ corresponds in this example to event class 5. Even
with the normalized form $X(y_{t})$ such spectrum ratios confuse two remaining
issues: jet production and selection bias. Contrast with Fig. 1 (second and
fourth) where selection bias can be studied in isolation.

Figure 3:  Left: Normalized $p_{t}$ spectra $X_{i}(y_{t})$ in ratio to
reference spectrum $X_{5}(y_{t})$. Right: Hard-component logarithmic
derivatives and power-law exponents $n$ for V0M and SPD.

A power-law exponent $n$ can be inferred directly from spectrum data in a
model-independent way (no power-law model fits required) via the logarithmic
derivative

$\displaystyle-\frac{d\ln[Y(p_{t})]}{d\ln(p_{t})}\rightarrow n~{}~{}\text{for
large $p_{t}$},$ (3)

where $Y(p_{t})=y_{t}Y(y_{t})/m_{t}p_{t}$. The result for data hard component
$Y(p_{t})$ must be the same as for intact data spectra because soft component
$\hat{S}_{0}(p_{t})$ is negligible at relevant $p_{t}$ values.

Figure 3 (third) shows log derivatives for V0M spectrum hard components. Above
$y_{t}=4$ ($p_{t}\approx 4$ GeV/c) inferred $n$ values stabilize near $n=5.6$
(for 13 TeV). Figure 3 (fourth) shows the major difference between V0M (fixed
$n$) and SPD (strongly varying $n$) trends: different responses to fluctuating
parton fragmentation from different $\eta$ intervals.

##  5 Statistical Significance and Z-scores

In comparisons of data with models it is essential to determine the
statistical significance of data-model deviations. The appropriate measure is
the Z-score [10]. The data/model ratio is a common method of data-model
comparison which however typically conveys a misleading representation of
model quality. The essential relations are conveyed by

$\displaystyle\frac{\text{data}}{\text{model}}-1\approx\frac{\text{data $-$
model}}{\text{error}}\times\frac{\text{error}}{\text{data}},$ (4)

where the approximation arises from replacing error/model by error/data. The
first factor at right defines Z-scores as data-model differences in ratio to
statistical uncertainties.

Figure 4:  First: Data/model ratios for V0M selection. Second: Corresponding
Z-scores. Third: Error/data ratios. Fourth: Z-scores for variable TCM in
relation to V0M data.

Figure 4 (first) shows V0M data/model ratios as in Fig. 1 (second). The second
panel shows the corresponding Z-scores while the third panel shows the
error/data ratios for V0M spectra corresponding to 60 million $p$-$p$ events.
The fourth panel shows Z-scores for a variable TCM in which the hard-component
model is adjusted to accommodate data (solely by varying the model width below
its mode). Surviving residuals are common to all event classes and both event
selection methods suggesting the presence of an artifact introduced during
data processing. Note the matching noise structure in Fig. 3 (third).

##  6 Alternative Spectrum Models

Certain models have been applied to $p_{t}$ spectra from small collision
systems in connection with claims of collectivity in those systems, for
instance Tsallis and blast-wave (BW) models. Those models were applied
recently to 13 TeV $p$-$p$ spectra as reported in Ref. [11].

###  6.1 Tsallis model – extent of equilibration

Figure 5 (first) shows Tsallis model fits (solid) compared with V0M spectra.
The second and third panels show corresponding data/model ratios and Z-scores.
The fourth panel shows Z-scores for the variable TCM on the same vertical
scale. The relation between the $\chi^{2}$ statistic and Z-scores is
$\chi^{2}=\sum_{i}Z_{i}^{2}$. Thus, $\chi^{2}$ values for the Tsallis model
are in this case $O(10,000)$ [2]. The Tsallis model is dramatically rejected
by spectrum data.

Figure 5:  First: Tsallis model (solid) vs V0M spectra (points). Second:
Tsallis data/model ratios. Third: Corresponding Z-scores. Fourth: Z-scores for
variable TCM.

###  6.2 Blast-wave model – hydrodynamic flows

Figure 6 (first) shows BW model fits (solid) to the same V0M spectra. The
second and third panels show corresponding data/model ratios and Z-scores.
Again one encounters $\chi^{2}$ values $O(10,000)$ for the BW model. The
fourth panel shows the same BW model applied to 5 TeV $p$-Pb data compared to
a TCM for identified hadrons reported in Ref. [6]. It is notable that whereas
the BW model (dashed) misses the $p$-Pb data to the same degree as in the
first panel the TCM for $K^{0}_{\text{S}}$ (solid) describes the data within
their uncertainties from 7 GeV/c down to zero, precluding any possibility of a
radial-flow contribution [2].

Figure 6:  First: Blast-wave model (solid) vs V0M spectra (points). Second:
Data/model ratios. Third: Corresponding Z-scores. Fourth: BW model applied to
5 TeV $p$-Pb spectra.

##  7 Conclusion

The same ensemble of 60 million 13 TeV $p$-$p$ collisions is partitioned in
two ways (V0M vs SPD). A two-component model (TCM) facilitates precise
separation of $p_{t}$ spectra into nonjet (soft) and jet-related (hard)
components over the entire $p_{t}$ acceptance. Complementary biases of the
jet-related hard component relative to the fixed TCM as reference are observed
reflecting the relation of different $\eta$ acceptances to jet production
mechanisms. A variable TCM (variation of hard-component model widths below or
above the mode) is adjusted to accommodate data. Z-scores (data-model
deviations in ratio to statistical uncertainties) are introduced to evaluate
deviation significance. Z-scores are then used to test model validity. The
variable TCM is observed to describe spectra within their point-to-point
uncertainties whereas Tsallis and blast-wave models, often invoked in the
context of “collectivity” and incomplete thermal equilibration, are
dramatically falsified with Z-scores equivalent to $\chi^{2}$ values
$O(10,000)$. Careful examination of spectrum evolution with charge-
multiplicity and pseudorapidity conditions reveals no significant evidence for
radial flow or for jet modification that might buttress claims of $p$-$p$
collectivity.

## References

  * [1] T. A. Trainor, arXiv:2104.08423 [hep-ph].
  * [2] T. A. Trainor, arXiv:2107.10899 [hep-ph].
  * [3] S. Acharya et al. (ALICE collaboration), Eur. Phys. J. C 79, no.10, 857 (2019).
  * [4] J. Adams et al. (STAR Collaboration), Phys. Rev. D 74, 032006 (2006).
  * [5] T. A. Trainor, Int. J. Mod. Phys. E 17, 1499 (2008).
  * [6] T. A. Trainor, J. Phys. G 47, no.4, 045104 (2020).
  * [7] T. A. Trainor, J. Phys. G 44, 075008 (2017).
  * [8] T. A. Trainor, Phys. Rev. C 80, 044901 (2009).
  * [9] T. A. Trainor, Phys. Rev. D 89, no.9, 094011 (2014).
  * [10] E. Kreyszig (1979). Advanced Engineering Mathematics (Fourth ed.), Wiley, p. 880, Eq. 5.
  * [11] R. Rath, A. Khuntia, R. Sahoo and J. Cleymans, J. Phys. G 47, no.5, 055111 (2020).