# Bifurcations and multistability in empirical mutualistic networks

Andrus Giraldo<EMAIL_ADDRESS>Deok-Sun Lee<EMAIL_ADDRESS>School
of Computational Sciences, Korea Institute for Advanced Study, Seoul 02455,
Korea

###### Abstract

Individual species may experience diverse outcomes, from prosperity to
extinction, in an ecological community subject to external and internal
variations. Despite the wealth of theoretical results derived from random
matrix ensembles, a theoretical framework still remains to be developed to
understand species-level dynamical heterogeneity within a given community,
hampering real-world ecosystems’ theoretical assessment and management. Here,
we consider empirical plant-pollinator mutualistic networks, additionally
including all-to-all intragroup competition, where species abundance evolves
under a Lotka-Volterra-type equation. Setting the strengths of competition and
mutualism to be uniform, we investigate how individual species persist or go
extinct under varying the interaction strengths. By employing bifurcation
theory in tandem with numerical continuation, we elucidate transcritical
bifurcations underlying species extinction and demonstrate that the Hopf
bifurcation of unfeasible equilibria and degenerate transcritical bifurcations
give rise to multistability, i.e., the coexistence of multiple attracting
feasible equilibria. These bifurcations allow us to partition the parameter
space into different regimes, each with distinct sets of extinct species,
offering insights into how interspecific interactions generate one or multiple
extinction scenarios within an ecological network.

††preprint: APS/123-QED

##  I Introduction

The stability of biodiversity and species abundance distribution in each
ecological community - whether maintaining current states or undergoing abrupt
changes, including the extinction of some or almost all species - is a crucial
question in ecological conservation efforts and also has been of central
interest to theoretical approaches in ecology [1, 2, 3]. Insights have been
gained from studying the properties of equilibrium states in an ensemble of
fully connected communities with random interaction strength, utilizing tools
e.g., from the spin glass physics and random matrix theory [4, 5,
*Allesina2015, 7, 8, 9, 10]. Yet, real-world communities are often sparsely
connected and heterogeneously structured [11, 12, 13, 14, 15], demanding
further investigations to adapt these theoretical frameworks to empirical
features [16, 17]. Moreover, while ensemble-averaged properties provide broad
insights, they often fall short of predicting specific future scenarios for
individual communities, which are crucial for effective assessment and
management. This gap underscores the importance of developing new frameworks
that are better suited to understand and predict the dynamics of specific
communities, based on their actual structures.

Under diverse environmental and internal perturbations, individual species can
experience variations in their abundances or face extinction, and various
theoretical approaches have been proposed to explore these dynamics [18, 19,
20, 21]. Simple dynamical models, such as Lotka-Volterra(LV)-type equations,
have been widely used because they effectively capture the evolution of
species abundance under nonlinear interspecific interactions [22, 23, 24, 25,
3, 26, 27, 28]. Recent empirical data-sets from real-world ecological networks
[29, 30, 31, 32] have spurred investigations into such model dynamics on these
networks, offering insights into the interplay between structure and dynamics
[33, 34, 35, 36].

Identifying equilibria, periodic solutions, and their stability from the
nonlinear dynamics in specific real-world networks can enable predictions of
possible specific scenarios for their biodiversity and species composition
[37]. This may be a key next-generation goal in the field of nonlinear
dynamics, especially in the era of big data. However, it can be challenging
because real-world systems are typically large and their dynamics are high
dimensional, making mathematical techniques of lower dimensional systems
harder to apply[38, 39, 40]. To proceed, various approximate methods have been
attempted, such as approximating the interspecific interaction matrix by a
low-rank matrix to gain insights into the attracting equilibrium, assumed to
be unique, in specific network structures [15, 36]. However, even small errors
in these predictions can lead to disastrous results when applied to real-world
ecosystems, making a methodology that is as exact as possible highly
desirable.

Given this background, we apply established methods from bifurcation theory
[41] and numerical continuation [42] to investigate the equilibria that are
both attracting and feasible, with the latter meaning no negative abundance,
and how they vary with system parameters in empirical plant-pollinator
networks governed by a LV-type equation. While keeping the real-world sparse
and heterogeneous structure of mutualistic partnership between two groups -
plant and pollinator - of species [11, 12, 13, 14, 15], we assume uniform
mutualism strength [2, 15, 36] and all-to-all uniform intragroup competition
that may arise from limited resources [13, 43].

We find that among exponentially many equilibria, only one is both feasible
and attracting unless interactions are sufficiently strong; As interaction
strengths increase, different equilibria successively become the only globally
attracting equilibrium through transcritical bifurcations each featuring no or
nonempty sets of extinct species. When interactions are sufficiently strong,
different initial conditions can lead to different equilibria, so a species
may survive in one but be extinct in another equilibrium. Our study reveals
that such multistability arises from two distinct types of bifurcations: Hopf
bifurcation and degenerate transcritical bifurcation. While feasible
equilibria cannot exhibit Hopf bifurcations in our studied networks due to
them having symmetric interaction matrices [44, 45], our study demonstrates
that the Hopf bifurcation of an unfeasible equilibrium, followed by a
transcritical bifurcation, can lead to the creation of a new feasible
attracting equilibrium. Also, degenerate transcritical bifurcations are shown
to create attracting equilibria with different components in a particular
number of surviving species from the original equilibrium, resulting in
similar biodiversity scenarios. In contrast, Hopf bifurcations can lead to
vastly different biodiversity scenarios. Our study thus elucidates the
geometrical mechanisms underlying species extinction and multistability in
empirical mutualistic networks. By identifying these bifurcations, we
determine the exact parameter regimes that display distinct sets of surviving
species. This can provide a platform for developing strategies to preserve and
control biodiversity in real ecological communities.

The paper is organized as follows. In Sec. II, we introduce the model and the
persistence diagrams illustrating the survival and extinction of individual
species and multistability in the empirical networks under study. In Sec. III,
we provide the mathematical properties of the model and the numerical
continuation methods to identify potential bifurcations. In Sec. IV.1, we
explore the transcritical bifurcations that underlie species extinction. The
two geometrical mechanisms for multistability are investigated in detail in
Sec. IV.2 and IV.3, respectively. In Sec. V, we present a phase diagram
illustrating how the bifurcations organize themselves as the strengths of
competition and mutualism are varied. Finally, we conclude in Sec. VI with a
summary of our findings and an outlook on open questions and future research.

##  II Model and overview

We consider the abundances $x_{1},x_{2},...,x_{s}$ of $s$ species, including
$n_{p}$ plant-group species and $n_{a}$ animal(pollinator)-group species. For
the remainder of the paper, we use indices $1$ to $n_{p}$ for plant species
and $n_{p}+1$ to $s=n_{p}+n_{a}$ for animals. Species interact via all-to-all
intragroup competition, and selective intergroup mutualism represented by the
mutualism adjacency matrix $A$ of dimensions $n_{p}\times n_{a}$ with
$A_{ij}=A_{ji}=1$ if plant species $i$ and animal species $n_{p}+j$ are in
mutualistic relation, or zero otherwise. These interspecific interactions
affect the species abundance evolving with time as

$\dfrac{d\mathbf{x}}{dt}=f(\mathbf{x})=\mathbb{X}\left(\boldsymbol{\alpha}+\mathbb{B}\mathbf{x}\right),$
(1)

where $\mathbf{x}:=(x_{1},x_{2},...,x_{s})$ is the vector of species
abundance, $\boldsymbol{\alpha}:=\alpha(1,1,...,1)$ is a vector of self-growth
rates with $\alpha$ a positive constant, and $\mathbb{X}:={\rm
diag}(x_{1},x_{2},...,x_{s})$ is a diagonal matrix with the species abundance.
The interaction matrix $\mathbb{B}$ represents the self-regulation, the
intragroup competition, and the intergroup mutualism by

$\displaystyle\mathbb{B}$
$\displaystyle:=-(1-c)\mathbb{I}-c\mathbb{J}+m\mathbb{A}$
$\displaystyle=\begin{pmatrix}-1&-c&\cdots&-c&&&&\\\ -c&-1&\cdots&-c&&&&\\\
\vdots&&&\vdots&&m\,A&&\\\ -c&\cdots&-1&-c&&&&\\\ -c&\cdots&-c&-1&&&&\\\
&&&&-1&-c&\cdots&-c\\\ &&&&-c&-1&\cdots&-c\\\ &m\,A^{T}&&&\vdots&&&\vdots\\\
&&&&-c&\cdots&-1&-c\\\ &&&&-c&\cdots&-c&-1\end{pmatrix},$ (2)

where $\mathbb{I}$ is the identity matrix of dimension $s\times s$, and

$\mathbb{J}:=\begin{pmatrix}\mathbf{J}_{n_{p}}&0\\\
0&\mathbf{J}_{n_{a}}\end{pmatrix}\quad{\rm
and}\quad\mathbb{A}:=\begin{pmatrix}0&A\\\ A^{T}&0\end{pmatrix}$

with $\mathbf{J}_{n_{p}}$ and $\mathbf{J}_{n_{a}}$ the all-ones matrices of
dimensions $n_{p}\times n_{p}$ and $n_{a}\times n_{s}$, respectively. Notice
that parameters $c$ and $m$ control the strength of competition and mutualism,
respectively. We set $\alpha=1$ as $\alpha$ only rescales time and abundance;
that is, if $\mathbf{x}(t)$ is a solution for a particular $\alpha$ then
$\alpha\mathbf{x}(\alpha t)$ is a solution to Eq. (1) with $\alpha=1$. Without
interspecific interactions, i.e., $c=m=0$, each species would be independent.

Figure 1: Abundance and persistence of individual species. (a1) Time evolution
of the abundances $x_{i}$’s of 107 species in network $A_{1}$ with $c=0.3$ and
$m=0.08$ for a randomly-selected initial condition. (a2) The same as (a1), but
with $m=0.25$. (a3) Persistence diagram, representing by horizontal lines,
green for plants and brown for animals, the ranges of $m$ for which species
survive when $c$ is fixed at $0.3$. This is obtained by numerical
continuation, starting from the full-coexistence equilibrium at $m=0$. (b1)
The same as (a2), but for a different initial condition. (b2) The same as
(b1), but with $m=0.485$. (b3) Persistence diagram starting from another
equilibrium at $m=m^{(2)}\approx 0.187267$. The species 53 represented by a
purple square has a negative abundance before a transcritical bifurcation
makes it zero at $m^{(2)}$. See Sec. IV.2. Vertical dashed lines in (a3) and
(b3) indicate the values of $m$ used in (a1),(a2), (b1), and (b2).

For the mutualism matrix $A$, we use two datasets from the Web of Life
database [29], one consisting of $43$ plants and $64$ pollinators, and the
other $29$ plants and $81$ pollinators, which we refer to as networks $A_{1}$
and $A_{2}$, respectively. The structure of these real-world mutualistic
communities is not uniform or random but structured such that they are often
nested [11, 46, 47] and individual species have widely different numbers of
mutualistic partners, resulting in different numbers of non-zero elements
across rows and columns in the matrix $A$ in Eq. (2) [13, 14]. The aim of the
present study is to investigate the equilibrium states of such structured
competitive-mutualistic networks.

Setting the initial abundance $x_{i}(0)$ of each species $i$ to be a random
number between $0$ and $1$, we numerically integrate system (1) with a fixed
competition strength $c=0.3$ and different values of the mutualism strength
$m$ to obtain $x_{i}(t)$’s for all $i$ and time $t$ as shown in Fig.
1(a1),(a2), (b1), and (b2). When the mutualistic interaction is sufficiently
weak, e.g., $m=0.08$ as in Fig. 1(a1), all the species coexist, which we call
the full-coexistence state. On the other hand, for $m=0.25$, only $17$ species
survive as shown in Fig. 1(a2); hence, strong mutualistic interactions may
break the full coexistence, resulting in some species becoming extinct while
others flourish.

The persistence of species under a gradual increase of $m$, modeling the
environmental or internal perturbations, exhibits a strong inequality; the
persistence line drawn in the range $[0,m^{*}_{F,i})$, for which each species
$i$ survives, is of quite different length. It is also remarkable that no line
extends infinitely, as the surviving species’ abundance grows without bounds
past the end point. These persistence lines are computed by starting a
numerical continuation scheme from the attracting equilibrium of system (1) at
$m=0$, and then studying its variation under increments of $m$; more details
are in the next section.

At large values of $m$, one can find a totally different equilibrium state, by
using different initial conditions, from that reached from the full-
coexistence state considered in Fig. 1(a3). This means that there can be
multiple attracting equilibria, for the same system parameters, each
representing a different scenario of species survival and extinction. For
instance, at $c=0.3$ and $m=0.3$, the time-evolution of individual species
abundance is different between Fig. 1(a2) and Fig. 1(b1). The final state in
Fig. 1(b2) at $m=0.485$ showcases how the new attracting equilibrium in Fig.
1(b1) exhibits different extinctions as $m$ increases compared to Fig. 1(a3).
In Fig. 1(b3), we show the persistence diagram of the new equilibrium, which
first appears at $m^{(2)}\approx 0.187267$ with some species already extinct
at that point.

The exact persistence diagrams in Fig. 1(a3) and (b3) extend the approximate
results in [36] for small values of $m$. In the subsequent sections, we
describe the numerical method used to generate these diagrams and investigate
the variation of the equilibrium states with $c$ and $m$ from the perspective
of bifurcation theory. The identified nature of the bifurcations underlying
species extinction and the appearance of multiple attracting equilibria
provides a geometrical picture of the equilibria and their stability.

##  III Numerical continuation method to trace the attracting and feasible
equilibria

The computational cost of integrating numerically system (1) with continuously
varying the system parameters, to obtain the persistence diagrams shown in
Fig. 1(a3) and (b3), is quite high; especially, if one would want a good
resolution in both $m$ and $c$. Instead, one can consider the equilibria ${\bf
x}^{*}$ of system (1) satisfying

$\mathbb{X}^{*}\left(\mathbf{1}_{s}+\mathbb{B}\mathbf{x}^{*}\right)=0$ (3)

and their stability. Once identifying an attracting equilibrium for particular
values of $c$ and $m$, we employ pseudo-arch length continuation techniques
[42] to study the variation in their components(abundances) and stability with
$c$ and $m$. In this way, we can pinpoint the bifurcations that the equilibria
might exhibit, allowing us to identify attracting equilibria for general
values of the system parameters without a lot of additional cost.

System (1) can generically have exponentially many equilibria given $2^{s}$
different scenarios of extinct species. With
$\mathcal{J}\subseteq\\{1,2,...,s\\}$ a subset of indices, the equilibrium
solution $\mathbf{x_{\mathcal{J}}^{*}}\in\mathbb{R}^{s}$ with the species in
$\mathcal{J}$ extinct is given by

$x_{\mathcal{J},j}^{*}=\left\\{\begin{array}[]{ll}0&\quad{\rm for}\quad
j\in\mathcal{J},\\\
-[\mathbb{B^{+}}^{-1}\mathbf{1}_{s-|\mathcal{J}|}]_{j}&\quad{\rm for}\quad
j\notin\mathcal{J},\end{array}\right.$ (4)

where $\mathbb{B^{+}}$ is the surviving species’ interaction matrix obtained
by eliminating in $\mathbb{B}$ the rows and columns that are in $\mathcal{J}$.
This equilibrium exists if and only if $\mathbb{B^{+}}$ is invertible;
otherwise, an equilibrium with only $\mathcal{J}$ species extinct does not
exist or there are infinitely many such equilibria. In the present study, it
should be noted that not all but only feasible and attracting equilibria are
physically meaningful. An equilibrium with $x_{i}\geq 0$ for all $i$ is
considered as feasible since species cannot have a negative abundance [2], and
as attracting, locally or globally, if certain or all initial conditions
converge to it. A major challenge is thus to identify feasible and attracting
equilibria and understand how they vary with the system parameters.

At $m=0$ for $0\leq c<1$, the full-coexistence equilibrium
$\mathbf{x_{\emptyset}^{*}}$ with no extinct species is globally attracting,
from Goh’s Theorem [24, 23], given that $\mathbb{B}$ is negative definite
[Appendix A]. As $\mathbb{B}$ is symmetric and $m\mathbb{A}$ can be seen as a
continuous symmetric perturbation made to $\mathbb{B}$ from $m=0$, the
eigenvalues of $\mathbb{B}$ vary continuously with $m$ and thus $\mathbb{B}$
remains negative definitive and $\mathbf{x_{\emptyset}^{*}}$ a global
attractor for sufficiently small $m$.

As $m$ becomes larger, two cases might arise: $\mathbb{B}$ stops being
negative definitive, or $\mathbf{x_{\emptyset}^{*}}$ finds one or more
components zero. In either case, Goh’s theorem cannot be applied anymore,
leading us to investigate the local attractivity of the full-coexistence or
other equilibria. If all the eigenvalues of the Jacobian
$Df|_{\mathbf{x^{*}_{\mathcal{J}}}}$ of system (1) at an equilibrium
$\mathbf{x^{*}_{\mathcal{J}}}$ have negative real parts, then
$\mathbf{x_{\mathcal{J}}}^{*}$ is stable (locally attracting). By studying the
variation of these eigenvalues as $m$ increases, one can find when the full-
coexistence equilibrium may become unstable and another equilibrium becomes
stable.

We apply numerical continuation techniques [48, 49] to trace changes in the
local stability of the equilibrium solutions satisfying Eq. (3) while
parameters are varied. As $m$ increases for fixed $c$, an attracting
equilibrium $\mathbf{x^{*}_{\mathcal{J}}}$ may become unstable at a certain
value $m_{*}$, which can be detected by monitoring if the determinant of the
Jacobian $Df|_{\mathbf{x^{*}_{\mathcal{J}}}}$ represented by [36]

$\det\left(Df|_{\mathbf{x^{*}_{\mathcal{J}}}}\right)=\prod_{i\in\mathcal{J}}\left(1+\sum_{\ell=1}^{s}B_{i\ell}x_{\mathcal{J},\ell}^{*}\right)\prod_{i\notin\mathcal{J}}\left(x^{*}_{\mathcal{J},i}\right)\det\left(\mathbb{B}^{+}\right)$
(5)

becomes zero. For system (1), we identify the following cases:

  1. 1.

The factor $\left(\prod_{i\notin\mathcal{J}}x^{*}_{\mathcal{J},i}\right)$
becomes zero, as one or more components not belonging to $\mathcal{J}$ become
zero, implying the extinction of the corresponding species, at $m_{*}$. If a
surviving species, say $j$, has zero abundance at $m_{*}$, this indicates a
generic transcritical bifurcation [41] as the equilibrium
$\mathbf{x^{*}_{\mathcal{J}}}$, which is feasible and attracting when $m$ is
close to but smaller than $m_{*}$, coincides with another equilibrium
$\mathbf{x}^{*}_{\mathcal{J}^{\prime}}$ with
$\mathcal{J}^{\prime}=\mathcal{J}\cup\\{j\\}$ at $m_{*}$. For $m>m_{*}$,
$\mathbf{x^{*}_{\mathcal{J}}}$ is unstable (also unfeasible since
$x^{*}_{\mathcal{J},j}<0$) and the new one
$\mathbf{x}^{*}_{\mathcal{J}^{\prime}}$ is attracting and feasible; the
species in $\mathcal{J}^{\prime}-\mathcal{J}$ are newly extinct 111If more
than one components, say the components in a set $\mathcal{K}$, become zero at
$m_{*}$, then the equilibrium $\mathbf{x^{*}_{\mathcal{J}}}$ coincides with
all the equilibria $\mathbf{x_{\mathcal{J}^{\prime}}^{*}}$ with all
$\mathcal{J}^{\prime}$ in $\mathcal{P}(\mathcal{K})$, where
$\mathcal{P}(\mathcal{K})$ is the power set of $\mathcal{K}$. This occurs if
Eq. (1) is equivariant [54] under a permutation of species, i.e., when a
collection of plants or animals have the same set of mutualistic partners..

  2. 2.

The factor $\det\left(\mathbb{B}^{+}\right)$ is zero. Then $\mathbb{B}^{+}$ is
not invertible, and the solution $\mathbf{x^{*}_{\mathcal{J}}}$ cannot be
given by Eq. (4); there are either infinitely many or no solutions at all. The
former corresponds to degenerate transcritical bifurcations, a mechanism
generating multistability detailed in Section IV.3. The latter case implies
that $\mathbf{x_{\mathcal{J}}^{*}}$ exhibits a transcritical bifurcation at
infinity [51, 52, 53]; Past this transition at $m_{*}$, one finds the new
equilibrium with abundances growing unbounded.

Increasing $m$ at fixed $c$, we compute which of the two cases occurs first.
If the first one occurs, then we consider the new attracting equilibrium
$\mathbf{x_{\mathcal{J}^{\prime}}^{*}}$ and check with it which of the two
cases occurs while increasing $m$ further. These procedures are repeated until
the second case occurs. Suppose that the second case occurs, and it is a
degenerate transcritical bifurcation. Then, we consider at least two different
equilibria, indicating multistability, and recursively study all the
attracting equilibria while increasing $m$ further. If the second case occurs
and there is a transcritical bifurcation at infinity, then we stop.

Since the considered networks are large, keeping track of all the equilibria
in system (1) and finding which one can create multistatibility is
challenging. To circumvent this problem, we numerically integrate system (1)
for particular values of $m$ and $c$ under different initial conditions and
investigate whether the solutions converge to more than one equilibria. If it
happens, then we apply our continuation routine to each of those distinct
equilibria in reverse, e.g., by decreasing $m$ at fixed $c$, to identify a
series of transcritical bifurcations leading to the attracting and feasible
equilibria with smaller sets of extinct species. Eventually an attracting but
unfeasible equilibrium is identified at $m^{\prime}>0$, which is found to have
been stabilized by Hopf bifurcation at $m$ smaller than $m^{\prime}$ as will
be detailed in Section IV.2. For each of these original equilibria, we produce
persistence diagram as shown in Fig. 1(a3) and (b3) to demonstrate
multistability.

These bifurcations disclose the geometrical mechanisms underlying species
extinction and multistability; Hopf bifurcations and degenerate transcritical
bifurcations lead to new attracting equilibria while transcritical
bifurcations transfer stability from one equilirium to another, maintaining
the number of attracting equilibria. The involved computations are performed
using the widely-used software package Auto07p [48, 49] that finds and traces
equilibria, periodic orbits, and suitable two-point boundary values problems
as functions of system parameters. Moreover, it allows for the identification
of different types of bifurcations, the monitoring of the stability of
equilibria, and the determination of when two or more solution branches meet
at a particular value and switch between them during continuation. The latter
is relevant when different equilibria coincide at the moment of a
transcritical bifurcation. A good introduction to such techniques and the use
of Auto07p can be found in [42].

##  IV Results

In this section, we present and analyze the results obtained by applying the
method described in Sec. III. Numerical continuation yields the variation of
an equilibrium, and the possible change of its stability as a parameter
varies, which informs us of the nature of bifurcations and allows us to
discover the geometrical mechanisms of species extinction and multistability.
The results of this section lead to the phase diagram in Sec. V.

###  IV.1 Transcritical bifurcations underlying species extinction

Figure 2: Transcritical bifurcations and species extinction. (a) The largest
eigenvalues of the equilibria $\mathbf{x^{*}_{\emptyset}}$ (blue line),
$\mathbf{x}^{*}_{\\{86\\}}$ (red) and $\mathbf{x}^{*}_{\\{86,76\\}}$ (purple)
as functions of $m$ for $c=0.3$ in network $A_{1}$. Their signs change at
$m_{\\{86\\}},m_{\\{86,76\\}}$ and $m_{\\{86,76,95,97\\}}$ represented by
triangles. (b) A magnification of (a), where shown are the two largest
eigenvalues of the equilibria $\mathbf{x}^{*}_{\\{86,76\\}}$ (purple lines)
and $\mathbf{x}^{*}_{\\{86,76,95,97\\}}$ (green lines). (c) Phase portrait in
the $(x_{95},x_{76},x_{86})$-space for $m<m_{\\{86\\}}$, where the full
coexistence equilibrium $\mathbf{x^{*}_{\emptyset}}$ (blue dot) is globally
attracting. Schematic trajectories (green lines) starting from different
initial conditions (black dots) are shown to converge to it; the arrows
represent the direction of time. (d) Two equilibria
$\mathbf{x^{*}_{\emptyset}}$ (blue dot) and $\mathbf{x}^{*}_{\\{86\\}}$ (red
dot) move along paths $\gamma_{\emptyset}(m)$ (blue curve) and
$\gamma_{\\{86\\}}(m)$ (red curve) with increasing $m$ until they meet in a
transcritical bifurcation (purple dot) at $m=m_{\\{86\\}}$. A solid(dashed)
line represents that the corresponding equilibrium is stable (unstable). (e)
As $m$ is further increased, $\mathbf{x}^{*}_{\\{86\\}}$ becomes stable and
$\mathbf{x^{*}_{\emptyset}}$ becomes unstable and unfeasible. At the same
time, another equilibrium $\mathbf{x}^{*}_{\\{86,76\\}}$ moves along a path
$\gamma_{\\{86,76\\}}(m)$ on the $x_{95}$ axis. (f) Another transcritical
bifurcation occurs at $m=m_{\\{86,76\\}}$ between $\mathbf{x}^{*}_{\\{86\\}}$
and $\mathbf{x}^{*}_{\\{86,76\\}}$. For $m>m_{\\{86,76\\}}$, the equilibrium
$\mathbf{x}^{*}_{\\{86,76\\}}$ is stable. (c’) The numerical counterpart of
(c) showing the solution trajectories from different initial conditions
converging to $\mathbf{x^{*}_{\emptyset}}$, which are obtained by numerically
integrating system (1) at $m=0$ and $c=0.3$ (f’) The numerical counterpart of
(f) showing the paths of the three equilibria undergoing two transcritical
bifurcations, obtained from the numerical continuation data at $c=0.3$.

Here we focus on the bifurcations generating the persistence diagram in Fig.
1(a3), which represents the survival and extinction of individual species when
the mutualism strength $m$ is increased slowly enough for the system to reach
equilibrium at each given $m$. We obtain the corresponding equilibria by
tracing the attracting feasible equilibria reached from the full-coexistence
equilibrium while increasing $m$ from zero within the numerical continuation
framework.

In Fig. 1(a3), at fixed $c=0.3$, the species labeled $86$ goes first extinct
at $m=m_{\\{86\\}}\approx 0.10471$ while all the $107$ species coexist for
$m<m_{\\{86\\}}$. As $m$ is increased further, the species $76$ next goes
extinct at $m_{\\{86,76\\}}\approx 0.10558$, and then two species $95$ and
$97$ go simultaneously extinct at $m_{\\{86,76,95,97\\}}\approx 0.10622$.

These extinction events are characterized by the exchanges of stability among
four equilibrium points,
$\mathbf{x_{\emptyset}}^{*},\mathbf{x}^{*}_{\\{86\\}},\mathbf{x}^{*}_{\\{86,76\\}}$,
and $\mathbf{x}^{*}_{\\{86,76,95,97\\}}$. Figure 2(a) shows the largest
eigenvalues, having the largest real parts, of the Jacobian at the equilibria
$\mathbf{x_{\emptyset}^{*}},\mathbf{x}_{\\{86\\}}$, and
$\mathbf{x}_{\\{86,76\\}}$. As $m$ increases passing $m_{\\{86\\}}$, the
largest eigenvalue of $\mathbf{x_{\emptyset}^{*}}$ becomes positive, while the
largest eigenvalues of $\mathbf{x}^{*}_{\\{86\\}}$ becomes negative, implying
that the former becomes unstable and the latter stable. Such exchange of
stability occurs also between $\mathbf{x}_{\\{86\\}}$ and
$\mathbf{x}_{\\{86,76\\}}^{*}$ at $m_{\\{86,76\\}}$, and between
$\mathbf{x}_{\\{86,76\\}}^{*}$ and $\mathbf{x}_{\\{86,76,95,97\\}}^{*}$ at
$m_{\\{86,76,95,97\\}}$. It should be noted that the two largest eigenvalues
of the Jacobian at $\mathbf{x}_{\\{86,76\\}}^{*}$ become positive as $m$
increases beyond $m_{\\{86,76,95,97\\}}$ [Fig. 2(b)], since they are equal due
to the symmetry of the interaction matrix $\mathbb{B}^{+}$ regarding two
species $95$ and $97$ that are connected to the same partner (plant) species
given the extinction of species $86$ and $76$.

All these transcritical bifurcations occur when the involved equilibria
collide in the phase space. In the phase space projected onto the abundance of
three species $(x_{95},x_{76},x_{86})$ for ease of representation, the full-
coexistence equilibrium $\mathbf{x}_{\emptyset}^{*}$, that is globally
attracting as represented by the inward flow in Fig. 2(c), moves towards the
$(x_{95},x_{76})$-plane, at which $x_{86}=0$, as $m$ increases until it
collides at $m=m_{\\{86\\}}$ with another equilibrium
$\mathbf{x}_{\\{86\\}}^{*}$ which is unstable as long as $m<m_{\\{86\\}}$
[Fig. 2(d)]. Their stability is exchanged at $m_{\\{86\\}}$ as described
above, indicating a transcritical bifurcation. For $m>m_{\\{86\\}}$, the
equilibrium $\mathbf{x}_{\\{86\\}}^{*}$ is stable moving in the
$(x_{95},x_{76})$-plane while the former equilibrium
$\mathbf{x}_{\emptyset}^{*}$ is unstable and unfeasible, moving in the space
of $x_{86}<0$ [Fig. 2(e)]. As $m$ is increased further,
$\mathbf{x}^{*}_{\\{86\\}}$ moves towards the $x_{95}$-axis until it coincides
with another equilibrium $\mathbf{x}^{*}_{\\{86,76\\}}$ at $m_{\\{86,76\\}}$,
which is another transcritical bifurcation. In Fig. 2(f), as $m$ further
increases, the equilibrium $\mathbf{x}^{*}_{\\{86,76\\}}$ moves towards the
equilibrium $\mathbf{x}^{*}_{\\{86,76,95,97\\}}$, located at the origin in
this projected phase space, and they collide and exchange stability at
$m_{\\{86,76,95,97\\}}$.

The described bifurcations are based on our numerical data showcased in Fig.
2(c’) and (f’). Particularly, Fig. 2(c’) illustrates the solution trajectories
that converge to $\mathbf{x^{*}_{\emptyset}}$ from different initial
conditions when $m=0$, the numerical counterpart of Fig. 2(c). Figure 2(f’)
shows the continuation paths of the three equilibria, the counterpart of the
ones schematically sketched in Fig. 2(f). Such transcritical bifurcations
continue until the attracting equilibrium escapes at infinity, as exemplified
by the termination of the persistence lines for the final 11 surviving species
at $m\approx m_{\infty}\approx 0.5109$ at $c=0.3$ in Fig. 1(a3).

###  IV.2 Multistability arising from Hopf bifurcations

Figure 3: Occurrence of Hopf bifurcation and multistability. (a) Three
equilibria in the $(x_{4},x_{9},x_{53})$-space for $m<m^{\rm(hopf)}$ with $c$
fixed at $0.3$. Schematic trajectories show that
$\mathbf{x}^{*}_{\mathcal{J}_{1}}$ (blue dot) is an attracting equilibrium
while $\mathbf{x}^{*}_{\mathcal{J}_{2}\cup\\{53\\}}$ (red dot) and
$\mathbf{x}^{*}_{\mathcal{J}_{2}}$ (red dot) are saddle.
$\mathbf{x}^{*}_{\mathcal{J}_{2}}$ is unfeasible. Also shown are the unstable
manifold $W^{u}(\mathbf{x}^{*}_{\mathcal{J}_{2}})$ (orange surface) and the
hyperplane $x_{53}=0$ (grey surface). Inset: Two largest eigenvalues
$\lambda$’s of the Jacobian at $\mathbf{x}^{*}_{\mathcal{J}_{2}}$ (black dots)
in the complex plane with the range $|\rm{re}(\lambda)|<0.8\cdot 10^{-4}$ and
$|\rm{imag}(\lambda)|<1.2\cdot 10^{-4}$. The purple and magenta curves
indicate how these eigenvalues will move as $m$ increases. (b) At
$m^{\rm(hopf)}$, a subcritical Hopf bifurcation occurs for
$\mathbf{x}^{*}_{\mathcal{J}_{2}}$. For $m\gtrsim m^{\rm(hopf)}$, an unstable
saddle periodic solution $\Gamma$ (green curve) is created, which has a
$(s-1)$-dimensional stable manifold $W^{s}(\Gamma)$ (cyan surface) and two-
dimensional unstable manifold $W^{u}(\Gamma)$ (red surface). Consequently the
equilibrium $\mathbf{x}^{*}_{\mathcal{J}_{2}}$ becomes stable (blue dot). The
equilibrium $\mathbf{x}^{*}_{\mathcal{J}_{2}\cup\\{53\\}}$ remains saddle.
Inset: The two largest eigenvalues of $\mathbf{x}^{*}_{\mathcal{J}_{2}}$ are
now complex conjugates with a negative real part. (c) Two equilibria
$\mathbf{x}^{*}_{\mathcal{J}_{2}\cup\\{53\\}}$ and
$\mathbf{x}^{*}_{\mathcal{J}_{2}}$ coincide (purple dot) when $m=m^{(2)}$,
which is a transcritical bifurcation. Inset:
$\mathbf{x}^{*}_{\mathcal{J}_{2}}$ has one zero eigenvalue. (d) For $m\gtrsim
m^{(2)}$, the equilibrium $\mathbf{x}^{*}_{\mathcal{J}_{2}}$ becomes a
feasible saddle equilibrium with a $(s-1)$-dimensional stable manifold
$W^{s}(\mathbf{x}^{*}_{\mathcal{J}_{2}})$ (blue surface) and one-dimensional
unstable manifold $W^{u}(\mathbf{x}^{*}_{\mathcal{J}_{2}})$ (orange curve). On
the other hand, $\mathbf{x}^{*}_{\mathcal{J}_{2}\cup\\{53\\}}$ becomes stable
(blue dot). Inset: $\mathbf{x}^{*}_{\mathcal{J}_{2}}$ has one positive
eigenvalue.

For large $m$, we find that multiple attracting equilibria may exist. One of
these corresponds to a continuation of increasing $m$ from the full-
coexistence equilibrium, while the remaining ones do not. Here we investigate
within the numerical continuation framework the underlying geometrical
mechanisms for the different equilibria as shown in Fig. 1(a3) and (b3).

We find a new feasible attracting equilibrium
$\mathbf{x}_{\mathcal{J}_{2}\cup\\{53\\}}^{*}$ shown in Fig. 1(b3) emerges at
$m=m^{\rm(2)}\approx 0.187267$, in addition to the one
$\mathbf{x}_{\mathcal{J}_{1}}^{*}$ of Fig. 1(a3), when $c$ is fixed at $0.3$.
The indices of extinct species in $\mathcal{J}_{1}$ and $\mathcal{J}_{2}$ are
quite different [Table 1 in Appendix B]. We find by numerical continuation
that $\mathbf{x}_{\mathcal{J}_{2}\cup\\{53\\}}^{*}$ acquires stability from an
unfeasible attracting equilibrium $\mathbf{x}_{\mathcal{J}_{2}}^{*}$, which
has a negative abundance for species $53$, through a transcritical bifurcation
at $m=m^{\rm(2)}$. The unfeasible equilibrium
$\mathbf{x}_{\mathcal{J}_{2}}^{*}$ is not stable for $m<m^{\rm(hopf)}\approx
0.187251$ but becomes stable by a Hopf bifurcation, which cannot occur for any
feasible equilibrium in system (1) with a symmetric interaction matrix
$\mathbb{B}$[44], at $m=m^{\rm(hopf)}$ smaller than $m^{(2)}$.

We sketch the phase portrait showing these two bifurcations in the projected
space of $x_{4},x_{9}$ and $x_{53}$ in Fig. 3. When $m<m^{\rm(hopf)}$, there
exists an attracting feasible equilibrium $\mathbf{x}^{*}_{\mathcal{J}_{1}}$.
In addition, there are two unstable saddle equilibria
$\mathbf{x}^{*}_{\mathcal{J}_{2}\cup\\{53\\}}$ and
$\mathbf{x}^{*}_{\mathcal{J}_{2}}$ at which the Jacobian has one and two
positive eigenvalues, respectively [Fig. 3(a) and the inset therein]. Here,
$\mathbf{x}^{*}_{\mathcal{J}_{2}\cup\\{53\\}}$ is a feasible equilibrium, and
$\mathbf{x}^{*}_{\mathcal{J}_{2}}$ is an unfeasible one with a negative
component for species $53$. The equilibrium $\mathbf{x}^{*}_{\mathcal{J}_{2}}$
possesses a two-dimensional unstable manifold
$W^{u}(\mathbf{x}^{*}_{\mathcal{J}_{2}})$ (orange surface) and a
$(s-2)$-dimensional stable manifold $W^{s}(\mathbf{x}^{*}_{\mathcal{J}_{2}})$
(not shown) which correspond to the set of initial conditions that converge
backward and forward in time to $\mathbf{x}^{*}_{\mathcal{J}_{2}}$,
respectively.

As $m$ increases past $m^{\rm(hopf)}$, a saddle periodic orbit $\Gamma$ (green
curve) surrounding $\mathbf{x}^{*}_{\mathcal{J}_{2}}$ (blue dot) is created,
which possesses a two-dimensional unstable manifold $W^{u}(\Gamma)$ (red
surface) and $(s-1)$-dimensional stable manifold $W^{s}(\Gamma)$ (cyan
surface), indicating a subcritical Hopf bifurcation. Therefore the unfeasible
equilibrium $\mathbf{x}^{*}_{\mathcal{J}_{2}}$ (blue dot) is now attracting
[Fig. 3(b)].

Then $\mathbf{x}^{*}_{\mathcal{J}_{2}}$ moves towards the
$(x_{4},x_{9})$-plane as $m$ increases further, and collides with the feasible
but currently unstable equilibrium
$\mathbf{x}^{*}_{\mathcal{J}_{2}\cup\\{53\\}}$ at $m=m^{\rm(2)}$ in a
transcritical bifurcation as shown in Fig. 3(c). Through this bifurcation, the
equilibrium $\mathbf{x}^{*}_{\mathcal{J}_{2}\cup\\{53\\}}$ becomes attracting,
and $\mathbf{x}^{*}_{\mathcal{J}_{2}}$ becomes a feasible saddle equilibrium
with one positive real eigenvalue, as shown in Fig. 3(d). Its $(s-1)$-stable
manifold $W^{s}(\mathbf{x}^{*}_{\mathcal{J}_{2}})$ becomes the boundary that
separates the initial conditions that converge to
$\mathbf{x}^{*}_{\mathcal{J}_{1}}$ and
$\mathbf{x}^{*}_{\mathcal{J}_{2}\cup\\{53\\}}$.

For the new feasible attracting equilibrium
$\mathbf{x}^{*}_{\mathcal{J}_{2}\cup\\{53\\}}$, transcritical bifurcations can
occur as $m$ increases representing the successive extinction of species as
described in Sec. IV.1 and shown in Fig. 1 (b3). It should be also noted that
more new equilibria can appear via the Hopf bifurcations of different
unfeasible equilibria; indeed, for the network $A_{1}$ we observe multiple
attracting equilibria arising through Hopf bifurcations as shown in Fig. 7 in
Appendix C (again for $c=0.3$).

###  IV.3 Multistability arising from Degenerate Transcritical bifurcation

Figure 4: Occurrence and consequences of a degenerate transcritical
bifurcation. (a1) Final abundances of individual species for any random
initial condition, corresponding to the only one attracting equilibrium, at
$m=0.28$ in networks $A_{2}$. (b1, b2, b3) Three examples among infinitely
many different equilibria observed at $m=0.3$. (c1, c2) Two attracting
equilibria at $m=0.31$. $c$ is fixed at $0.4$ in all panels. Figure 5:
Degenerate transcritical bifurcation and multistability. (a) Phase portrait in
the $(x_{4},x_{6})$-space for $m<m^{\rm(d)}$, where shown are an attracting
equilibrium $\mathbf{x}^{*}_{\mathcal{J}_{3}}$ (blue dot), two saddle
equilibria $\mathbf{x}^{*}_{\mathcal{J}_{3}\cup\\{4\\}}$ and
$\mathbf{x}^{*}_{\mathcal{J}_{3}\cup\\{6\\}}$ (red dots), each one having only
one positive eigenvalue, and $\mathbf{x}^{*}_{\mathcal{J}_{3}\cup\\{4,6\\}}$
(black dot) with two positive eigenvalues of the Jacobian. (b) At
$m=m^{\rm(d)}$, a degenerate transcritical bifurcation occurs to form a one-
dimensional line of nonhyperbolic equilibria (line with black dots). (c) For
$m\gtrsim m^{\rm(d)}$, $\mathbf{x}^{*}_{\mathcal{J}\cup\\{4\\}}$ and
$\mathbf{x}^{*}_{\mathcal{J}\cup\\{6\\}}$ are stable and
$\mathbf{x}^{*}_{\mathcal{J}}$ becomes an unstable saddle equilibrium. (d) At
$m=m_{\mathcal{J}_{3}\cup\\{4,6\\}}$, four equilibria coincide. (e) For
$m\gtrsim m_{\mathcal{J}_{3}\cup\\{4,6\\}}$, the equilibrium
$\mathbf{x}^{*}_{\mathcal{J}\cup\\{4,6\\}}$ becomes stable and the other three
equilibria become unfeasible. (a’, b’, c’) The numerical counterparts of (a),
(b) and (c), respectively, showing the solution trajectories of system (1)
with different initial conditions at $c=0.4$ and $m=0.28$ for (a’), $m=0.3$
for (b’), and $m=0.31$ for (c’).

Our study demonstrates that multistability can be generated also by degenerate
transcritical bifurcations for the LV dynamics in system (1). Degenerate
transcritical transitions occur when the reduced interaction matrix
$\mathbb{B^{+}}$ for surviving species has zero determinant, as described in
Sec. III. By considering the network $A_{2}$ as a representative example, we
showcase how this mechanism arises in a real network.

As $m$ increases from zero at fixed $c=0.4$ in network $A_{2}$, the full-
coexistence equilibrium $\mathbf{x_{\emptyset}^{*}}$ exhibits a series of
transcritical bifurcations until $m$ reaches $m^{\rm(d)}=0.3$. At $m=0.29$,
the final abundances of species are obtained as given in Fig. 4 (a1) without
regard to the initial conditions, implying the existence of only one
attracting equilibrium. At $m=m^{\rm(d)}=0.3$, in contrast, one finds the
final abundance of species $4$ and $6$ taking infinitely many different values
depending on the initial conditions, three of which are shown in Fig. 4(b1) to
(b3) as an example. When $m$ is slightly larger than $m^{\rm(d)}$, which we
will denote in this paper as $m\gtrsim m^{\rm(d)}$, we find only two final
abundances shown in Fig. 4(c), one with species $4$ surviving and the other
with species $6$ surviving; they cannot coexist simultaneously. The survival
and extinction of other species are not different between the two solutions.
We remark that during this transition, species $4$ and $6$ do not have the
same surviving mutualistic partners; such equivariance under the permutation
of their labels [54] is not necessary for this transition to occur.

The geometrical mechanism for this type of multistability, emerging through a
degenerate transcritical bifurcation, is illustrated in Fig. 5. At
$m<m^{\rm(d)}$, there exists an attracting feasible equilibrium
$\mathbf{x}^{*}_{\mathcal{J}_{3}}$ (blue dot), with $\mathcal{J}_{3}$ given in
Table 2 in Appendix B, and two saddle feasible equilibria
$\mathbf{x}^{*}_{\mathcal{J}_{3}\cup\\{4\\}}$ and
$\mathbf{x}^{*}_{\mathcal{J}_{3}\cup\\{6\\}}$ (red dots) [Fig. 5(a)]. At
$m=m^{\rm(d)}$ [Fig. 5(b)] a degenerate transcritical bifurcation occurs such
that $\mathbf{x}^{*}_{\mathcal{J}_{3}}$,
$\mathbf{x}^{*}_{\mathcal{J}_{3}\cup\\{4\\}}$ and
$\mathbf{x}^{*}_{\mathcal{J}_{3}\cup\\{6\\}}$ all become nonhyperbolic (having
an imaginary eigenvalue of the Jacobian). They are connected through a line of
nonhyperbolic equilibria (purple line) [Fig. 5(b)] corresponding to the one-
dimensional nullspace of the reduced interaction matrix $\mathbb{B^{+}}$ for
$\mathbf{x}^{*}_{\mathcal{J}}$. Different initial conditions converge to
different equilibria on the line. When $m$ is larger than $m^{\rm(d)}$, the
line of equilibria disappears, and the stability of equilibria is
interchanged; $\mathbf{x}^{*}_{\mathcal{J}_{3}}$ becomes saddle, and
$\mathbf{x}^{*}_{\mathcal{J}_{3}\cup\\{4\\}}$ and
$\mathbf{x}^{*}_{\mathcal{J}_{3}\cup\\{6\\}}$ become attracting [Fig. 5(c)].
Notably, $\mathbf{x}^{*}_{\mathcal{J}_{3}}$ has a stable manifold (blue line)
that separates the basin of attraction of the two new attracting equilibria;
thus, certain initial conditions can converge to
$\mathbf{x}^{*}_{\mathcal{J}_{3}\cup\\{4\\}}$ while others to
$\mathbf{x}^{*}_{\mathcal{J}_{3}\cup\\{6\\}}$. Therefore, we see two new
attracting equilibria where a species is extinct in one and surviving in the
other. Figure 5(a’), (b’) and (c’) showcase with actual numerical simulations
how the solutions converge to different equilibria during this transition as
indicated in their schematic sketches counterparts in Fig. 5(a), (b), and (c).

As $m$ is further increased after the degenerate transcritical bifurcation,
the two attracting equilibria can exhibit their respective extinction
sequences, keeping multistability or multistability may disappear soon, the
latter of which is the case in our example. At
$m=m_{\mathcal{J}_{3}\cup\\{4,6\\}}$, as shown in Fig. 5(d), all the
equilibria $\mathbf{x}^{*}_{\mathcal{J}_{3}\cup\\{4\\}}$,
$\mathbf{x}^{*}_{\mathcal{J}_{3}\cup\\{6\\}}$ and
$\mathbf{x}^{*}_{\mathcal{J}_{3}}$ coincide with the equilibrium
$\mathbf{x}^{*}_{\mathcal{J}_{3}\cup\\{4,6\\}}$. When $m$ is further
increased, $\mathbf{x}^{*}_{\mathcal{J}_{3}\cup\\{4,6\\}}$ is the only
attracting equilibrium [Fig. 5(e)]. Given that the persistence diagrams of
these two equilibria are equivalent up to $m^{\rm(d)}$ where this bifurcation
occurs (blue line) [Fig. 8 in Appendix D], it is only for
$m^{\rm(d)}<m<m_{\mathcal{J}_{3}\cup\\{4,6\\}}$ that both diagrams differ,
meaning multistability. For $m>m_{\mathcal{J}_{3}\cup\\{4,6\\}}$, both species
4 and 6 are extinct and multistability disappears.

It is worth mentioning that multistability from Hopf bifurcations can still
occur in network $A_{2}$ as $m$ increases. Also, we should remark that the
equilibria connected to the full-coexistence equilibrium in network $A_{1}$
undergo a degenerate transcritical bifurcation though it is not shown in Fig.
1(a3).

##  V Phase diagram

Figure 6: Bifurcation diagram for network $A_{1}$ (a), and $A_{2}$ (b). Shown
are representative lines of transcritical bifurcations $\mathbf{T}$ (brown and
orange), Hopf bifurcation $\mathbf{H}$ (green), degenerate transcritical
bifurcation $\mathbf{DT}$ (blue and cyan) and the boundary of bounded
asymptotic behavior $\mathbf{U}$ (black). The line $\mathbf{B}$ (purple)
represents the smallest values of $m$ and $c$, where the full interaction
matrix $\mathbb{B}$ in Eq. (2) stops being negative definitive. The subindex
indicates the extinct species of the equilibrium exhibiting the bifurcation.
See Fig. 7 in Appendix C for persistence diagram of
$\mathbf{x}^{*}_{\mathcal{J}_{4}}$ and $\mathbf{x}^{*}_{\mathcal{J}_{5}}$.

The numerical continuation approach allows one to obtain the loci of
bifurcations in the parameter space and thereby identify the parameter regimes
or phases characterized by distinct sets of extinct species and
multistability. Such a phase diagram is presented in Fig. 6(a) for network
$A_{1}$, in which the first transcritical bifurcation line labeled as
$\mathbf{T}_{\\{86\\}}$ (light brown) is drawn, distinguishing the full-
coexistence phase $\mathbf{x^{*}_{\emptyset}}$ and the phase with species 86
extinct $\mathbf{x}_{\\{86\\}}^{*}$. Above the line $\mathbf{T}_{\\{86\\}}$ in
the $(c,m)$-plane, a series of transcritical bifurcation lines are located,
representing further extinctions, some of which we present by darker brown
lines. We also illustrate representative Hopf bifurcations, represented by
green lines, which correspond to the mechanism of creating multistability;
indeed, really close to these lines we find transcritical bifurcations (not
shown) that mark the onset of multistability in system (1), above which
multiple attracting feasible equilibria coexist.

As stated earlier, the full coexistence equilibrium
$\mathbf{x}^{*}_{\emptyset}$ is globally attracting if all of its components
are positive and the full interaction matrix $\mathbb{B}$ is negative
definite. To understand when this global property is lost, we obtain the
boundary where $\mathbb{B}$ stops being negative definitive, which is the line
labeled as $\mathbf{B}$ in Fig. 6(a). Below this line, $\mathbb{B}$ is
negative definite. As it turns out that the line $\mathbf{B}$ is located above
the first transcritical bifurcation line $\mathbf{T}_{\\{86\\}}$, the full-
coexistence equilibrium $\mathbf{x}_{\emptyset^{*}}$ is a global attractor and
there is no multistability in the orange colored region in Fig. 6(a).

We find that all these bifurcation lines (phase boundaries) emanate from the
point $(c,m)=(1,0)$. System (1) is highly degenerate at $(c,m)=(1,0)$ as it is
reduced to ${dx_{i}\over
dt}=x_{i}(1-\sum_{p^{\prime}=1}^{n_{p}}x_{p^{\prime}})$ for all plants $i$,
and ${dx_{j}\over dt}=x_{j}(1-\sum_{a^{\prime}=1}^{n_{a}}x_{a^{\prime}})$ for
all animals $j$, leading to a family of nonhyperbolic equilibria in the
$(s-2)$-dimensional hyperplane given by the intersection of the hyperplanes
$\sum_{p}x_{p}=1$ and $\sum_{a}x_{a}=1$.

We also compute the boundary for unbounded abundances $\mathbf{U}$, where the
equilibria disappear at infinity. In Fig. 6(a), it is shown that the
transcritical bifurcation lines, which correspond to the successive
extinctions of the full coexistence equilibrium, terminate at the black line
$\mathbf{U}$. In the persistence diagrams in Figs. 1(a3), 1(b3), 7(a), and
7(c), the last surviving species find their abundances diverging at certain
values of $m$. To obtain this boundary numerically, we first need to identify
the set of the last surviving species that cannot undergo further
transcritical bifurcations. To do so, we fix $m$ at a value before the
equilibrium diverges, and then we apply continuation on $c$ until it is equal
to one. Then, we fix $c$ at one and then continue on $m$. During this
procedure, we jump to new attracting equilibria if further transcritical
bifurcation occurs. In this way, we can find the last set of surviving species
when $c$ is one. After identifying this set, we can restart the continuation
procedure by varying $m$ with $c=0$. Since only mutualism occurs,
$\mathbf{x^{*}_{\emptyset}}$ cannot exhibit extinction, and it will grow
unbounded as $m$ increases. We then choose a value of $m^{\rm(u)}$ for which
one of the abundances of the set of the last surviving species becomes large,
e.g., $x_{\rm surv}=10^{4}$. We then obtain the equilibrium by solving the
zero-problem in Eq. (3) under the constraint $x_{\rm surv}=10^{4}$, and trace
by numerical continuation this equilibrium and its transcritical bifurcations
while simultaneously varying $c$ and $m$. This leads us to approximate the
terminus of the transcritical bifurcation lines numerically. We remark that
for each attracting equilibrium that emerges through the Hopf bifurcation
mechanism, there exists an associated boundary of unbounded abundance (not
shown). This demonstrates that compactification techniques can be used to
study these bifurcations [55, 56, 51, 52, 53].

The phase diagram for network $A_{1}$ in Fig. 6 (a) implies much. As the
transcritical bifurcation lines, emanating from the degenerate point $(1,0)$,
have negative slopes in the $(c,m)$-plane, one can expect that an increase in
mutualism or competition can lead to species extinction; If the competition
(mutualism) strength is increased, then the mutualism (competition) strength
should be decreased sufficiently to avoid species extinction. It also suggests
that as competition in the network intensifies, multistability can arise even
with weak mutualism. This indicates that biodiversity scenarios in highly
competitive uniform networks are quite fragile, even when weak nonuniform
mutualistic dynamics are introduced.

The phase diagram for network $A_{2}$ is also shown in Fig. 6(b), in which we
present i) the first transcritical bifurcation line (brown) for the extinction
of species 109, ii) the degenerate transcritical bifurcation line (blue), and
iii) the transcritical bifurcation (orange) that distinguish the regime of
either species 4 or 6 surviving and that of both extinct. Like for network
$A_{1}$, all these lines emanate from the degenerate point $(c,m)=(1,0)$,
implying its role as the organizing center for the transcritical, Hopf, and
degenerate transcritical bifurcation lines that create extinction and
multistability.

Interestingly, multistability associated with the survival and extinction of
species $4$ and $6$ occurs only when the competition strength is large enough
for the degenerate transcritical bifurcation to occur at a smaller value of
$m$ than that for the transcritical bifurcation. For instance, at a small
value of competition strength, $c=0.250$, no multistability arises; the phase
portrait transits from Fig. 5(a) to (e) via (d), without (b) or (c), at
$m\approx 0.3645$. That is, species $4$ and $6$ become extinct simultaneously,
and there is no range of $m$ where one of them exclusively exists. We also
present in Fig. 6(b) the locus where the equilibria disappear at infinity
(black line).

##  VI Summary and Discussion

In this work, we considered a model community of two species groups— plants
and pollinators—under uniform intragroup competition and empirical
heterogeneous intergroup mutualism. Through a geometrical and numerical
continuation approach, we provided an exact numerical method to systematically
pinpoint the loci of the extinction of individual species in the parameter
space for two empirical networks. Our approach has pushed us to understand
emergent phenomena in structured communities, whereas most theoretical and
numerical approaches have so far studied the random, unstructured case.

Most importantly, through our approach, we were able to identify the specific
bifurcations responsible for the creation of multistability in the empirical
networks we studied. These include subcritical Hopf bifurcations and
degenerate transcritical bifurcations, which serve as mechanisms for creating
multiple attracting equilibria. The multistability that they generate exhibit
different characteristics. Specifically, the Hopf bifurcations lead to the
formation of new attractive equilibria where the composition of surviving
species may not overlap with that of pre-existing equilibria. Conversely,
degenerate transcritical bifurcations induce multistability by generating new
attracting equilibria that differ only in the finite set of species that
survive.

Given that the present study is limited to selected large empirical networks
and simple model dynamics, it is desirable to extend our findings to a broader
range of interacting communities. This extension should include a thorough
investigation of the structural and dynamical conditions for the emergence of
multistability, specifically through mechanisms such as Hopf and degenerate
transcritical bifurcations. A deeper theoretical exploration of how empirical
structural characteristics of ecological communities, such as degree
heterogeneity, nestedness, non-uniform interaction strengths, and higher-order
interactions, affect bifurcations and multistability is important.
Furthermore, a comprehensive analysis of the stability of all equilibrium
points in smaller systems, utilizing numerical continuation methods developed
in this study, will enhance our understanding of these complex dynamics.

###### Acknowledgements.

This work was supported by KIAS Individual Grants [Nos CG086101 (A.G.) and
CG079902 (D.-S.L.)] at Korea Institute for Advanced Study.

## References

  * May [1972] R. M. May, Will a Large Complex System be Stable?, Nature 238, 413 (1972).
  * Roberts [1974] A. Roberts, The stability of a feasible random ecosystem, Nature 251, 607 (1974).
  * Takeuchi [1996] Y. Takeuchi,  _Global Dynamical Properties of Lotka-Volterra Systems_ (WORLD SCIENTIFIC, 1996).
  * Rieger [1989] H. Rieger, Solvable model of a complex ecosystem with randomly interacting species, J. Phys. A 22, 3447 (1989).
  * Allesina and Tang [2012] S. Allesina and S. Tang, Stability criteria for complex ecosystems, Nature 483, 205 (2012).
  * Allesina and Tang [2015] S. Allesina and S. Tang, The stability–complexity relationship at age 40: a random matrix perspective, Population Ecology 57, 63 (2015).
  * Bunin [2017] G. Bunin, Ecological communities with Lotka-Volterra dynamics, Phys. Rev. E 95, 1 (2017).
  * Garcia Lorenzana and Altieri [2022] G. Garcia Lorenzana and A. Altieri, Well-mixed Lotka-Volterra model with random strongly competitive interactions, Phys. Rev. E 105, 1 (2022).
  * Gibbs  _et al._ [2018] T. Gibbs, J. Grilli, T. Rogers, and S. Allesina, Effect of population abundances on the stability of large random ecosystems, Phys. Rev. E 98, 1 (2018).
  * Mambuca  _et al._ [2022] A. M. Mambuca, C. Cammarota, and I. Neri, Dynamical systems on large networks with predator-prey interactions are stable and exhibit oscillations, Phys. Rev. E 105, 14305 (2022).
  * Bascompte  _et al._ [2003] J. Bascompte, P. Jordano, C. J. Melián, and J. M. Olesen, The nested assembly of plant–animal mutualistic networks, Proc. Natl. Acad. Sci. U.S.A. 100, 9383 (2003).
  * Thébault and Fontaine [2010] E. Thébault and C. Fontaine, Stability of ecological communities and the architecture of mutualistic and trophic networks, Science 329, 853 (2010).
  * Suweis  _et al._ [2013] S. Suweis, F. Simini, J. R. Banavar, and A. Maritan, Emergence of structural and dynamical properties of ecological mutualistic networks, Nature 500, 449 (2013).
  * Maeng  _et al._ [2012] S. E. Maeng, J. W. Lee, and D.-S. Lee, Interspecific competition underlying mutualistic networks, Phys. Rev. Lett. 108, 108701 (2012).
  * Stone [2016] L. Stone, The google matrix controls the stability of structured ecological and biological networks, Nat. Commun. 7, 12857 (2016).
  * Park  _et al._ [2024] J. I. Park, D.-S. Lee, S. H. Lee, and H. J. Park, Incorporating heterogeneous interactions for ecological biodiversity (2024), arXiv:2403.15730 [physics.bio-ph] .
  * Poley  _et al._ [2024] L. Poley, T. Galla, and J. W. Baron, Interaction networks in persistent lotka-volterra communities (2024), arXiv:2404.08600 [q-bio.PE] .
  * Fowler and MacMahon [1982] C. W. Fowler and J. A. MacMahon, Selective extinction and speciation: Their influence on the structure and functioning of communities and ecosystems, Am. Nat. 119, 480 (1982).
  * Assaf and Meerson [2006] M. Assaf and B. Meerson, Spectral theory of metastability and extinction in birth-death systems, Phys. Rev. Lett. 97, 200602 (2006).
  * Volkov  _et al._ [2003] I. Volkov, J. R. Banavar, S. P. Hubbell, and A. Maritan, Neutral theory and relative species abundance in ecology, Nature 424, 1035 (2003).
  * Suweis  _et al._ [2023] S. Suweis, F. Ferraro, S. Azaele, and A. Maritan, Generalized lotka-volterra systems with time correlated stochastic interactions (2023), arXiv:2307.02851 [q-bio.PE] .
  * de Oca and Zeeman [1996] F. de Oca and M. Zeeman, Extinction in nonautonomous competitive Lotka-Volterra systems, Proc. Am. Math. Soc. 124, 3677 (1996).
  * Baigent [2017] S. Baigent, Lotka–Volterra Dynamical Systems, in  _Dynamical Complex Systems_ (WORLD SCIENTIFIC (EUROPE), 2017) pp. 157–188.
  * Goh [1979] B. S. Goh, Stability in Models of Mutualism, Am. Nat. 113, 261 (1979).
  * Smale [1976] S. Smale, On the differential equations of species in competition, J. Math. Biol. 3, 5 (1976).
  * Zeeman [1995] M. L. Zeeman, Extinction in Competitive Lotka-Volterra Systems, Proc. Am. Math. Soc. 123, 87 (1995).
  * Zeeman and Zeeman [2002a] E. C. Zeeman and M. L. Zeeman, An n-dimensional competitive Lotka-Volterra system is generically determined by the edges of its carrying simplex, Nonlinearity 15, 2019 (2002a).
  * Zeeman and Zeeman [2002b] E. Zeeman and M. Zeeman, From local to global behavior in competitive Lotka-Volterra systems, Trans. Am. Math. Soc. 355, 713 (2002b).
  * [29] https://www.web-of-life.es/.
  * Montoya and Solé [2002] J. M. Montoya and R. V. Solé, Small world patterns in food webs, J. Theor. Biol. 214, 405 (2002).
  * Dunne  _et al._ [2002] J. A. Dunne, R. J. Williams, and N. D. Martinez, Food-web structure and network theory: the role of connectance and size, Proc. Natl. Acad. Sci. U.S.A. 99, 12917 (2002).
  * Pascual and Dunne [2006] M. Pascual and J. A. Dunne,  _Ecological networks: linking structure to dynamics in food webs_ (Oxford University Press, 2006).
  * Bastolla  _et al._ [2009] U. Bastolla, M. A. Fortuna, A. Pascual-García, A. Ferrera, B. Luque, and J. Bascompte, The architecture of mutualistic networks minimizes competition and increases biodiversity, Nature 458, 1018 (2009).
  * Rohr  _et al._ [2014] R. P. Rohr, S. Saavedra, and J. Bascompte, On the structural stability of mutualistic systems, Science 345, 1253497 (2014).
  * Maeng  _et al._ [2019] S. E. Maeng, J. W. Lee, and D.-S. Lee, Competition-induced increase of species abundance in mutualistic networks, J. Stat. Mech.: Theory Exp 2019, 033502 (2019).
  * Lee  _et al._ [2022] H. W. Lee, J. W. Lee, and D.-S. Lee, Stability and selective extinction in complex mutualistic networks, Phys. Rev. E 105, 014309 (2022).
  * Marcus  _et al._ [2024] S. Marcus, A. M. Turner, and G. Bunin, Extinctions as a vestige of instability: the geometry of stability and feasibility (2024), arXiv:2405.11360 [q-bio.PE] .
  * Arancibia-Ibarra  _et al._ [2021] C. Arancibia-Ibarra, P. Aguirre, J. Flores, and P. van Heijster, Bifurcation analysis of a predator-prey model with predator intraspecific interactions and ratio-dependent functional response, Appl. Math. Comput. 402, 126152 (2021).
  * Aguadé-Gorgorió  _et al._ [2024] G. Aguadé-Gorgorió, J.-F. Arnoldi, M. Barbier, and S. Kéfi, A taxonomy of multiple stable states in complex ecological communities, Ecol. Lett. 27, e14413 (2024).
  * Strogatz [2018] S. H. Strogatz,  _Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering_ (CRC press, 2018).
  * Kuznetsov [2004] Y. A. Kuznetsov,  _Elements of applied bifurcation theory_ , 3rd ed., Applied Mathematical Sciences, Vol. 112 (Springer-Verlag, New York, 2004).
  * Doedel [2007] E. J. Doedel, Lectures Notes on Numerical Analysis of Nonlinear Equations, in  _Numerical Continuation Methods for Dynamical Systems: Path Following and Boundary Value Problems_ , edited by B. Krauskopf, H. M. Osinga, and J. Galán-Vioque (Springer, The Netherlands, 2007) pp. 1–50.
  * Pascual-García and Bastolla [2017] A. Pascual-García and U. Bastolla, Mutualism supports biodiversity when the direct competition is weak, Nat. Commun. 8, 14326 (2017).
  * MacArthur [1970] R. MacArthur, Species Packing and Competitive, Theor. Popul. Biol. 11, 1 (1970).
  * Ros  _et al._ [2023] V. Ros, F. Roy, G. Biroli, G. Bunin, and A. M. Turner, Generalized lotka-volterra equations with random, nonreciprocal interactions: The typical number of equilibria, Phys. Rev. Lett. 130, 257401 (2023).
  * Lee  _et al._ [2012] D.-S. Lee, S. E. Maeng, and J. W. Lee, Scaling of nestedness in complex networks, J. Korean Phys. Soc. 60, 648 (2012).
  * Payrató-Borràs  _et al._ [2019] C. Payrató-Borràs, L. Hernández, and Y. Moreno, Breaking the spell of nestedness: The entropic origin of nestedness in mutualistic systems, Phys. Rev. X 9, 031024 (2019).
  * Doedel [1981] E. J. Doedel, AUTO: A program for the automatic bifurcation analysis of autonomous systems, Congr. Numer. 30, 265 (1981).
  * Doedel and Oldeman [2010] E. J. Doedel and B. E. Oldeman,  _AUTO- 07p: Continuation and Bifurcation Software for Ordinary Differential Equations_ (Department of Computer Science, Concordia University, Montreal, Canada, 2010).
  * Note [1] If more than one components, say the components in a set $\mathcal{K}$, become zero at $m_{*}$, then the equilibrium $\mathbf{x^{*}_{\mathcal{J}}}$ coincides with all the equilibria $\mathbf{x_{\mathcal{J}^{\prime}}^{*}}$ with all $\mathcal{J}^{\prime}$ in $\mathcal{P}(\mathcal{K})$, where $\mathcal{P}(\mathcal{K})$ is the power set of $\mathcal{K}$. This occurs if Eq. (1) is equivariant [54] under a permutation of species, i.e., when a collection of plants or animals have the same set of mutualistic partners.
  * Giraldo  _et al._ [2017] A. Giraldo, B. Krauskopf, and H. M. Osinga, Saddle invariant objects and their global manifolds in a neighborhood of a homoclinic flip bifurcation of case B, SIAM J. Appl. Dyn. Syst. 16, 640 (2017).
  * Giraldo  _et al._ [2020] A. Giraldo, B. Krauskopf, and H. M. Osinga, Computing connecting orbits to infinity associated with a homoclinic flip bifurcation, J. Comput. Dyn. 7, 489 (2020).
  * Matsue [2017] K. Matsue, On Blow-Up Solutions of Differential Equations with Poincaré-Type Compactifications, SIAM J. Appl. Dyn. Syst. 17 (2017).
  * Golubitsky and Schaeffer [1985] M. Golubitsky and D. G. Schaeffer,  _Singularities and Groups in Bifurcation Theory_ , Vol. 51 (Springer New York, 1985).
  * Messias [2009] M. Messias, Dynamics at infinity and the existence of singularly degenerate heteroclinic cycles in the Lorenz system, J. Phys. A 42 (2009).
  * Messias [2011] M. Messias, Dynamics at infinity of a cubic Chua’s system, Internat. J. Bifur. Chaos Appl. Sci. Engrg. 21, 333 (2011).

##  Appendix A Global attractiveness of the full coexistence equilibrium in
case of no mutualism

######  Theorem A.1 (Goh’s theorem [24]).

Suppose that the Lotka-Volterra system
$\dot{\mathbf{x}}=\mathbb{X}\left(\boldsymbol{r}+\mathbb{C}\mathbf{x}\right)$,
for $\boldsymbol{r}\in\mathbb{R}^{s}$ and $\mathbb{C}\in M_{s}(\mathbb{R})$,
has a unique interior equilibrium
$x^{*}=-\mathbb{C}^{-1}\boldsymbol{r}\in\mathbb{R}^{s}_{>0}$. Then this
equilibrium is globally attracting on $\mathbb{R}^{s}_{>0}$ if there exists a
diagonal matrix $\mathbb{D}>0$ (every entry is positive) such that
$\mathbb{C}\mathbb{D}+\mathbb{D}\mathbb{C}^{T}$ is negative definite.

Thus, we have the following:

######  Corollary A.2.

The matrix $\mathbb{B}_{0}:=-(1-c)\mathbb{I}-c\mathbb{J}$ is negative
definite. Thus, the equilibrium $\mathbf{x^{*}_{\emptyset}}$ is globally
attracting equilibrium in $\mathbb{R}^{s}_{>0}$ for system (1) when $m=0$ and
$0\leq c<1$.

###### Proof.

For an arbitrary vector $\mathbf{x}=(x_{1},x_{2},...,x_{s})\in\mathbb{R}^{s}$,
it follows from

$\mathbb{B}_{0}=\begin{pmatrix}-(1-c)\mathbf{I}_{n_{p}}-c\mathbf{J}_{n_{p}}&\mathbf{0}_{n_{p},n_{a}}\\\\[2.84526pt]
\mathbf{0}_{n_{a},n_{p}}&-(1-c)\mathbf{I}_{n_{a}}-c\mathbf{J}_{n_{a}}\end{pmatrix}$
(6)

that $\mathbf{x}^{T}\mathbb{B}_{0}\mathbf{x}<0$ (negative definiteness) is
equivalent to

$\mathbf{x}_{p}^{T}\left[(1-c)\mathbf{I}_{n_{p}}+c\mathbf{J}_{n_{p}}\right]\mathbf{x}_{p}+\mathbf{x}_{a}^{T}\left[(1-c)\mathbf{I}_{n_{a}}+c\mathbf{J}_{n_{a}}\right]\mathbf{x}_{a}>0,$

where $\mathbf{x}_{p}=(x_{1},x_{2},...,x_{n_{p}})$ and
$\mathbf{x}_{a}=(x_{n_{p}+1},x_{n_{p}+2},...,x_{n_{a}})$. Notice that

$\displaystyle\mathbf{x}_{p}^{T}\left[(1-c)\mathbf{I}_{n_{p}}+c\mathbf{J}_{n_{p}}\right]\mathbf{x}_{p}$
$\displaystyle=(1-c)\mathbf{x}_{p}^{T}\mathbf{x}_{p}+c\mathbf{x}_{p}^{T}\mathbf{J}_{n_{p}}\mathbf{x}_{p}$
$\displaystyle=(1-c)\sum^{n_{p}}_{i=1}x^{2}_{i}+c\sum^{n_{p}}_{i=1}\sum^{n_{p}}_{j=1}x_{i}x_{j}=(1-c)\sum^{n_{p}}_{i=1}x^{2}_{i}+c\left(\sum^{n_{p}}_{i=1}x_{i}\right)^{2}.$
(7)

It follows that Eq. (7) is strictly positive since $0\leq c<1$. By a
symmetrical argument, it is clear that
$\mathbf{x}_{a}^{T}\left[(1-c)\mathbf{I}_{n_{a}}+c\mathbf{J}_{n_{a}}\right]\mathbf{x}_{a}$
is also strictly positive, thus we have that

$\mathbf{x}_{p}^{T}\left[(1-c)\mathbf{I}_{n_{p}}+c\mathbf{J}_{n_{p}}\right]\mathbf{x}_{p}+\mathbf{x}_{a}^{T}\left[(1-c)\mathbf{I}_{n_{a}}+c\mathbf{J}_{n_{a}}\right]\mathbf{x}_{a}>0$

holds for any arbitrary vector $\mathbf{x}$; hence, $\mathbb{B}_{0}$ is
negative definite. Now by making $\mathbb{D}=\mathbb{I}$ and
$\mathbb{C}=\mathbb{B}_{0}$ then, by Goh’s theorem,
$\mathbf{x^{*}_{\emptyset}}$ is globally attracting equilibrium in
$\mathbb{R}^{s}_{>0}$ for system (1) when $m=0$ and $0\leq c<1$. ∎

##  Appendix B Sets of extinction species $\mathcal{J}_{1},\mathcal{J}_{2}$,
and $\mathcal{J}_{3}$

In this section, we provide in Table 1 and 2 the list of extinct species
$\mathcal{J}_{1},\mathcal{J}_{2}$, and $\mathcal{J}_{3}$ for the equilibria
studied in the main text.

| 6, 13, 14, 17, 18, 19, 21, 24, 25, 27, 28, 29, 30, 31, 32, 34, 35, 36, 37,
38, 39, 40, 41, 42, 43, 50, 51, 52, 54,
---|---
$\mathcal{J}_{1}$ | 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 70, 72, 73, 74, 75, 76, 77, 80, 82, 84, 85, 86, 87, 88, 89, 90,
| 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 102, 104, 105, 106, 107
| 13, 14, 18, 19, 22, 24, 25, 26, 27, 28, 30, 32, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 51, 56, 57, 60, 62, 67, 68, 69,
$\mathcal{J}_{2}$ | 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 88, 89, 90, 92, 93, 94, 95, 96, 97, 99, 100, 101,
| 102, 103, 104, 105, 106, 107
Table 1: Sets of extinct species $\mathcal{J}_{1}$ and $\mathcal{J}_{2}$ of the equilibria shown in Fig. 3.  | 3, 5, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 2931, 34, 39, 41, 42, 43, 46, 47, 48, 49,
---|---
$\mathcal{J}_{3}$ | 51, 53, 54, 55, 56, 57, 60, 61, 62, 63, 64, 65, 66, 67, 69, 70, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84,
| 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97,98, 99, 100, 101, 102,
103, 104, 105, 106, 107, 108, 109, 110
Table 2: Sets of extinct species $\mathcal{J}_{3}$ of the equilibrium
exhibiting the degenerate transcritical bifurcation for network $A_{2}$ at
$c=0.4$ and $m=0.3$ shown in Fig. 5

##  Appendix C Additional equilibria arising through Hopf bifurcations in
network $A_{1}$

In this section, we present the persistence diagram of two additional
attracting equilibria that emerge from Hopf bifurcations in network $A_{1}$.
Both families of equilibria shown in Fig. 7 exist for a narrow range of $m$
and disappear at infinity as $m$ increases.

Figure 7:  Additional multistability driven by Hopf bifurcations in network
$A_{1}$. Panels (a) and (c) show the persistence diagrams computed starting
from unfeasible equilibria, $\mathbf{x}^{*}_{\mathcal{J}_{4}}$ and
$\mathbf{x}^{*}_{\mathcal{J}_{5}}$ appearing in Fig. 6, respectively, that
become stable through Hopf bifurcations. The purple squares indicate the
species whose abundances are negative until transcritical bifurcations occur
turning feasible equilibria stable. The vertical dashed lines in panels (a)
and (c) indicate the values of $m$ selected to create panels (b1), (b2),(d1),
and (d2). Panels (b1) and (d1) show the time trajectories of the abundances
$x_{i}$’s for the initial conditions converging to the indicated equilibria in
panels (a) and (c), respectively, while panels (b2) and (d2) illustrate the
corresponding final abundances of species.

##  Appendix D Persistence diagrams displaying degenerate transcritical
bifurcations in network $A_{2}$

In this section, we present the persistence diagrams displaying a degenerate
transcritical bifurcation in network $A_{2}$.

Figure 8: Two different persistence diagrams, (a) and (b), at $c=0.4$ in
network $A_{2}$. They are computed starting from the full-coexistence
equilibrium $\mathbf{x_{\emptyset}^{*}}$ at $m=0$. The vertical lines are
drawn at $m=m^{\rm(d)}$ (blue) and $m=m_{\mathcal{J}_{3}\cup\\{4,6\\}}$
(orange), representing the degenerate transcritical bifurcation that creates
multistability, and the transcritical bifurcation that terminates it as
described in Sec. IV.3.