# Self-segregation in heterogeneous metapopulation landscapes

Jean-François de Kemmeter1, Timoteo Carletti1, Malbor Asllani2,3 $1$ naXys,
Namur Institute for Complex Systems, & Department of Mathematics, University
of Namur, rue Grafé, 2 B5000, Belgium
$2$ School of Mathematics and Statistics, University College Dublin, Belfield,
Dublin 4, Ireland
$3$ Department of Mathematics, Florida State University, 1017 Academic Way,
Tallahassee, FL 32306, United States of America

###### Abstract

Complex interactions are at the root of the population dynamics of many
natural systems, particularly for being responsible for the allocation of
species and individuals across apposite niches of the ecological landscapes.
On the other side, the randomness that unavoidably characterises complex
systems has increasingly challenged the niche paradigm providing alternative
neutral theoretical models. We introduce a network-inspired metapopulation
individual-based model (IBM), hereby named _self -segregation_, where the
density of individuals in the hosting patches (local habitats) drives the
individuals spatial assembling while still constrained by nodes’ saturation.
In particular, we prove that the core-periphery structure of the networked
landscape triggers the spontaneous emergence of vacant habitat patches, which
segregate the population in multistable patterns of isolated (sub)communities
separated by empty patches. Furthermore, a quantisation effect in the number
of vacant patches is observed once the total system mass varies continuously,
emphasising thus a striking feature of the robustness of population stationary
distributions. Notably, our model reproduces the patch vacancy found in the
fragmented habitat of the Glanville fritillary butterfly _Melitaea cinxia_ ,
an endemic species of the Åland islands. We argue that such spontaneous
breaking of the natural habitat supports the concept of the highly contentious
(Grinnellian) niche vacancy and also suggests a new mechanism for the
endogeneous habitat fragmentation and consequently the peripatric speciation.

##  I Introduction

The distinct feature of complex systems is their critical reliance on the
interactions between the constituent entities. Such strong dependence has
given birth to the discipline of network science [1, 2, 3, 4]. From this
perspective, nodes represent either single individuals or groups thereof and
edges, the respective interactions among them. According to the positive,
negative or neutral nature of such interactions, they can identify mutualism,
competition, or commensalism, respectively [5, 6]. There is a long tradition
in population dynamics, further strengthened by complexity science, to
attribute collective behaviors such as self-organisation or adaption, in
biological or social systems, to interactions between entities (e.g., Lotka-
Volterra models) [5, 1, 6]. In ecology, the coexistence of a large variety of
species is thought to be possible because of their ability to occupy specific
niches, namely the set of environmental factors (resources, prey/predators,
competitors, abiotic components, etc.) matched with a given species, as a
consequence of a long and highly selective evolution process [7, 8].

An alternative approach adopted in recent years to tackle the question of
species diversity is that of neutral theory, according to which organised
patterns can occur exclusively due to the randomness of several factors that
decisively condition species survivability [9, 7, 8, 10, 11, 12]. Stochastic
processes are the common ground around which the modern neutral theory of
biodiversity is erected [9, 8, 10]. In this regard, intra- and inter-specific
interactions are considered simultaneously relevant, consisting in a striking
difference of the neutral model from niche ones [13]. Thus elementary
processes such as births/deaths or dispersal (emigration and immigration)
result crucial in describing the collective dynamics emerging from individual-
based models (IBM). In the context of neutral theory, there have been several
attempts to better understand the segregation of species in ecological niches
in terms of the randomness of the decisive factors in action during the
evolutionary process [9, 6]. For instance, there is a consolidated debate
whether the species fully occupy all the niches to overcome Gause’s law of
competitive exclusion [5, 8, 10] or either vacant niches can exist due to an
optimal fitness landscape [6]. According to the latter thesis, the fitness
function that describes the adaption of the species to the niche domain, is
rather rugged and multipeaked. This way, the randomness associated with the
initial conditions will allow the individuals assembling around some (local)
optima and eventually allow the emergence of vacant patches around the
landscape minima. On the other side, ecological landscapes found in nature are
highly fragmented and are often conditioned by different geographic, abiotic,
or human factors [8, 7, 14, 15, 16, 17, 18, 19, 20]. Understanding the
mechanisms of fragmentation of the ecological landscapes can give valuable
insight into the the evolution of animals, plants, and other living beings in
their respective habitats [21, 22, 23]. By anticipating on the following let
us emphasize that the proposed model is capable to reproduce the behavior we
can observe in some empirical ecosystems, e.g., the Glanville fritillary
butterfly _Melitaea cinxia_ in the Åland islands and the giant kelp
_Macrocystis pyrifera_ in southern California [16, 17, 18, 19, 20, 15, 24]. In
both cases, the landscape is naturally fragmented in several available habitat
patches, however empirical observations show that most of them remain empty
because animals spontaneously select a subset of them where to live.

In the fashion of neutral theory we here consider a single species model (or
of indistinguishable individuals) capable to reproduce the emergence of empty
niches hereby strictly intended as vacant (local) habitats following
Grinnellian formalism [25], and the resulting fragmentation of the available
habitat into smaller isolated ones [26, 27]. To this aim, we focus exclusively
on a dispersal process, deliberately neglecting other types of dynamics (e.g.,
birth/death processes). Doing so, we assume that any aspatial intra- or
interspecific dynamics occur considerably faster than the spatial one.
Furthermore, such local dynamics are assumed not to be affected by exogenous
factors, thus overall conserving the total mass of the species involved (see
Sec. V for further details). The model we propose can help to better
understand the potential mechanisms that lay the ground for phenomena such as
peripatric speciation or endogenous habitat fragmentation. It is important to
emphasize that in the latter case, the fragmentation occurs exclusively due to
the population dynamics of the species rather than exogenous factors that
influence the habitat landscape. Our approach is reminiscent of the celebrated
Schelling’s model [28, 29, 30], according to which self-organised patterns of
segregation spontaneously emerge due to the interactions between (at least)
two antagonistic types of individuals. At variance with such paradigm, and in
the spirit of neutral theory, we hereby assume all individuals to be
equivalent (neutral paradigm) [13]. Within this framework, we show that the
fragmentation of the habitat is now exclusively due to the spatial
heterogeneity of the networked landscape, particularly the presence of
multiple cores linked together through bridge (peripheral) nodes. The paper is
organised as follows: in Sec. II, we introduce the model and derive the mean-
field equations describing the dispersal of individuals across the patches. We
investigate the long-term distribution of the population sizes and the
emergence of vacant patches in Sec. III. In Sec. IV, we illustrate our results
using real and synthetic networks of habitat patches as a substrate of the
dynamical process. In particular, we show how our model can result useful in
understanding the habitat fragmentation of _Glanville fritillary_ butterfly, a
species endemic of Åland islands. We then conclude and further discuss the
perspective impact of our work.

##  II Self-segregation process

We consider a metapopulation network, whose nodes are supposed to represent
the existing spatial patches or the local habitats connected through edges
mimicking physical paths bridging distant patches. In our model, the only
random process is the dispersal of individuals, assumed to be
indistinguishable, driven by the affinity that individuals of different but
adjacent nodes exhibit for each other. We thus assume an overall carrying
capacity that sets the total population size at an (average) fixed value. As
anticipated earlier, this assumption holds true when the timescale of the
local dynamics on each of the spatial patches is disentangled from the global
one of the dispersal across the different patches [31]. For instance, if we
assume the local birth/death processes to be very slow, the total population
can be safely assumed to be constant during the whole diffusion process. As we
will illustrate in the following through an empirical example, another
plausible alternative is to consider that such local dynamics are much faster
than the dispersal one. Indeed, the last case can allow splitting the dynamics
of the model into two steps, the first (very short) phase where the local
populations reach a stationary value in each node, followed by a second where
individuals move among the habitats. The interactions among individuals of
adjacent patches are grounded on the mutual benefit that individuals have when
sharing the same biogeographical area or territory. Examples span from more
extreme cases such as the clumping of individuals of the same species to
prevent the Allee effect [32, 33] to less drastic ones such as the sharing of
habitat of territorial animals during the mating season [34]. We will denote
by _self- segregation_ the positive influence that the indistinguishable
individuals exert on each other. On the other hand, to model/introduce the
habitat saturation conditioned by limited (abiotic) resources or negative
interactions such as intra- or inter-species competition, etc., we assume that
the nodes have a maximum carrying capacity of possible individuals to
allocate. Hence, populated spatial patches are highly likely to attract new
agents as long as they do not become too crowded. The resulting collective
behavior is the emergence of different configurations of isolated
(sub)communities  \- each (sub)community consisting of a set of connected and
occupied spatial patches - surrounded by vacant patches of the habitat. The
terminology (sub)community, as used here, is borrowed from network science
[1], where it emphasizes the presence of strongly connected (sub)networks
clearly distinguished from each other. Let us nevertheless mention that an
alternative definition in ecology would be a separated/isolated
(sub)population [13]. However, for reasons of wording compactness, throughout
this paper, we have preferred to use the term (sub)community instead.

We start by considering a simple (i.e., without self-loops or multi-edges)
connected and undirected graph made of $\Omega$ nodes $\\{v_{i}\\}_{\Omega}$
whose structure is encoded in its adjacency matrix entries, $A_{ij}=1$ if
there is a link between nodes $v_{i}$ and $v_{j}$, and zero otherwise. A
fundamental assumption of our model is that the network is characterised by a
heterogeneous degree distribution $\\{p(k_{i})\\}_{\Omega}$ where the degree
of node $v_{i}$ is defined as $k_{i}=\sum_{j=1}^{\Omega}A_{ij}$. The
population size at time $t$ within node $v_{i}$ is written $n_{i}(t)$. If we
denote the state of the system at time $t$ by
$\mathbf{n}(t)\triangleq(n_{1}(t),n_{2}(t),\cdots,n_{\Omega})$ and the
probability of observing such state by $P(\mathbf{n},t)$, then the
mathematical formalism that governs the individuals dynamics is described by
the master equation [35, 36]:

$\frac{\text{d}P(\mathbf{n},t)}{\text{d}t}=\sum_{\mathbf{n^{\prime}}\neq\mathbf{n}}\Big{[}T(\mathbf{n}|\mathbf{n^{\prime}})P(\mathbf{n^{\prime}},t)-T(\mathbf{n^{\prime}}|\mathbf{n})P(\mathbf{n},t)\Big{]},\,$
(1)

where $T(\mathbf{n^{\prime}}|\mathbf{n})$ stands for the transition
probability from state $\mathbf{n}$ to state $\mathbf{n^{\prime}}$. Following
the idea of the finite carrying capacity of the habitat patches (and similarly
to [37, 38]), we impose to each node $v_{i}$ a maximum number of agents $1\leq
n_{i}\leq N$ that can be hosted. Also to focus exclusively on the effect of
the attractiveness that agents exert on each other and at the same time their
intention to avoid overcrowding, we assume that individuals are not allowed
for other actions. With a slight abuse of notation, the transition from node
$v_{i}$ to node $v_{j}$ is written as:

$T(n_{i}-1,n_{j}+1|n_{i},n_{j})=\frac{A_{ij}}{k_{i}}\frac{n_{i}}{N}g\left(\frac{n_{j}}{N}\right),$
(2)

where the function $g(\cdot)$ represents the probability for the agents to
settle in the chosen host node and quantifies the reciprocal benefit among
individuals while still considering the finite size of the nodes. A
significant difference with other random processes, e.g., the biased random
walks [39], is that the probability of choosing the host node is independent
of the densities of individuals in the other neighbour nodes 11footnotetext:
In a biased random walk, the motion of individuals depends on the relative
attribute of the hosting node compared to all the neighbour ones [39, 51].
[40]. This subtle but crucial feature makes sense for the spatial
interpretation of the habitat support [41, 11]: the individual cannot a priori
choose the most suitable node before first “testing the ground”. Said
differently, the agents perceive the pressure of limited resources or other
individuals of the same kind once they reside on a given patch, but not
earlier. Throughout this paper, we will make a straightforward assumption
regarding the function $g(\cdot)$: the hosting node will recruit individuals
proportionally to the density of agents already present therein, and at the
same time, the probability for a recruited individual to settle in the
selected node, is proportional to the available free space (capacity
constraint). In formula we have

$g\left(\frac{n_{j}}{N}\right)=\frac{n_{j}}{N}\frac{N-n_{j}}{N},$ (3)

which constitutes the prominent logistic function and has plentiful
applications in many areas of science and in particular in ecology [5]. Since
looking for an exact solution of Eq. (1) results challenging in general, we
apply a standard procedure to obtain the deterministic description of the
self-segregation problem. We start by multiplying both sides of Eq. (1) by
$n_{i}$, and summing over $n_{i}$ for all $i$, we obtain the evolution
equation for the average density

$\langle
n_{i}(t)\rangle=\sum_{\mathbf{n}}n_{i}P(\mathbf{n},t)=\sum_{n_{i}}n_{i}\sum_{\mathbf{n}\setminus
n_{i}}P(\mathbf{n},t)=\sum_{n_{i}}n_{i}P(n_{i},t)$

recalling that $P(\mathbf{n},t)$ is the joint probability distribution
$P(n_{1},\dots,n_{\Omega},t)$. Regarding the r.h.s. of Eq. (1) instead, we
need to point out that the only available state reachable from
$\mathbf{n}(t)=(n_{1},n_{2},\cdots,n_{i},\dots,n_{j},\dots,n_{\Omega})$ in the
time interval $\Delta t$ is $\mathbf{n}(t)=(n_{1},n_{2},\cdots,n_{i}\mp
1,\dots,n_{j}\pm 1,\dots,n_{\Omega})$ where $v_{i}$ and $v_{j}$ are two
adjacent nodes. Thus by shifting the two indices of the two sums of Eq. (1)
respectively by $+1$ and $-1$ we obtain

$\frac{\text{d}\langle n_{i}\rangle}{\text{d}t}=\sum_{j=1}^{\Omega}\langle
T(n_{i}+1,n_{j}-1|n_{i},n_{j})\rangle-\langle
T(n_{i}-1,n_{j}+1|n_{i},n_{j})\rangle$

where for sake of simplicity we have ommitted the inactive terms $n_{l}$ and
where again $\langle
f(\mathbf{n},t)\rangle=\sum_{\mathbf{n}}f(\mathbf{n},t)P(\mathbf{n},t)$. To
obtain the deterministic equation, we proceed by dividing both sides by the
carrying capacity $N$ and rescale the time $t\mapsto t/N$. The thermodynamic
limit ${N\rightarrow+\infty}$ yields the mean-field (MF) equations:

$\frac{\text{d}\rho_{i}}{\text{d}t}=\sum_{j=1}^{\Omega}\mathcal{L}_{ij}\rho_{i}\rho_{j}\left[1-\rho_{i}-\dfrac{k_{j}}{k_{i}}(1-\rho_{j})\right],\,\forall
i,$ (4)

where $\rho_{i}\equiv\lim_{N\rightarrow+\infty}\langle n_{i}\rangle/N$ is the
population density within node $v_{i}$ and
$\mathcal{L}_{ij}=A_{ij}/k_{j}-\delta_{ij}$ corresponds to the random walk
(RW) Laplacian [1, 4]. Notice also that in the limit of large $N$ we drop any
correlation among different nodes i.e., $\langle n_{i}n_{j}\rangle\sim\langle
n_{i}\rangle\langle n_{j}\rangle$ based on the van Kampen ansatz
22footnotetext: Based on van Kampen ansatz, we can write the (rescaled)
discrete variable as $n_{i}/N=\rho_{i}+\xi_{i}/\sqrt{N}$, valid for large $N$,
where $\rho_{i}$ is the deterministic variable describing the density of node
$v_{i}$ and $\xi_{i}$ the stochastic variable. Thus $\langle
n_{i}n_{j}\rangle/N^{2}=\langle\rho_{i}\rho_{j}+\xi_{i}\xi_{j}/N+\left(\xi_{j}\rho_{i}+\xi_{i}\rho_{j}\right)/\sqrt{N}\rangle\sim\langle
n_{i}\rangle\langle n_{j}\rangle/N^{2}$ in the limit for large $N$ while the
densities are kept constant. [36, 42]. From Eq. (4), it can be easily shown
that the total population size is constant over time, i.e.,
$\sum_{i}\frac{\text{d}\rho_{i}}{\text{d}t}=0$, a consequence of the fact that
only dispersal processes are considered. We will define the average density
$\beta=\sum_{i}\rho_{i}/\Omega$ taking values in the interval $[0,1]$.

Figure 1: Self-segregation and the ensuing (sub)communities in a heterogeneous
network. In all the panels, nodes size and colour represent, respectively,
their degree and the asymptotic mass density (the larger the size the higher
the degree and the darker the colour the larger the density); moreover each
system has been initialised with a uniform density of individuals per node,
i.e. $\rho_{i}(0)=\beta$. When the total mass is sufficiently large (leftmost
panel), all the nodes are occupied and a single large community prevails
(emphasised with the shaded oval). However, once the average density $\beta$
starts decreasing (from the left panel to the right one), nodes are depleted
accordingly to their degree, the smallest the first, leaving only nodes with
higher degrees occupied. Among the nodes of degree $2$, the first one to
become empty is the one with the largest bridgeness [43] (see Appendix B). In
such a scenario, (connected) subgraphs with filled nodes give birth to
separated (sub)communities (emphasised with the shaded ovals) separated by of
empty nodes. For instance, when $\beta=0.15$, the mass segregates into two
subcommunities $\mathcal{C}_{1}$ ad $\mathcal{C}_{2}$ of respective sizes $4$
and $2$. In the subcommunity $\mathcal{C}_{2}$, the two nodes have distinct
stationary densities since they have different degrees. Notice also that while
the average density $\beta$ varies, the number of vacant nodes does not
necessarily change accordingly, a phenomenon we refer to as “quantisation” of
the occupancy (see Fig. 6 for further details).

##  III Fragmentation, multistability and resilience of (sub)communities

Starting from the mean-field Eq. (4) we can determine the equilibrium states
$\\{\rho_{i}^{*}\\}_{\Omega}$ and based on that we can afterwards perform a
(linear) stability analysis. Let us first note that any node might potentially
become vacant as $\rho_{i}^{*}=0$ is a fixed point of system (4)
$(i=1,\cdots,\Omega)$. Obviously, there is a limited number of nodes that can
become vacant as the total mass must remain constant while satisfying the
limiting carrying capacities of the nodes. As empty nodes emerge, the mass can
segregate into $M$ isolated subcommunities, the $m$-th ($m=1,\cdots,M$) one
being a subnetwork of $\Omega_{m}$ connected and filled nodes. Such
subcommunities are separated by empty patches that prevent individuals from
migrating to another subcommunity. A straightforward computation (the
interested reader can refer to the Appendix A) allows to obtain the population
density in the $m$-th one at equilibrium to be given by

$\rho_{i}^{*}=0\qquad\mathrm{or}\qquad\rho_{i}^{*}=1-\frac{C_{m}}{k_{i}},\,\forall
i,$ (5)

where $C_{m}=\left(1-\beta_{m}\right)/\langle 1/k\rangle_{\Omega_{m}}$ is a
conservation constant that depends on the initial configuration
$\\{\rho_{i}(0)\\}_{\Omega}$ and $\beta_{m}=\sum_{j\in
M_{m}}\rho_{j}^{*}/\Omega_{m}$ is the average stationary population density of
the $m$-th subcommunity $M_{m}$ with size $\Omega_{m}$. Notice that, in
general, the interpretation of the system behavior depends on the knowledge of
the conservation constant $C_{m}$ which cannot be a priori inferred. Eq. (5)
results thus an implicit equation for the stationary nodes’ densities. Let us
point out also that in equation (5), $k_{i}$ is the degree of node $v_{i}$
considering the whole network, not only the subcommunity to which it belongs.
In particular, in the case of a single community, $\beta_{1}=\beta$ is the
global density, i.e., the total number of agents divided by the number of non
empty nodes. Observe that the average $\langle 1/k\rangle_{\Omega_{m}}$, has
been performed over the $\Omega_{m}$ nodes forming the $m$-th subcommunity.
From relation (5), it is straightforward that nodes resulting occupied at
equilibrium should have a degree $k_{i}>C_{m}$ or will be vacant otherwise.
Because $C_{m}$ decreases with the total mass, the last observation suggests
that acting on the mass we can induce the agents to segregate at the nodes
with higher degrees, possibly leaving the ones with a lower degree vacant. To
elucidate the robustness of the (sub)communities emergence we will consider
the case when a single empty node, adjacent to a given (sub)community, is
slightly perturbed 33footnotetext: Obviously, an isolated, empty node
(surrounded by other empty ones) is always unstable since the mass added to it
will never leave the node.[44] and postpone a complete stability analysis to
Appendix A. We start by assuming the system settles on some equilibrium and we
focus on the $m$-th subcommunity, whose nodes are characterised by
$\rho_{j}^{*}=1-C_{m}/k_{j}>0$. Suppose it exists an empty node $v_{i}$
connected exclusively to the nodes of the $m$-th subcommunity to which we
allocate a small amount of mass, such that this new configuration still
satisfies the same constraints ($m$ and $\beta_{m}$ remain unchanged assuming
$\Omega_{m}$ is large enough to neglect the small amount of mass we have to
remove from each node but the $i$–th one to satisfy the constraint), that is
$\rho_{i}^{*}\mapsto\delta$ and $\rho_{j}^{*}=1-C_{m}/k_{j}$ for $j\neq i$. A
direct computation allows to obtain the linearised dynamics that governs node
$v_{i}$ at early times

$\frac{\text{d}\delta}{\text{d}t}\approx\delta\left(1-\frac{C_{m}}{k_{i}}\right)\sum_{\begin{subarray}{c}j=1\\\
\left(j\neq
i\right)\end{subarray}}^{\Omega_{m}}\mathcal{L}_{ij}\rho_{j}^{*}\,.$ (6)

From here we can conclude that the $v_{i}$-th node will increase its mass
(acquiring it from other $\Omega_{m}$ nodes) if $C_{m}<k_{i}$ and will
transfer it to the other nodes of the (sub)community $m$ otherwise. Similarly,
if we add a small mass $\delta$ to an occupied node, the system is always
stable. In fact, if we slightly perturb only the $i$-th node of the $m$-th
subcommunity, $\rho_{i}^{*}\mapsto\rho_{i}^{*}+\delta$, with $\delta>0$, a
straightforward computation allows writing
$\frac{\text{d}\delta}{\text{d}t}\approx-\delta\rho_{i}^{*}\sum_{j=1}^{\Omega}\mathcal{L}_{ij}\rho_{j}^{*}\,,$
thus proving the stability of the subcommunity. This is reasonable since the
added mass to some pre-existing (sub)community cannot escape from it being
surrounded by empty nodes acting as movement barriers.

The local analysis emphasises the role of the degree distribution
heterogeneity in the existence of stable states with empty nodes. For
instance, in a regular graph, the only possible (non-trivial) stable state is
the one being uniformly occupied since $C_{m}<k$ for all the nodes
44footnotetext: This can be easily obtained from the formula
$C_{m}=\left(1-\beta_{m}\right)k$ since the average
$\langle\cdot\rangle_{\Omega}$ over the degrees drops being the graph regular.
From here, it is immediate that $C_{m}/k=\left(1-\beta_{m}\right)<1$, which
justifies the stability of the fully occupied nodes state.[45]. To illustrate
the outcomes of the self-segregation dynamics, we refer to Fig. 1, where the
system has been initialised with a uniform distribution of nodes’ densities
for all the panels. From the previous analysis, if the average density is
sufficiently large, i.e., $\beta>1-{1}/{\Omega}$, all nodes are occupied at
equilibrium, resulting in a single community. However, as $\beta$ is gradually
reduced, some empty nodes (the least connected ones and those with the largest
bridgeness centrality [43]) start to emerge. Decreasing the total mass would
result in increasing the constant $C_{m}$, which in turn forces all the nodes
for which $k_{j}>C_{m}$ is no longer satisfied, to deplete. The mass
accumulates on fewer nodes in such a way to satisfy the new equilibrium
condition (5). Also from Fig. 1, one can notice that for (sufficiently) small
variations of the mass, the set of empty/occupied nodes remains the same,
showing that the system is resilient in switching to a state with more (less)
vacant nodes while the mass decreases (increases). Such “quantisation”
phenomenon is exclusively due to the discrete distribution of the network
support and can be understood starting from the conservation constant
$C_{m}=\left(1-\beta_{m}\right)/\langle 1/k\rangle_{\Omega_{m}}$ which changes
continuously with $\beta_{m}$, if the number of (sub)communities and empty
patches does not vary, that in turn depends on the total average density
$\beta$. Since the degrees are discrete and the network is finite, density
intervals, where $C_{m}$ does not overtake the next lowest degree of occupied
nodes, will certainly exist. Furthermore, the length of such intervals must
decrease when the differences between the successive degrees become smaller,
suggesting that a broader degree distribution will show a less pronounced
quantisation effect. Nevertheless, in the following, we will see that, for
non-regular networks, a node can be either empty or occupied depending both on
the initial conditions and the network structural features. The latter yield
to stable configurations with different numbers and distributions of vacant
nodes even for a fixed value of the metapopulation size, a phenomenon we
hereby refer to as _multistability_. The latter yield stable configurations
with different numbers and distributions of vacant nodes for a fixed value of
the metapopulation size, a phenomenon we hereby refer to as multistability.
More precisely, we will show the existence of regions of the phase space,
determined by $\beta$, i.e., the hyperplane
$\beta_{\Omega}=\sum_{i}\rho_{i}(0)$, containing several attractors, each one
with a relatively small stability basin thus even slightly changed initial
conditions can determine very different system outcomes, in particular with a
substantial change in the number of empty nodes.

##  IV Role of topology in the formation of (sub)communities

So far, we have shown that the heterogeneous degree distribution
characterizing the network supports is crucial in the fragmentation of the
habitat in isolated (sub)communities separated by empty patches. In this
section, we explore in detail the structural features responsible for the
breaking of the contiguous habitat and the emergence of vacant ones in both
synthetic and empirical ecological networks. In particular, we show that low
degree bridge nodes are responsible for the breaking of the original community
in smaller ones by being the first to deplete. Let us remember the definition
of bridge nodes, the latter being nodes with high betweenness centrality and
connecting different and distant regions in the network. They thus act as
global bridges allowing separate parts of the network to communicate, in
opposition with local bridges that connect nearby groups of nodes. Removal of
global bridges will thus easily disconnect the network. In the next two
sections we will consider the behavior of the above defined model on a
synthetic spatial network, the random geometric network, and on the empirical
network arising from the ecosystem of the Glanville fritillary butterfly in
the Åland islands. We further investigate the role of network topology on the
metapopulation breaking for other synthetic and real networks in Appendix B.

###  IV.1 Random geometric networks

We start by considering a family of random geometric graphs as a test bed for
numerically verifying the theory developed throughout this paper. Such spatial
networks have resulted useful for modelling the dispersal of individuals in
ecological settings [46, 47, 48]. Starting from a set of uniformly distributed
nodes in the Cartesian space, we connect every two nodes whose Euclidian
distance is lower than some threshold, the rationale being that individuals
are more likely to migrate to spatial patches that are not too far away.
Panels $a)$ and $b)$ of Fig. 2 correspond to random geometric (connected)
graphs obtained by drawing $\Omega=200$ points in the unit square and
connecting them if their euclidean distance is lower than $r=0.175$ or
$r=0.275$ for panel $a)$ and $b)$, respectively. This procedure results in two
distinct graphs on top of which we ran the dynamical system given by (4),
considering an average node population density $\beta=0.05$. In both cases, we
emphasize the distinct emerged (sub)communities by means of shaded areas.
While a single giant community is observed when $r=0.275$, multiple isolated
(sub)communities are found when $r=0.175$. In panel $c)$, we show the fraction
of empty nodes as a function of $\beta$, for random geometric graphs with
distinct values of the parameter $r$. We observe a quantisation effect in the
fraction of empty nodes. As $r$ increases, this quantisation effect becomes
more diluted since the degree distribution becomes sharper. The results were
averaged over different configuration of the initial distributions of the
densities $\rho_{i}(0)$, and the min-max deviation is reported by the shaded
area. Hence, the system is multistable, namely for the same value of $\beta$
and a fixed network, the asymptotic distribution of mass can vary because it
depends on the initial conditions.

|  |
---|---|---
(a) | (b) | (c)
Figure 2:  (a) Emergence of isolated (sub)communities separated by empty nodes
in a random geometric graph made of $\Omega=200$ nodes randomly drawn in the
unit square. Every pair of nodes is connected if their Euclidean distance is
less than $r=0.175$. White nodes correspond to empty nodes while colored nodes
correspond to filled ones. Each shaded area corresponds to a subset of
connected filled nodes. (b) The same but for $r=0.275$. (c) Fraction of empty
nodes in a random geometric graph made of $\Omega=200$ nodes, as a function of
$\beta$. Each curve corresponds to a distinct value of $r$.

###  IV.2 The Glanville fritillary metapopulation

In this section, we will switch our attention to an empirical spatial
ecological network representing the natural habitat of the Melitaea cinxia
species. Commonly known as Glanville fritillary butterfly, it is an endemic
species of the Åland islands in southwestern Finland [18]. Its habitat
consists of a fragmented landscape of meadows, pastures and rocky coastal
areas separated one from the other by uninhabitable areas [20]. Butterflies
disperse from patch to patch with an average migration distance of a few
hundred meters to a few kilometers [20, 16, 19]. Interestingly, an important
fraction ($75-80\%$) of these suitable habitat patches are left vacant, i.e.,
are unoccupied by the butterflies [20]. Although birth and death processes as
well as extinction and colonisations events occur at the level of the patches,
the size of the metapopulation has remained relatively stable over longer
timescales, supporting the idea of a fixed global carrying capacity [17, 18].
So far, around $4500$ habitat patches suitable for hosting this butterfly
species have been catalogued [20]. For the sake of clarity, we here restrict
our attention to a subset of them found within the black rectangle shown in
Fig. 3 $(a)$. We draw connections between nodes assuming an exponential
dispersal kernel; for every pair of distinct nodes, we connect them if their
geodesic distance (in kilometers) is less than a random number drawn from the
exponential distribution $f(x)=\lambda e^{-\lambda x}$, with $\lambda=0.65$.
This corresponds to an average dispersal distance of approximately
$\dfrac{1}{\lambda}\approx 1.5$km, in adequation with biological data [20, 16,
19]. As expected, when $\beta$ is low enough, the mass segregates into
numerous subnetworks separated by vacant patches, see panel $(b)$ obtained for
$\beta=0.08$. In this case, there are $222$ empty nodes and most of the mass
accumulates into three (sub)communities of respective sizes $57,56$ and $24$
nodes that correspond to the main cores of the network, as shown in the inset
of panel $(d)$. Panel $(c)$ shows similar results for $\beta=0.14$, where this
time a slightly greater total mass results in fewer vacant nodes ($170$) and
there is a single large community and few smaller ones. Panel $(d)$ shows the
fraction of empty nodes as a function of $\beta$ with the min-max deviation.

|
---|---
(a) | (b)
|
(c) | (d)
Figure 3:  $(a)$ Map of the Åland islands (southwestern Finland). The boundary
of Åland islands was retrieved from [49]. The suitable habitat patches [50] of
the Glanville fritillary butterfly are indicated by the red points. We here
restrict our attention to the subset of $395$ habitat patches found within the
black rectangle. $(b)$ Emergence of isolated (sub)communities (shaded areas)
separated by empty patches (white nodes) for $\beta=0.08$ in the random
geometric graph corresponding to the subset of habitat patches (see the main
text for the construction of this graph). $(c)$ The same but for $\beta=0.14$.
In this case, there is a single large community in addition to a few smaller
ones. $(d)$ Fraction of empty nodes as a function of $\beta$ with the min-max
deviation. The inset highlights the core-periphery structure of the graph
showing the adjacency matrix.

Let us mention that other species were also shown to live in a fragmented
habitat and to disperse between the patches. This is for instance the case of
African wild dogs in South Africa [14] or the giant kelp _Macrocystis
pyrifera_ in southern California [24, 15].

##  V Discussion and Conclusions

In this paper, we have introduced a self-segregation process and have shown
that it can induce a (meta)population network to split in (sub)communities of
contiguously occupied patches, separated by empty ones. We base our analysis
on analytical study and numerical validation using both random geometric
networks and empirical ecological ones. The stationary process reveals the
role of the network degree heterogeneity in the emergence of vacant nodes as
novel states whose fraction is quantised during the node filling process and
exhibits multistable configurations. Such features are absent in other random
processes such as random walks [37, 38, 51, 39, 52, 53] or Schelling’s
segregation [28, 29, 30, 54, 55]. From the ecological perspective, the
outcomes of our model stress that rich patterns endowed with vacant patches of
the habitat can exist as a natural outcome of trade-off between first
principles such as positive intra- or inter-species interactions e.g.,
cooperation and negative ones e.g., competition. Such theoretical evidence
supports even more the vacant niche paradigm as a genuine consequence of
population dynamics [6]. Furthermore, the quantisation effect illustrates the
robustness of the occupancy configurations to changes in the density of
individuals. The emergence of multistable states demonstrates that segregation
can be affected by the randomness of the process, in line with the paradigm of
the fitness landscape [6]. On the other side, our model proposes a robust
mechanism for the fragmentation of the habitat in disconnected spatial
patterns. We emphasize that such phenomenon is driven by endogenous factors,
i.e., facilitations between individuals and most importantly by the structure
of the metapopulation network. Topological features of the spatial support
have been recently shown to enable the patterning dynamics in minimal models
of single species systems [56]. Habitat fragementation has recently raised as
a strong candidate for the speciations of animals or plants [21, 22, 23]. In
particular, according to the peripatric speciation theory, the segregation of
individuals in isolated macroscopic patches can promote the development of the
new species as a consequence of isolated evolution [57].

The model we have proposed here focuses exclusively on a dispersal process
which puts the emphasis on the spatial interactions between individuals in
adjacent patches of the habitat rather than the local intrinsic dynamics of
the species considered here, thus deliberatively neglecting the birth/death
process or other intra- or inter-species interactions such as competition or
mutualism. Such considerations hold true when the local and global timescales
of the dynamics disentangle [31] to allow a fast relaxation at the system’s
equilibrium at the patch level, consequently allowing overall conservation of
mass by the slower segregation process. As an example and future potential
extension of our model, one can consider the family of Lotka-Volterra models,
which in the classic prey-predator case can produce periodic oscillations [5]
or fixed point when a density-dependent is introduced [58]. However, as
already emphasized, when such local dynamics tend to relax considerably faster
to the equilibrium compared to the tendency to migrate, the total mass is kept
constant on average, thus imposing an intrinsic fixed global carrying
capacity.

The relevance of vacant spatial patches can extend beyond the ecological
setting and help better understand how segregation in human dynamics that has
energetically flourished in many post-industrialised societies can lead to the
spontaneous creation of urban prairies [59]. Similarly to what we have
obtained in this work, a drop in the population density, usually triggered by
deurbanization or high crime levels, has been responsible for the emergence of
abandoned urban land that has reverted to green space consequently. Instead of
taking advantage of the abundantly available cheap land, people have decided
to segregate into (sub)communities where the population can enhance their
fitness to better use of resources (e.g., infrastructure, public transport,
education, etc.) [59]. In this regard, we are confident that the segregation
paradigm we present here can naturally extend from microscopic (bacteria,
viruses, etc.) to macroscopic scale (e.g., human communities).

##  VI Acknowledgments

JFDK is supported by a FNRS Aspirant Fellowship under the Grant FC38477. Part
of the results were obtained using the computational resources provided by the
“Consortium des Equipements de Calcul Intensif” (CECI), funded by the Fonds de
la Recherche Scientifique de Belgique (F.R.S.-FNRS) under Grant No. 2.5020.11
and by the Walloon Region.

##  Appendix A The stability analysis of the fixed points

We start by first emphasising that for the simple connected network made of
$\Omega$ nodes the system of ODEs (7)

$\frac{\text{d}\rho_{i}}{\text{d}t}=\sum_{l=1}^{\Omega}\mathcal{L}_{il}\rho_{i}\rho_{l}\left[1-\rho_{i}-\dfrac{k_{l}}{k_{i}}(1-\rho_{l})\right],\,\forall
i,$ (7)

of the main text generalises to:

$\displaystyle\frac{\text{d}\rho_{i}}{\text{d}t}$
$\displaystyle=\sum_{l=1}^{\Omega}\mathcal{L}_{il}\Big{[}\rho_{l}g(\rho_{i})-\frac{k_{l}}{k_{i}}\rho_{i}g(\rho_{l})\Big{]}$
$\displaystyle=\sum_{l=1}^{\Omega}A_{il}\Big{[}\frac{\rho_{l}g(\rho_{i})}{k_{l}}-\frac{\rho_{i}g(\rho_{l})}{k_{i}}\Big{]},\,\forall
i,$ (8)

where we have kept a general form for the function $g(n_{i}/N)$. We assume $g$
to be a $C^{1}$ concave function on the open interval $(0,1)$ with
$g(0)=0=g(1)$ and $g(x)>0$ for all $x\in(0,1)$. Consequently, the stationary
density of node $v_{i}$ ($i=1,\cdots,\Omega$) is either $0$ or given by the
unique non-zero solution (provided it exists) of the equation:

$\frac{\rho_{i}^{*}}{C_{m}k_{i}}=g(\rho_{i}^{*}),$ (9)

with the constant $C_{m}$ such that $\sum_{l\in
M_{m}}\rho_{l}^{*}=\Omega_{m}\beta_{m}$ where the sum runs over the nodes
belonging to the same community $M_{m}$ as node $v_{i}$ and
$\Omega_{m}\beta_{m}$ is the total mass in this (sub)community. Fig. 4
represents graphically the solution of Eq. 9 with $\zeta=C_{m}k$, for two
distinct values of $k$. Theorem $1$ establishes general results for the local
stability of the fixed points.

Figure 4:  Geometrical interpretation of the stationary densities solution of
Eq. (9). We show the graph of a generic function $g$ satisfying the required
assumptions, together with two straight lines associated to two values of the
node degree. Either, there will be only one intersection point, corresponding
to $\rho^{*}=0$ or, for sufficiently large degrees, there can be two
intersection points, namely $\rho^{*}=0$ and $\rho^{*}>0$. In the latter case,
only the non-zero solution is stable, as shown in Theorem $1$.

###### Theorem 1

Let $\mathcal{G}$ be a simple connected and undirected network made of
$\Omega$ nodes $v_{1},v_{2},\cdots,v_{\Omega}$ such that their corresponding
degrees are sorted in increasing order, i.e. $k_{1}\leq k_{2}\leq\cdots\leq
k_{\Omega}$. Consider the following dynamical system:

$\frac{\text{d}\rho_{i}}{\text{d}t}=\sum_{j=1}^{\Omega}\mathcal{L}_{ij}\Big{[}\rho_{j}g(\rho_{i})-\frac{k_{j}}{k_{i}}\rho_{i}g(\rho_{j})\Big{]},$
(10)

with initial conditions $\rho_{i}(0)=\rho_{i}(0)$ and
$\sum_{i=1}^{\Omega}\rho_{i}(0)=\beta\Omega$, for some $\beta\in(0,1)$.
Suppose that $g$ is a $C^{1}$ concave function on the open interval $(0,1)$
with $g(0)=0=g(1)$ and $g(x)>0\,\forall x\in(0,1)$. Assume moreover that all
the non-empty nodes in the steady state form a single community. Then:

  1. 1.

The following fixed point:

$\left\\{\begin{array}[]{ll}\rho_{i}^{*}=0~{}~{}\forall
i:k_{i}<\frac{1}{cg^{\prime}(0)}\\\ \rho_{i}^{*}>0~{}~{}\forall
i:k_{i}>\frac{1}{cg^{\prime}(0)}\end{array}\right.,$ (11)

with $\rho_{i}^{*}>0$ the non-zero solution of (9) is a locally stable fixed
point of the above dynamical system, with the constant $c$ given by:

$c=\frac{\beta\Omega}{\sum_{j=1}^{\Omega}k_{j}g(\rho_{j}^{*})}\,,$ (12)

and it coincides with $C_{1}$ defined from (9).

  2. 2.

Any fixed point with $\rho_{i}^{*}=0$ for some
$k_{i}>\frac{1}{cg^{\prime}(0)}$ is locally unstable.

The local stability analysis will be carried out by computing the jacobian
matrix $J$ associated to the previous system:

$\displaystyle J_{ij}$
$\displaystyle=\frac{\partial}{\partial\rho_{j}}\left(\frac{\text{d}\rho_{i}}{\text{d}t}\right)$
$\displaystyle=\sum_{l=1}^{\Omega}\frac{A_{il}}{k_{l}}\left[\delta_{lj}g(\rho_{i})+\rho_{l}g^{\prime}(\rho_{i})\delta_{ij}-\delta_{ij}g(\rho_{l})-\rho_{i}g^{\prime}(\rho_{l})\delta_{ij}\right]$
$\displaystyle=A_{ij}\left[\frac{g(\rho_{i})}{k_{j}}-\frac{\rho_{i}g^{\prime}(\rho_{j})}{k_{i}}\right]+\delta_{ij}\sum_{l}A_{il}\left[\frac{\rho_{l}g^{\prime}(\rho_{i})}{k_{l}}-\frac{g(\rho_{l})}{k_{i}}\right].$

Let us remember that the nodes are ranked in increasing order of their
degrees, i.e. $k_{1}\leq k_{2}\leq k_{\Omega}$. Let
$s\in\\{1,2,\cdots,\Omega\\}$ be such that:

$\displaystyle\left\\{\begin{array}[]{ll}k_{i}cg^{\prime}(0)<1{}&\text{for
}i=1,\cdots,s\\\ k_{i}cg^{\prime}(0)>1{}&\text{for
}i=s+1,\cdots,{\Omega}\end{array}\right..$ (15)

Let us now evaluate the jacobian matrix at the fixed point given by Eq. (11).
We will denote by $M$ the jacobian matrix $J$ evaluated at this fixed point.
The fixed point will be locally stable if and only if all the eigenvalues of
$M$ are negative.

  1. 1.

For $i,j=1,\cdots,s$, we have $\rho_{i}^{*}=0$ and $\rho_{j}^{*}=0$ and thus:

$\displaystyle M_{ij}$
$\displaystyle=\delta_{ij}\sum_{l=s+1}^{\Omega}A_{il}\left[\frac{\rho_{l}^{*}g^{\prime}(0)}{k_{l}}-\frac{g(\rho_{l}^{*})}{k_{i}}\right]$
$\displaystyle=\delta_{ij}\sum_{l=s+1}^{\Omega}\frac{A_{il}\rho_{l}^{*}}{k_{l}}\left[g^{\prime}(0)-\frac{1}{ck_{i}}\right].$
(16)

Since $k_{i}cg^{\prime}(0)<1$ for $i=1\cdots,s$, we deduce that $M_{ii}<0$ for
$i=1,\cdots,s$. We also deduced that any fixed point with $\rho_{i}^{*}=0$ and
$k_{i}cg^{\prime}(0)>1$ will be unstable.

  2. 2.

For $i=1,\cdots,s$ and $j=s+1,\cdots,\Omega$, we have $M_{ij}=0$. Consequently
$M$ is a block matrix and it remains to investigate the stability of the
submatrix $(M_{ij})_{ij}$ with $i,j=s+1\cdots,\Omega$.

  3. 3.

For $i,j=s+1,\cdots,\Omega$, we have :

$\displaystyle M_{ij}$
$\displaystyle=\frac{A_{ij}\rho_{i}^{*}}{k_{i}}\left[\frac{1}{ck_{j}}-g^{\prime}(\rho_{j}^{*})\right]+$
$\displaystyle+\delta_{ij}\sum_{l=s+1}^{\Omega}\frac{A_{il}\rho_{l}^{*}}{k_{l}}\left[g^{\prime}(\rho_{i}^{*})-\frac{1}{ck_{i}}\right].$
(17)

Let us observe that:

$\displaystyle\sum_{i=s+1}^{\Omega}M_{ij}$
$\displaystyle=\sum_{i=s+1}^{\Omega}\frac{A_{ij}\rho_{i}^{*}}{k_{i}}\left[\frac{1}{ck_{j}}-g^{\prime}(\rho_{j}^{*})\right]+$
$\displaystyle+\sum_{l=s+1}^{\Omega}\frac{A_{jl}\rho_{l}^{*}}{k_{l}}\left[g^{\prime}(\rho_{j}^{*})-\frac{1}{ck_{j}}\right]=0,$
(18)

since the network is undirected. Moreover, $\sum_{\begin{subarray}{c}i=s+1\\\
i\neq j\end{subarray}}^{\Omega}M_{ij}=-M_{jj}$. Assuming
$g^{\prime}(\rho_{i}^{*})-\frac{1}{ck_{i}}<0~{}\forall i=s+1\cdots,\Omega$,
one can make use of the Gershgorin theorem to deduce the stability of the
submatrix $(M_{i,j})_{ij=s+1}^{\Omega}$. Indeed, according to this theorem,
the spectrum of this submatrix is contained in the set of disks

$\bigcup\limits_{j=s+1}^{\Omega}D(M_{jj},R_{j})$ (19)

centred in $M_{jj}<0$ and with radius

$R_{j}=\sum_{\begin{subarray}{c}i=s+1\\\ i\neq
j\end{subarray}}^{\Omega}|M_{ij}|=-M_{jj}.$ (20)

Consequently, we deduce that the spectrum of the submatrix
$(M_{i,j})_{ij=s+1}^{\Omega}$ is contained in the left half plane, hence
showing that all the eigenvalues of the matrix $M$ have a negative real part.

##  Appendix B Detailed analysis of metapopulations breaking

In this Appendix, we further investigate the role of the network topology on
the emergence of (sub)communities. Let us recall that the dynamical system is
given by:

$\frac{\text{d}\rho_{i}}{\text{d}t}=\sum_{j=1}^{\Omega}\mathcal{L}_{ij}\rho_{i}\rho_{j}\left[1-\rho_{i}-\dfrac{k_{j}}{k_{i}}(1-\rho_{j})\right]\,,\,\forall
i=1,\dots,\Omega\,,$

with $g(x)=x(1-x)$. For sufficiently large values of the average density
$\beta$ and generic initial conditions without empty nodes, there will be no
asymptotically empty nodes, in which case the stationary nodes densities are
given by

$\rho_{i}^{*}=1-\frac{1-\beta}{k_{i}\langle 1/k\rangle},$

where $\langle 1/k\rangle=\sum_{j=1}^{\Omega}\frac{1}{k_{j}}$. Empty nodes
emerge as soon as there is at least a node of degree $k_{i}$ such that
$k_{i}{=}\frac{1-\beta}{\langle 1/k\rangle}$. Hence, the first node(s) to
become empty are the ones with the lowest degree $k_{\mathrm{min}}$ and this
happens when $\beta=1-k_{\mathrm{min}}\langle 1/k\rangle$. As an example let
us consider the network model of Fig. 1 in the main text that we here
repropose as Fig. 5 for which $k_{\mathrm{min}}=2$. There are three nodes with
degree $2$ whose density could potentially vanish for some $\hat{\beta}$. A
simple computation returns $\hat{\beta}=1-\frac{20}{27}\sim 0.259$.

Figure 5: As the reaction-diffusion process takes place, the first node(s) to
become empty is (are) the one(s) with the lowest degree, namely
$k_{\mathrm{min}}=2$, for the above graph.

To better infer the role of network topology on both the metapopulation
breaking and the quantisation of the fraction of vacancies, in Fig. 6 we
analyse a set of both synthetic and real networks with different structural
characteristics. In panel $a)$, we consider Erdős-Rényi graphs and show that
as the size of the random networks increases, the quantisation effect in
occupying, respectively depleting the nodes becomes more diluted. The
importance of the heterogeneity of the degree distribution in the formation of
(sub)communities is further accentuated in diassortative networks where nodes
are naturally arranged in interacting groups of highly and loosely connected
entities. In panel $b)$, we consider a synthetic and a real core-periphery
network, where at variance with the Erdős-Rényi topology (of panel $a)$),
besides quantisation, we observe the emergence of many different
(sub)communities (panels $d)$ and $e)$). This scenario is a neat manifestation
of multistability, i.e., different stable states for different initial
conditions for the same system parameter $\beta$, as attested by the shaded
area in panel $b)$. In fact, for two initial nodes densities conditions, with
same total mass, randomly chosen from a uniform distribution, we have two
distinct configurations of (sub)communities with a different fraction of
vacant nodes. There is a clear difference between the synthetic and the
empirical network since the latter is provided with “bridge” nodes that
connect two or more cores. Let us recall that the bridgeness centrality of
node $v_{i}$ is the fraction of shortest paths starting and ending not in the
neighbourhood of node $v_{i}$ and running through it [43]. This peculiarity
makes the real core-periphery network prone to split in different
(sub)communities once the density $\beta$ starts to decrease, at variance with
the synthetic network where multiple (sub)communities appear only when the
density takes low values. In panels $c)$ and $f)$, we consider the USAir97
network, an undirected network whose links represent the US air flights in
$1997$ [60]. The network was shown to be a core-periphery network with a
single densely connected core [61]. As the parameter $\beta$ decreases, most
of the mass concentrates in the core of the network, as shown in panel $f)$,
leaving a single large community of occupied nodes. Such difference between
the behaviors arising from the two networks can be justified by the presence
of bridge nodes that allow a broader repartition of the mass $\beta_{m}$.

Figure 6:  (Upper panels) Fraction of empty nodes $f$ vs. the average density
$\beta$. The shaded blue and red areas (in panels b) and c)) show the
multistability phenomenon as measured by the min-max node occupancy resulting
from several independent simulations. a) Erdos-Renyi networks with different
network sizes and same probability to have a link among two nodes ($p=1/3$).
As the number of nodes increases, the quantisation becomes weaker while the
multistability is almost absent. b) The London Tube network [62] (red curve)
vs. a synthetic core-periphery one (blue curve). The shaded area corroborates
the multistability claim, and the quantisation appears neatly manifested by
the plateaus. (Inset) Number of (sub)communities $N$ vs. the fraction of empty
nodes $f$ showing that to the same fraction of empty nodes corresponds a
different number of (sub)communities. c) The disassortative USAir97 network
[60]. (Lower panels) (Sub)communities visualisation (shaded coloured areas)
and empty nodes (white circles): d) and e) results for two different densities
initialisations for the same value of $\beta=0.15$ close to the abrupt
transition in $f$ vs $\beta$ (see panel b)) in the London Tube network. f) the
same result of panels d) and e) for the USAir97 network for $\beta=0.15$.

To further elucidate this assertion, in Fig. 7, we analyse the impact of the
network modularity on the fraction of empty nodes $f$. The red curve shows $f$
vs. $\beta$ for a disassortative modular network made of $20$ densely
connected subgraphs formed by a number of nodes comprised between $5$ and
$10$, linked together through bridge nodes, i.e. nodes whose removal splits
the network into disconnected subgraphs. Such networks can be obtained
starting from a ring network and then replacing every two nodes by a densely
connected subgraph. The corresponding adjacency matrix is shown in the inset
on the lower left corner and confirms the modularity of the network (each
pixel represents a connection between $2$ nodes). As the average density
$\beta$ decreases, we observe that the mean fraction of empty nodes increases
in a continuous way up to reaching a first plateau corresponding to the case
where (almost all the) bridge nodes are empty and the network is split into
(almost) $20$ subcommunities. For smaller values of $\beta$ the fraction of
empty nodes increases again until it stabilises for low values of $\beta$,
i.e. smaller than $\sim 10^{-2}$, when also nodes into the subgraphs are
empty. The blue curve shows $f$ vs. $\beta$ for a network obtained using a
rewiring degree preserving process on the former one. The new adjacency matrix
is now represented in the inset on the upper right corner. Due to the rewiring
procedure, the network now looks like a random network but the degree sequence
is unchanged. The corresponding behaviour of $f$ vs. $\beta$ is now quite
different and exhibits a step-like behaviour. In both cases, the results were
averaged over $100$ distinct network realisations. The min-max deviation in
the resulting value of $f$ is represented by the dashed region.

Figure 7: Fraction of empty nodes $f$ vs. the average density $\beta$. The
results were averaged over $100$ distinct network realisations. The shaded
area shows the multistability as measured by the min-max occupancy resulting
from all the realisations.

To understand these different behaviours we analyse the number of
(sub)communities emerging in the modular network as well as in the rewired one
(see Fig. 8). For a very large range of values of $\beta$ the rewired network
exhibits one (sub)community (see the blue curve in Fig. 8), associated to a
unique constant $C_{1}\equiv C_{1}(\beta)$. Consequently, nodes densities are
essentially constrained by their degree, as shown in Fig. 9 where we reported
the stationary nodes densities $\rho_{i}^{*}$ vs. their degree (nodes
densities were initialized randomly such that the average density satisfies
$\beta=0.15$). As a consequence, nodes with a degree larger than $C_{1}$ are
filled while all the others become empty. Since the degree sequence is finite,
this implies the observed jumps on the blue curve given in Fig. 7. On the
other hand, for the diassortative modular network, the mass first distributes
among the distinct modules ($20$ in this case) leaving the bridges nodes empty
(see the red curve in Fig. 8). To each of these modules corresponds a constant
$C_{m}=\left(1-\beta_{m}\right)/\langle 1/k\rangle_{\Omega_{m}}$. Inside the
$m$-th module, all the nodes with a degree $k_{i}<C_{m}$ will then become
empty. As $C_{m}$ is of order $C_{1}/l$, with $l$ the number of modules, nodes
with smaller degree can also be filled, leading to a lower value of the
fraction of empty nodes, except when $\beta$ is very large or very small (see
Fig. 7). Moreover, as the constants $C_{m}$ differ from each other, distinct
stationary densities can be observed for nodes with the same degree (see Fig.
10). Since the way the mass distributes among the modules depends (in a highly
non-trivial way) on the initial densities, averaging over sufficiently many
configurations leads in fine to a smooth evolution of $f$ vs. $\beta$ (see the
red curve in Fig. 7).

Figure 8: Number $N$ of (sub)communities vs. the fraction of empty nodes $f$.
The shaded area shows the multistability as measured by the min-max deviation
resulting from all the realisations. Figure 9: Stationary nodes densities
$\rho_{i}^{*}$ as a function of their degree $k$, for the rewired version of
the modular disassortative network used in Fig. 7. The initial conditions are
given by $\rho_{i}(0)=0.15$. Nodes densities were initialized randomly such
that the average density satisfies $\beta=0.15$. Figure 10: Stationary nodes
densities $\rho_{i}^{*}$ as a function of their degree $k$, for the modular
disassortative network used in Fig. 7. Nodes densities were initialized
randomly such that the average density satisfies $\beta=0.15$.

To highlight the impact of node bridgeness in the emergence of
(sub)communities, we show in Fig. 12 the results obtained using a small
modular network with a varying fraction of bridge nodes of degree $2$ between
the modules and we computed the number of (sub)communities $N$ vs. $\beta$. To
construct such network, we first generate $m$ cliques of size $k$ ($k\geq 3$).
We then connect the $i$-th clique to the $(i+1)$-th clique ($i=1,\cdots,m-1$)
by means of a path of length $1$ or $2$. The procedure is illustrated in Fig.
11, for $m=4$ and $k=4,5$. As $\beta$ is reduced, the bridge nodes become
empty, leaving the mass in the disconnected modules. The larger the number of
bridge nodes, the larger the number of (sub)communities.

Figure 11: Modular network with cliques of size $k=4,5$ connected through
paths of length $1$ or $2$. The removal of a (blue) node belonging to a path
of length $2$ breaks the connectivity of the network.  Figure 12: Number of
(sub)communities $N$ as a function of the average node density $\beta$, for a
modular network with a number $i$ of bridge nodes. The network was built using
the procedure described in the text with $10$ cliques of size $4\leq k\leq 7$.
For each value $i$, the data were averaged over $10$ independent initial
configurations.

We end up this Appendix by considering random geometric graphs obtained by
starting with a set of uniformly distributed points in the unit square and
connecting them if their euclidean distance is less than a given thresold $r$.
The larger $r$ and the smaller the fraction of empty nodes as shown in Fig.
13. In particular, for $r=\sqrt{2}$, we obtain a complete graph, for which
there is no empty node.

Figure 13: Fraction of empty nodes in random geometric graphs obtained by
connecting every pair of (uniformly distributed) points whose euclidean
distance is smaller than a given thresold $r$. The larger $r$ and the smaller
$f$.

## References

  * Newman [2017] M. E. J. Newman,  _Networks : an introduction_ (Oxford University Press, 2017).
  * Boccaletti  _et al._ [2006] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D. Hwang, Physics Reports 424, 175 (2006).
  * Porter and Gleeson [2016] M. A. Porter and J. P. Gleeson,  _Dynamical systems on networks : a tutorial_ (Springer, 2016).
  * Barrat  _et al._ [2008] A. Barrat, M. Barthelemy, and A. Vespignani,  _Dynamical processes on complex networks_ (Cambridge university press, 2008).
  * Murray [2002] J. D. Murray,  _Mathematical biology_ , 3rd ed., Interdisciplinary applied mathematics (Springer, New York, 2002).
  * Kauffman [1993] S. A. Kauffman,  _The origins of order_ (Oxford University Press, 1993).
  * MacArthur and Wilson [2001] R. H. MacArthur and E. O. Wilson,  _The theory of island biogeography_ (Princeton University Press, Princeton, 2001).
  * Hubbell [2001] S. P. Hubbell,  _The unified neutral theory of biodiversity and biogeography_ , Monographs in population biology No. 32 (Princeton University Press, Princeton, 2001).
  * Caswell [1976] H. Caswell, Ecological Monographs 46, 327 (1976).
  * Azaele  _et al._ [2016] S. Azaele, S. Suweis, J. Grilli, I. Volkov, J. R. Banavar, and A. Maritan, Reviews of Modern Physics 88, 035003 (2016).
  * Alonso  _et al._ [2006] D. Alonso, R. Etienne, and A. Mckane, Trends in Ecology & Evolution 21, 451 (2006).
  * Black and McKane [2012] A. J. Black and A. J. McKane, Trends in Ecology & Evolution 27, 337 (2012).
  * Leibold  _et al._ [2004] M. A. Leibold, M. Holyoak, N. Mouquet, P. Amarasekare, J. M. Chase, M. F. Hoopes, R. D. Holt, J. B. Shurin, R. Law, D. Tilman, _et al._ , Ecology letters 7, 601 (2004).
  * Nicholson  _et al._ [2020] S. K. Nicholson, D. G. Marneweck, P. A. Lindsey, K. Marnewick, and H. T. Davies-Mostert, African Journal of Wildlife Research 50, 8 (2020).
  * Castorani  _et al._ [2015] M. C. Castorani, D. C. Reed, F. Alberto, T. W. Bell, R. D. Simons, K. C. Cavanaugh, D. A. Siegel, and P. T. Raimondi, Ecology 96, 3141 (2015).
  * Hanski  _et al._ [1994] I. Hanski, M. Kuussaari, and M. Nieminen, Ecology 75, 747 (1994).
  * Orsini  _et al._ [2008] L. Orsini, J. Corander, A. Alasentie, and I. Hanski, Molecular ecology 17, 2629 (2008).
  * Hanski [2011] I. A. Hanski, Proceedings of the National Academy of Sciences 108, 14397 (2011).
  * Hanski  _et al._ [2017] I. Hanski, T. Schulz, S. C. Wong, V. Ahola, A. Ruokolainen, and S. P. Ojanen, Nature Communications 8, 1 (2017).
  * Opedal  _et al._ [2020] Ø. H. Opedal, O. Ovaskainen, M. Saastamoinen, A.-L. Laine, and S. van Nouhuys, Ecology 101, e03186 (2020).
  * Weir and Schluter [2004] J. T. Weir and D. Schluter, Proceedings of the Royal Society of London. Series B: Biological Sciences 271, 1881 (2004).
  * Dias  _et al._ [2013] M. S. Dias, J.-F. Cornu, T. Oberdorff, C. A. Lasso, and P. A. Tedesco, Ecography 36, 683 (2013).
  * Thompson  _et al._ [2018] K. A. Thompson, L. H. Rieseberg, and D. Schluter, Trends in Ecology & Evolution 33, 815 (2018).
  * Bell  _et al._ [2015] T. W. Bell, K. C. Cavanaugh, and D. A. Siegel, Remote Sensing of Environment 167, 218 (2015).
  * Grinnell [1917] J. Grinnell, The Auk 34, 427 (1917).
  * Goodnight [2011] C. J. Goodnight, Philosophical Transactions of the Royal Society B: Biological Sciences 366, 1401 (2011).
  * Wilson [1992] D. S. Wilson, Ecology 73, 1984 (1992).
  * Schelling [1969] T. C. Schelling, The American Economic Review 59, 488 (1969).
  * Schelling [1971] T. C. Schelling, The Journal of Mathematical Sociology 1, 143 (1971).
  * Rogers and McKane [2011] T. Rogers and A. J. McKane, Journal of Statistical Mechanics: Theory and Experiment 2011, P07006 (2011).
  * Asllani  _et al._ [2020] M. Asllani, B. R. Da Cunha, E. Estrada, and J. P. Gleeson, New Journal of Physics 22, 063037 (2020).
  * Allee and Bowen [1932] W. C. Allee and E. S. Bowen, Journal of Experimental Zoology 61, 185 (1932).
  * Jørgensen and Fath [2008] S. E. Jørgensen and B. D. Fath, eds.,  _Encyclopedia of ecology_ , 1st ed. (Elsevier, Amsterdam ; Boston, 2008) oCLC: ocn173240026.
  * Hixon [1987] M. A. Hixon, American Zoologist 27, 229 (1987).
  * Gardiner [2004] C. W. Gardiner,  _Handbook of Stochastic Methods_, edited by H. Haken, Springer Series in Synergetics, Vol. 13 (Springer Berlin Heidelberg, Berlin, Heidelberg, 2004).
  * Kampen [2007] N. G. v. Kampen,  _Stochastic processes in physics and chemistry_ , 3rd ed., North-Holland personal library (Elsevier, Amsterdam ; Boston, 2007) oCLC: ocm81453662.
  * Asllani  _et al._ [2018] M. Asllani, T. Carletti, F. Di Patti, D. Fanelli, and F. Piazza, Physical Review Letters 120, 158301 (2018).
  * Carletti  _et al._ [2020] T. Carletti, M. Asllani, D. Fanelli, and V. Latora, Physical Review Research 2, 033012 (2020).
  * Fronczak and Fronczak [2009] A. Fronczak and P. Fronczak, Physical Review E 80, 016107 (2009).
  * Note [1] In a biased random walk, the motion of individuals depends on the relative attribute of the hosting node compared to all the neighbour ones [39, 51].
  * Fanelli and McKane [2010] D. Fanelli and A. J. McKane, Physical Review E 82, 021113 (2010).
  * Note [2] Based on van Kampen ansatz, we can write the (rescaled) discrete variable as $n_{i}/N=\rho_{i}+\xi_{i}/\sqrt{N}$, valid for large $N$, where $\rho_{i}$ is the deterministic variable describing the density of node $v_{i}$ and $\xi_{i}$ the stochastic variable. Thus $\langle n_{i}n_{j}\rangle/N^{2}=\langle\rho_{i}\rho_{j}+\xi_{i}\xi_{j}/N+\left(\xi_{j}\rho_{i}+\xi_{i}\rho_{j}\right)/\sqrt{N}\rangle\sim\langle n_{i}\rangle\langle n_{j}\rangle/N^{2}$ in the limit for large $N$ while the densities are kept constant.
  * Jensen  _et al._ [2016] P. Jensen, M. Morini, M. Karsai, T. Venturini, A. Vespignani, M. Jacomy, J.-P. Cointet, P. Mercklé, and E. Fleury, Journal of Complex Networks 4, 319 (2016).
  * Note [3] Obviously, an isolated, empty node (surrounded by other empty ones) is always unstable since the mass added to it will never leave the node.
  * Note [4] This can be easily obtained from the formula $C_{m}=\left(1-\beta_{m}\right)k$ since the average $\langle\cdot\rangle_{\Omega}$ over the degrees drops being the graph regular. From here, it is immediate that $C_{m}/k=\left(1-\beta_{m}\right)<1$, which justifies the stability of the fully occupied nodes state.
  * Grilli  _et al._ [2015] J. Grilli, G. Barabás, and S. Allesina, PLoS computational biology 11, e1004251 (2015).
  * Gross  _et al._ [2020] T. Gross, K. T. Allhoff, B. Blasius, U. Brose, B. Drossel, A. K. Fahimipour, C. Guill, J. D. Yeakel, and F. Zeng, Philosophical Transactions of the Royal Society B 375, 20190455 (2020).
  * Ryser  _et al._ [2019] R. Ryser, J. Häussler, M. Stark, U. Brose, B. C. Rall, and C. Guill, Proceedings of the royal society B 286, 20191177 (2019).
  * Hijmans and University of California [2015] R. J. Hijmans and B. M. o. V. Z. University of California, Boundary, Åland islands, https://geodata.lib.utexas.edu/catalog/stanford-vr234jm5375 (2015).
  * Øystein Opedal [2020] Øystein Opedal, oysteiop/alandinteractions, Publication release (1.0). Zenodo. https://doi.org/10.5281/zenodo.3956435 (2020).
  * Gómez-Gardeñes and Latora [2008] J. Gómez-Gardeñes and V. Latora, Physical Review E 78, 065102 (2008).
  * Lambiotte  _et al._ [2014] R. Lambiotte, J.-C. Delvenne, and M. Barahona, IEEE Transactions on Network Science and Engineering 1, 76 (2014).
  * Masuda  _et al._ [2017] N. Masuda, M. A. Porter, and R. Lambiotte, Physics Reports 716-717, 1 (2017).
  * Gandica  _et al._ [2016] Y. Gandica, F. Gargiulo, and T. Carletti, Chaos, Solitons & Fractals 90, 46 (2016).
  * Gargiulo  _et al._ [2017] F. Gargiulo, Y. Gandica, and T. Carletti, Advances in Complex Systems 20, 1750001 (2017).
  * Asllani and Carletti [2018] M. Asllani and T. Carletti, Physical Review E 97, 042302 (2018).
  * Coyne and Orr [2004] J. A. Coyne and H. A. Orr,  _Speciation_ (Sinauer Associates, Sunderland, Mass, 2004) oCLC: ocm55078441.
  * Kot [2001] M. Kot,  _Elements of mathematical ecology_ , edited by Cambridge (Cambridge University Press, 2001).
  * Prener  _et al._ [2020] C. G. Prener, T. H. Braswell, and D. J. Monti, Journal of Urban Affairs 42, 371 (2020).
  * Rossi and Ahmed [2015] R. A. Rossi and N. K. Ahmed, in  _AAAI_ (2015).
  * Fasino and Rinaldi [2020] D. Fasino and F. Rinaldi, Symmetry 12, 94 (2020).
  * De Domenico  _et al._ [2014] M. De Domenico, A. Sole-Ribalta, S. Gomez, and A. Arenas, Proceedings of the National Academy of Sciences 111, 8351 (2014).