# Stability and selective extinction in complex mutualistic networks

Hyun Woo Lee  Jae Woo Lee<EMAIL_ADDRESS>Department of Physics, Inha
University, Incheon 22212, Korea  Deok-Sun Lee<EMAIL_ADDRESS>School
of Computational Sciences and Center for AI and Natural Sciences, Korea
Institute for Advanced Study, Seoul 02455, Korea

###### Abstract

We study species abundance in the empirical plant-pollinator mutualistic
networks exhibiting broad degree distributions, with uniform intra-group
competition assumed, by the Lotka-Volterra equation. The stability of a fixed
point is found to be identified by the signs of its non-zero components and
those of its neighboring fixed points. Taking the annealed approximation, we
derive the non-zero components to be formulated in terms of degrees and the
rescaled interaction strengths, which lead us to find different stable fixed
points depending on parameters, and we obtain the phase diagram. The selective
extinction phase finds small-degree species extinct and effective interaction
reduced, maintaining stability and hindering the onset of instability. The
non-zero minimum species abundances from different empirical networks show
data collapse when rescaled as predicted theoretically.

##  I Introduction

A community of randomly interacting species can become unstable as the number
of species and their interaction connectivity together go beyond a threshold
May (1972). Such a random-interaction model can be informative with the help
of the random matrix theory, and it has been instrumental in the theoretical
study of ecological systems, illuminating their features from the perspective
of stability Goh (1979); Allesina and Tang (2015); Stone (2016); Grilli  _et
al._ (2017); Bunin (2017); Cui  _et al._ (2020); Pettersson  _et al._ (2020a);
*Pettersson2020. Recently available data-sets point out the complex
organization of interspecific interactions, neither completely random nor
ordered Montoya and Solé (2002); Dunne  _et al._ (2002); Bascompte  _et al._
(2003); Jordano  _et al._ (2003); Montoya  _et al._ (2006); Bascompte and
Jordano (2007); Guimarães  _et al._ (2007); Olesen  _et al._ (2007); Thébault
and Fontaine (2010); Maeng and Lee (2011), and they have drawn the attention
of researchers to the origins and implications of over-represented network
structural features Bastolla  _et al._ (2009); Suweis  _et al._ (2013);
Saavedra  _et al._ (2016); Yan  _et al._ (2017).

Contrary to the unstructured communities in which every species is subject to
the identical randomness in its interaction profile, individual species can be
in fundamentally different states under structured interactions. For instance,
the mutualistic partnership between flowering plants and pollinating bees is
characterized by different numbers of partners, called degrees, from species
to species. Considering also the intrinsic competitions among plants and among
pollinators due to limited resources Maeng  _et al._ (2012); Suweis  _et al._
(2013); Pascual-García and Bastolla (2017); Barbier  _et al._ (2018); Gracia-
Lázaro  _et al._ (2018); Wang  _et al._ (2021); Maeng  _et al._ (2019); Cai
_et al._ (2020), one finds that the abundance of a species increase due to the
benefit from mutualistic partners but also decrease due to the cost from
competition, the imbalance of which may lead some species to flourish but
others to become extinct Gracia-Lázaro  _et al._ (2018); Wang  _et al._
(2021). The mechanism driving such different fates across species remains to
be elucidated James  _et al._ (2012); Saavedra and Stouffer (2013); Allesina
and Tang (2012); Stone (2020).

Here we investigate the different abundances and different likelihood of
extinction of individual species in heterogeneous mutualistic networks from
the perspective of stability. We consider the Lotka-Volterra-type (LV)
equation for species abundance on plant-animal mutualistic networks, with the
mutualistic interaction constructed from an empirical dataset web , which is
heterogeneous, and all-to-all intra-group competition assumed Maeng  _et al._
(2012); Suweis  _et al._ (2013); Stone (2016); Pascual-García and Bastolla
(2017); Barbier  _et al._ (2018); Gracia-Lázaro  _et al._ (2018); Wang  _et
al._ (2021); Maeng  _et al._ (2019); Cai  _et al._ (2020). The strengths of
mutualism and competition are set to be uniform. We restrict ourselves to the
stationary state and study the stable fixed points. Exponentially many fixed
points exist with zero components at different species, but only the stable
one is relevant to the stationary state.

To find the stable fixed point, we first show that stability can be assessed
by the signs of the non-zero components of the considered fixed point and its
neighboring fixed points. Next, approximating the adjacency matrix to be in
factorized form, we derive the non-zero components of each fixed point to be
formulated in terms of degrees and the rescaled interaction strengths. Using
these results, we devise an algorithm to classify species into surviving and
extinct ones and thereby formulate the stable fixed point, which turns out to
work well as supported by good agreement with numerical solutions. The
extinction or the diverging abundances of selected species happens depending
on parameters, the analytic understanding of which allows us to obtain the
phase diagram, including the full coexistence, selective extinction, and
unstable phase. In the selective extinction phase, small-degree species go
extinct, which results in reducing the effective interaction among the
surviving species and suppressing the onset of instability. Our study enables
a principled discrimination between surviving and extinct species and the
prediction of the abundances of the surviving species, helping us to
understand the interplay of stability, species abundance, and extinction in
structured ecological communities.

##  II Model

We consider a system of $N^{\rm(P)}$ flowering plant species and $N^{\rm(A)}$
pollinating animal species. Their abundances $x_{i}$’s evolve with time under
the LV equation as

${dx_{i}(t)\over dt}=x_{i}\left(\alpha_{i}+\sum_{j=1}^{S}B_{ij}x_{j}\right),$
(1)

where $S=N^{\rm(P)}+N^{\rm(A)}$ is the total number of species, $\alpha_{i}=1$
is the intrinsic growth rate, and $\mathbb{B}=(B_{ij})$ is the $S\times S$
interaction matrix

$\mathbb{B}\equiv-\mathbb{I}-c(\mathbb{J}^{(0)}-\mathbb{I})+m\mathbb{A}.$ (2)

Here $\mathbb{I}\equiv\begin{pmatrix}{\bf I}^{\rm(PP)}&{\bf 0}\\\ {\bf 0}&{\bf
I}^{\rm(AA)}\end{pmatrix}={\bf I}^{\rm(PP)}\oplus{\bf I}^{\rm(AA)}$ is the
identify matrix representing intraspecific regulation,
$\mathbb{J}^{(0)}\equiv\begin{pmatrix}{\bf J}^{\rm(PP)}&{\bf 0}\\\ {\bf
0}&{\bf J}^{\rm(AA)}\end{pmatrix}={\bf J}^{\rm(PP)}\oplus{\bf J}^{\rm(AA)}$
consists of the matrices of $1$’s ($J_{pp^{\prime}}=J_{aa^{\prime}}=1$ for all
$p,p^{\prime},a,a^{\prime}$) representing all-to-all competition among plants
and among pollinators along with strength $0<c<1$, and
$\mathbb{A}\equiv\begin{pmatrix}{\bf 0}&{\bf A}^{\rm(PA)}\\\ {\bf
A}^{\rm(AP)}&{\bf 0}\end{pmatrix}={\bf A}^{\rm(PA)}\oplus{\bf A}^{\rm(AP)}$ is
the symmetric adjacency matrix ($A_{pa}=A_{ap}=0,1$) with ${\bf
A}^{\rm(AP)}={{\bf A}^{\rm(PA)}}^{\intercal}$ representing the mutualistic
interaction along with the mutualism strength $m>0$. The useful properties of
the matrices of $1$’s, which are given in Appendix A, enable various analytic
calculations. There are $L\equiv\sum_{p,a}A_{pa}$ mutualistic partner pairs.

We select a real-world community Arroyo  _et al._ (1982) in a database web  to
construct the adjacency matrix $\mathbb{A}$ and use it to define $\mathbb{B}$
by Eq. (2), and build all our theoretical framework, which will be applied
later to other communities. Notice that the elements of $\mathbb{B}$ are not
random numbers but represent the interaction relationships among different
species with a uniform interaction strength $c$ or $m$. The whole interaction
network encoded in $\mathbb{B}$ and the distributions of the mutualism
degrees, $k_{p}\equiv\sum_{a}A_{pa}$ of plants and $k_{a}\equiv\sum_{p}A_{pa}$
of animals, are presented in Figs. 1(a) and 1(b), respectively. Different
degrees of species are a fundamental heterogeneity in their mutualistic
interaction profiles, which have been neglected in the random-interaction
model assuming the interaction strength between each pair of species to be an
independent and identically distributed random variable May (1972); Goh
(1979); Allesina and Tang (2015); Stone (2016); Grilli  _et al._ (2017); Bunin
(2017); Cui  _et al._ (2020); Pettersson  _et al._ (2020a); *Pettersson2020,
but they are of main concern in the present study.

Figure 1: Interaction network and species abundance of the selected community.
(a) Interaction network of $N^{\rm(P)}=43$ plants species (green triangle) and
$N^{\rm(A)}=64$ animal species (red circle) connected by $L=196$ mutualistic
interaction links (blue) and all-to-all intra-group competition links (gray).
Nodes of light red and light green represent extinct species for $c=0.1$ and
$m=0.2$. (b) Mutualism degree distributions for plants (triangle) and animals
(circle). (c) Abundances of individual species (different lines) for $c=0.1$
and $m=0.2$. (d) The stationary-state abundance vs. degree for animal species
with $c=0.1$ and selected $m$’s.

Integrating Eq. (1) up to $T=10^{3}$ with the initial condition $x_{i}(0)=1$,
we find, as shown in Fig. 1 (c), that $x_{i}(t)$’s for different species $i$
exhibit different behaviors as functions of time. They become stationary in
the long-time limit and we approximate the stationary-state abundance by the
species abundance at the final time step $T$,

$x_{i}^{\rm(st)}\equiv\lim_{t\to\infty}x_{i}(t)\simeq x_{i}(T).$ (3)

Some species show $x_{i}^{\rm(st)}=0$, implying their extinction. Also, as
shown in Fig. 1 (d), $x_{i}^{\rm(st)}$ tends to grow with degree $k_{i}$. This
correlation will be clarified in the next sections. Another remarkable feature
is that the effect of mutualism on the species abundance can be drastically
different depending on species Gracia-Lázaro  _et al._ (2018); Cai  _et al._
(2020); Wang  _et al._ (2021); The species with large degrees find their
abundances increasing with $m$ but the abundances of the species having small
degrees decrease with $m$ [Fig. 1 (d)]. In the next sections, we develop the
analytic approach to understand the nature and origin of such heterogeneity in
the species abundance depending on parameters.

##  III Stability

The state of a dynamical system, like Eq. (1), is expected to converge to a
stable fixed point in the long-time limit. If a stable fixed point exists and
is unique, its components will give the stationary-state abundances that we
obtain numerically.

Depending on which components are zero, there are $2^{S}$ different fixed
points of Eq. (1); A fixed point $\vec{x}^{*}=(x_{i}^{*})$ has components

$x_{i}^{*}=\begin{cases}0&\ {\rm for}\ i\in\mathcal{S}^{(0)},\\\
-\sum_{j\in\mathcal{S}^{(+)}}((B^{(+)})^{-1})_{ij}&\ {\rm for}\
i\in\mathcal{S}^{(+)},\end{cases}$ (4)

where $\mathcal{S}^{(0)}$ and $\mathcal{S}^{(+)}$ are the set of the species
with zero and non-zero components in the considered fixed point, respectively,
and $((B^{(+)})^{-1})_{ij}$ is the inverse of the effective interaction matrix
$\mathbb{B}^{(+)}$ obtained by eliminating the rows and columns of the species
of $\mathcal{S}^{(0)}$ in $\mathbb{B}$ Pettersson  _et al._ (2020a);
*Pettersson2020. We keep the indices $i$ or $j$ of the original interaction
matrix $\mathbb{B}$ such that $B^{(+)}_{ij}=B_{ij}$ as long as
$i,j\in\mathcal{S}^{(+)}$.

The fixed point in Eq. (4) is stable if a small perturbation $\delta
x_{i}=x_{i}-x_{i}^{*}$ does not grow persistently with time but vanishes in
the long-time limit. The time-evolution of the perturbation is given by
${d\over dt}\delta x_{i}=H_{ij}\delta x_{j}$, which involves the Jacobian
matrix at the fixed point $\vec{x}^{*}=(x_{i}^{*})$ given by

$H_{ij}=\delta_{ij}\left[1+\sum_{\ell=1}^{S}B_{i\ell}x_{\ell}^{*}\right]+x_{i}^{*}B_{ij}.$
(5)

For $i\in\mathcal{S}^{(0)}$, it holds that $x_{i}^{*}=0$, and therefore one
can see that
$H_{ij}=\delta_{ij}\left[1+\sum_{\ell=1}^{S}B_{i\ell}x_{\ell}^{*}\right]$. For
$i\in\mathcal{S}^{(+)}$, it holds that
$1+\sum_{\ell}B_{i\ell}x_{\ell}^{*}=1+\sum_{\ell\in\mathcal{S}^{(+)}}B^{(+)}_{i\ell}x_{\ell}^{*}=0$,
leading to $H_{ij}=x_{i}^{*}B_{ij}$. If all the eigenvalues $\lambda_{i}$’s of
$\mathbb{H}=(H_{ij})$ have negative real parts, then the small perturbation
will die out and the fixed point can be considered as stable.

We derive the approximate expression for the eigenvalues of $\mathbb{H}$. For
$i\in\mathcal{S}^{(0)}$, let us consider a neighboring fixed point
$\vec{x}^{*\prime}=(x_{\ell}^{*\prime})$ with
$\mathcal{S}^{(+)\prime}=\mathcal{S}^{(+)}\cup\\{i\\}$, which satisfies
$1+\sum_{\ell}B_{i\ell}x_{\ell}^{*\prime}=1+\sum_{\ell\in\mathcal{S}^{(+)}}B_{i\ell}x_{\ell}^{*\prime}+B_{ii}x_{i}^{*\prime}=0$.
Assuming that $x_{\ell}^{*\prime}\simeq x_{\ell}^{*}$ for
$\ell\in\mathcal{S}^{(+)}$ and using $B_{ii}=-1$, we find
$H_{ij}\simeq\delta_{ij}x_{i}^{*\prime}$. Therefore, with rows and columns
rearranged, the Jacobian matrix $\mathbb{H}$ contains one zero submatrix, say,
$H^{(0+)}_{ij}=0$, and three non-zero block submatrices $H^{(00)}_{ij}\simeq
x_{i}^{*\prime}\delta_{ij},H^{(+0)}_{ij}=x_{i}^{*}B_{ij}$, and
$H^{(++)}_{ij}=x_{i}^{*}B^{(+)}_{ij}$ such that

$\mathbb{H}=\begin{pmatrix}{\bf H}^{(00)}&{\bf 0}\\\ {\bf H}^{(+0)}&{\bf
H}^{(++)}\end{pmatrix}=\begin{pmatrix}(\delta_{ij}x_{i}^{*\prime})&{\bf 0}\\\
(x_{i}^{*}B_{ij})&(x_{i}^{*}B^{(+)}_{ij})\end{pmatrix}.$ (6)

Given the zero block submatrix, all the eigenvalues of $\mathbb{H}$ come from
the diagonal blocks, $H^{(00)}_{ij}\simeq x_{i}^{*\prime}\delta_{ij}$ and
$H^{(++)}_{ij}=x_{i}^{*}B^{(+)}_{ij}$. ${\bf H}^{(00)}$ is already
diagonalized, with $x_{i}^{*\prime}$’s as its eigenvalues. To obtain the
eigenvalues of ${\bf H}^{(++)}$, we decompose it as
$H^{(++)}_{ij}=-x_{i}^{*}\delta_{ij}+V_{ij}$ with $V_{ij}\equiv-
cx_{i}^{*}(1-\delta_{ij})+mx_{i}^{*}A_{ij}$ and apply the perturbation theory
with $V_{ij}$ taken as a perturbation to obtain the approximate eigenvalues
$-x_{i}^{*}$’s for small $c$ and $m$ as described in Appendix B and Ref. Stone
(2020). Therefore we find that the eigenvalues $\lambda_{i}$’s of $\mathbb{H}$
are approximately

$\lambda_{i}\simeq\begin{cases}x_{i}^{*\prime}&\ {\rm for}\
i\in\mathcal{S}^{(0)},\\\ -x_{i}^{*}&\ {\rm for}\
i\in\mathcal{S}^{(+)}.\end{cases}$ (7)

A concrete example of constructing the Jacobian matrix and deriving Eq. (7)
for a small community is presented in Appendix B.

Using Eq. (7), we can see that all the eigenvalues are negative and the fixed
point in Eq. (4) is stable when the following conditions are met: i) every
species $i$ that would have a negative fixed-point abundance
($x_{i}^{*\prime}<0$) if it were added to $\mathcal{S}^{(+)}$ has zero
abundance ($x_{i}^{*}=0$) and is in $\mathcal{S}^{(0)}$, and ii) every species
$i$ in $\mathcal{S}^{(+)}$ has a positive fixed-point abundance
($x_{i}^{*}>0$). If the components of a fixed point are known, one can use
these stability conditions to predict which species go extinct and which
species survive. In the next section, we derive the approximate analytic
formula for the components of a fixed point and use it along with Eq. (7) to
infer the stable fixed point.

##  IV Analytic approaches to the fixed point

In this section we assume that all species survive,
$\mathcal{S}^{(0)}=\emptyset$, and obtain the components of the corresponding
fixed point in Eq. (4). While the inverse of a matrix is not available in a
closed form in general, here we first consider the case of zero or weak
mutualism and then take an approximation for the adjacency matrix to derive
the components of the fixed point. The obtained results will be generalized
straightforwardly to the case of $\mathcal{S}^{(0)}\neq\emptyset$ in Sec. V.2.

###  IV.1 No mutualism

Let us consider the case of no mutualism but competition only. The interaction
matrix for $m=0$ is given in a simple form as

$\mathbb{B}=\mathbb{B}_{0}\equiv-(1-c)\mathbb{I}-c\mathbb{J}^{\rm(0)}.$ (8)

Trying $\mathbb{B}_{0}^{-1}=-{1\over 1-c}({\bf I}^{\rm(PP)}\oplus{\bf
I}^{\rm(AA)}\ominus b^{\rm(P)}{\bf J}^{\rm(PP)}\ominus b^{\rm(A)}{\bf
J}^{\rm(AA)})$ as its inverse and inserting it into
$\mathbb{B}_{0}\mathbb{B}_{0}^{-1}=\mathbb{I}$, we find that
$b^{\rm(P)}={\tilde{c}^{\rm(P)}\over N^{\rm(P)}}$ and
$b^{\rm(A)}={\tilde{c}^{\rm(A)}\over N^{\rm(A)}}$, where the rescaled
competition strengths $\tilde{c}^{\rm(P)}$ and $\tilde{c}^{\rm(A)}$ are
defined as

$\tilde{c}^{\rm(G)}\equiv{cN^{\rm(G)}\over cN^{\rm(G)}+1-c}$ (9)

with G representing either P or A. The properties of the matrix of $1$’s ,
such as ${\bf J}^{\rm(PP)}{\bf J}^{\rm(PP)}=N^{\rm(P)}{\bf J}^{\rm(PP)}$, are
used for derivation. The rescaled competition $\tilde{c}^{\rm(G)}$ ranges
between $0$ and $1$, and grows with $c$ and $N^{\rm(G)}$ as long as $0<c<1$.

The inverse matrix $\mathbb{B}_{0}^{-1}$ is represented in a compact form as

$\mathbb{B}_{0}^{-1}=(1-c)^{-1}\left(-\mathbb{I}+\mathbbm{\tilde{c}}\mathbb{\widetilde{J}}^{(0)}\right)$
(10)

with $\mathbbm{\tilde{c}}\equiv\tilde{c}^{\rm(P)}{\bf
I}^{\rm(PP)}\oplus\tilde{c}^{\rm(A)}{\bf I}^{\rm(AA)}$ and
$\mathbb{\widetilde{J}}^{(0)}\equiv{{\bf J}^{\rm(PP)}\over
N^{\rm(P)}}\oplus{{\bf J}^{\rm(AA)}\over N^{\rm(A)}}$. Then the component of
the fixed point in Eq. (4) with $\mathcal{S}^{(0)}=\emptyset$ is given by

$x_{i}^{*(0)}=x_{0}^{({\rm G}_{i})}\equiv{1-\tilde{c}^{({\rm G}_{i})}\over
1-c}={1\over cN^{({\rm G}_{i})}+1-c},$ (11)

where ${\rm G}_{i}$ is the group the species $i$ belongs to, either P or A.
The superscript $(0)$ means the zeroth-order approximation. The first-order
correction will be presented as well in the next subsection and then we move
to the approximation for general $m$. For $0<c<1$, Eq. (11) is positive for
all $i$. Therefore all the eigenvalues of the Jacobian in Eq. (7) with
$\mathcal{S}^{(0)}=\emptyset$ are negative, and the fixed point with Eq. (11)
for all $i$ is the stable fixed point. The increase of the competition
strength $c$ or of the number of species $N^{({\rm G}_{i})}$ leads to the
decrease of the abundance $x_{i}^{*(0)}$.

###  IV.2 Weak mutualism

Suppose that the mutualism strength $m$ is not zero but small. Expanding the
inverse $\mathbb{B}^{-1}$ as

$\mathbb{B}^{-1}=\left(\mathbb{B}_{0}+m\mathbb{A}\right)^{-1}=\mathbb{B}_{0}^{-1}\sum_{n=0}^{\infty}\left(-m\mathbb{A}\mathbb{B}_{0}^{-1}\right)^{n}$
(12)

with $\mathbb{B}_{0}$ given in Eq. (8), and utilizing the relations like ${\bf
A}^{\rm(PA)}{\bf J}^{\rm(AA)}={\bf K}^{\rm(PP)}{\bf J}^{\rm(PA)}$ with
$K_{pp^{\prime}}=k_{p}\delta_{pp^{\prime}}$ a block in the degree matrix

$\mathbb{K}\equiv{\bf K}^{\rm(PP)}\oplus{\bf K}^{\rm(AA)},$ (13)

one can evaluate the first order term in $m$ in Eq. (12) to obtain the first-
order approximation

$x_{i}^{*(1)}\simeq x_{0}^{({\rm G}_{i})}\left[1+{m\over
1-c}{1-\tilde{c}^{({\rm\bar{G}}_{i})}\over 1-\tilde{c}^{({\rm
G}_{i})}}\left(k_{i}-\tilde{c}^{({\rm G}_{i})}\langle k\rangle^{({\rm
G}_{i})}\right)\right]$ (14)

with the mean degree $\langle k\rangle^{\rm(G)}\equiv{L\over N^{\rm(G)}}$.
This first-order approximation works if $||\mathbb{A}\mathbb{B}_{0}^{-1}||$ is
sufficiently small. See Appendix C for more details.

The formula in Eq. (14) allows us to understand the origin of the ambivalent
effects of mutualism on the species abundance as observed in Fig. 1 (d). The
two terms in the parentheses in Eq. (14) represent the mutualistic benefit of
a species $i$ (in group ${\rm G}_{i}$) from its $k_{i}$ mutualistic partners
in group $\bar{\rm G}_{i}$, and the competition with other species in the same
group ${\rm G}_{i}$ that also benefit from mutualism, respectively. Their
difference may be positive or negative depending on degree $k_{i}$. It is the
species with $k_{i}>\tilde{c}^{\rm(G)}\langle k\rangle^{({\rm G}_{i})}$ that
finds abundance increasing with increasing $m$; the abundance of the species
with $k_{i}<\tilde{c}^{\rm(G)}\langle k\rangle^{({\rm G}_{i})}$ decreases with
$m$, as its mutualistic benefit is overwhelmed by the competition with the
species in the same group to the extent proportional to $m$ in the first-order
approximation.

One caveat is that Eq. (14) can be negative depending on parameters, which
suggests that the fixed point with all species surviving is unstable and that
some species will turn out to have zero abundance in the stable fixed point.
This will be explored in Sec. V.2.

###  IV.3 Annealed approximation for general $m$

Each term for $n\geq 1$ in Eq. (12) represents the sum of the influences of
other species on the abundance of a species built up over the pathways
involving $n$ mutualistic pairs. To analytically track such higher-order
contributions, we consider the annealed adjacency matrix
$\tilde{A}_{pa}\equiv{k_{p}k_{a}\over L}$, meaning the probability to connect
$p$ and $a$ in the network ensemble for a given degree sequence Lee  _et al._
(2009), and equivalently

$\mathbb{\widetilde{A}}=L^{-1}\mathbb{K}\mathbb{J}^{(1)}\mathbb{K},$ (15)

where $\mathbb{K}$ is the degree matrix introduced in Eq. (13) and
$\mathbb{J}^{(1)}$ contains the matrices of $1$’s at the off-diagonal blocks
as $\mathbb{J}^{(1)}\equiv{\bf J}^{\rm(PA)}\oplus{\bf J}^{\rm(AP)}$ with
$J_{pa}=J_{ap}=1$ for all $p$ and $a$. Then, after some algebra utilizing the
properties of the matrices of $1$’s as detailed in Appendix D, we find each
term $\mathbb{B}_{0}^{-1}(-m\mathbb{\widetilde{A}}\mathbb{B}_{0}^{-1})^{n}$
reduced to

$\displaystyle\mathbb{B}_{0}^{-1}\left(-m\mathbb{\widetilde{A}}\mathbb{B}_{0}^{-1}\right)^{n}$
$\displaystyle=\begin{cases}-{\tilde{m}^{n}\over
1-c}\mathbb{\widetilde{K}}\mathbb{\widetilde{J}}^{(1)}\mathbb{\widetilde{K}}&\
{\rm for}\ n=1,3,5,\ldots,\\\ -{\tilde{m}^{n}\over
1-c}\mathbb{\widetilde{K}}\mathbb{\widetilde{J}}^{(0)}\mathbb{\widetilde{K}}&\
{\rm for}\ n=2,4,6,\ldots,\end{cases}$ (16)

where we introduced $\mathbb{\widetilde{J}}^{(1)}\equiv{{\bf
J}^{\rm(PA)}\over\sqrt{N^{\rm(P)}N^{\rm(A)}}}\oplus{{\bf
J}^{\rm(AP)}\over\sqrt{N^{\rm(P)}N^{\rm(A)}}}$,
$\mathbb{\widetilde{K}}\equiv{{{\bf K}^{\rm(PP)}\over\langle
k\rangle^{\rm(P)}}-\tilde{c}^{\rm(P)}{\bf
I}^{\rm(PP)}\over\sqrt{\xi^{\rm(P)}-\tilde{c}^{\rm(P)}}}\oplus{{{\bf
K}^{\rm(AA)}\over\langle k\rangle^{\rm(A)}}-\tilde{c}^{\rm(A)}{\bf
I}^{\rm(AA)}\over\sqrt{\xi^{\rm(A)}-\tilde{c}^{\rm(A)}}}$, and the rescaled
mutualism strength $\tilde{m}$

$\tilde{m}\equiv\frac{m}{1-c}\sqrt{\langle k\rangle^{(P)}\langle
k\rangle^{(A)}(\xi^{(P)}-\tilde{c}^{(P)})(\xi^{(A)}-\tilde{c}^{(A)})}$ (17)

with the ratio of the first two moments of the mutualism degree

$\xi^{\rm(G)}\equiv{\langle k^{2}\rangle^{\rm(G)}\over{\langle
k\rangle^{\rm(G)}}^{2}}$ (18)

quantifying the heterogeneity of degree Yan  _et al._ (2017). The rescaled
mutualism $\tilde{m}$ is the key parameter governing the species abundance,
capturing the effects of network structural heterogeneity on the species
abundance.

All the terms for $n\geq 1$ in Eq. (12) are proportional to either
$\mathbb{\widetilde{J}}^{(0)}$ or $\mathbb{\widetilde{J}}^{(1)}$, with
$\tilde{m}^{n}$ in the coefficient. Consequently, the sum of the influences of
interspecific interactions over all possible pathways in Eq. (12) is reduced
to two infinite geometric series, manifesting the advantage of the annealed
approximation. Then the inverse matrix is expressed in a closed form as

$\mathbb{\widetilde{B}}^{-1}=\mathbb{B}_{0}^{-1}-{1\over 1-c}{\tilde{m}\over
1-\tilde{m}^{2}}\mathbb{\widetilde{K}}\left(\tilde{m}\mathbb{\widetilde{J}}^{(0)}+\mathbb{\widetilde{J}}^{(1)}\right)\mathbb{\widetilde{K}}.$
(19)

Substituting Eq. (19) in Eq. (4), we obtain the fixed point

$\displaystyle\tilde{x}^{*}_{i}=x_{0}^{({\rm
G}_{i})}\left(1+\bar{k}_{i}\tilde{m}{\tilde{m}+\eta^{({\rm G}_{i}\bar{\rm
G}_{i})}\over 1-\tilde{m}^{2}}\right)$ (20)

with the rescaled degree

$\bar{k}_{i}\equiv{{k_{i}\over\langle k\rangle^{({\rm
G}_{i})}}-\tilde{c}^{({\rm G}_{i})}\over\xi^{({\rm G}_{i})}-\tilde{c}^{({\rm
G}_{i})}},$ (21)

and the asymmetry factor

$\eta^{\rm(PA)}\equiv{1-\tilde{c}^{\rm(A)}\over
1-\tilde{c}^{\rm(P)}}\sqrt{\langle
k\rangle^{\rm(P)}(\xi^{\rm(P)}-\tilde{c}^{\rm(P)})\over\langle
k\rangle^{\rm(A)}(\xi^{\rm(A)}-\tilde{c}^{\rm(A)})}={1\over\eta^{\rm(AP)}}.$
(22)

The formula in Eq. (20) is the main result of the present work, representing
the abundance of individual species under heterogeneous mutualistic
interactions and uniform intra-group competition. It is the cornerstone of the
results that follow in the next sections. The abundance is given by a non-
linear function of the rescaled mutualism $\tilde{m}$ in Eq. (17), revealing
how the higher-order contributions of interspecific interactions are combined
with the network structure. The increase of mean connectivity $\langle
k\rangle^{\rm(P,A)}$ or the increase of the degree heterogeneity
$\xi^{\rm(P,A)}$ enhances the rescaled mutualism strength. As $\tilde{m}$
increases, $\tilde{x}_{i}^{*}$ may increase or decrease, depending on the sign
of the rescaled degree $\bar{k}_{i}$. The rescaled degree quantifies the
imbalance of the mutualism benefit and the competition cost;
$\tilde{x}_{i}^{*}$ increases (decreases) with $\tilde{m}$ if $\bar{k}_{i}$ is
positive (negative) as long as $0<\tilde{m}<1$.

One can notice that $\tilde{x}_{i}^{*}$ in Eq. (20) can be negative for some
species $i$ depending on parameters and degree, implying then that Eq. (20) is
not the stable fixed point and invoking the necessity to classify correctly
surviving and extinct species. The divergence of Eq. (20) at $\tilde{m}=1$
suggests the onset of instability, which can be suppressed up to a larger
value of $\tilde{m}$ than one, along with the extinction of selected species
as we will see.

##  V Phase diagram

Some of the formulated abundances in Eq. (20) can be negative depending on
parameters. Then, by Eq. (7), some species may have to have zero abundance and
the remaining surviving species should have positive abundances different from
Eq. (20) as the interspecific interaction among the surviving species should
be considered to formulate their abundances. In this section, we investigate
the stable fixed point depending on parameters by identifying the species to
go extinct, if any, and recalculating the abundance of the surviving species,
and we obtain the phase diagram.

###  V.1 Full coexistence phase

Let us call it full coexistence if there is no extinct species, i.e., if
$x_{i}^{\rm(st)}>0$ for all $i$. For sufficiently small $\tilde{m}$, the
numerically and analytically obtained values for the stationary-state
abundance, $x_{i}^{\rm(st)}$ and $\tilde{x}_{i}^{*}$ in Eq. (20), are in good
agreement. This agreement implies that in the full-coexistence phase i) the
annealed approximation works, $x_{i}^{\rm(st)}\simeq\tilde{x}_{i}^{\rm(st)}$,
and ii) Eq. (20) is stable, $\tilde{x}_{i}^{\rm(st)}=\tilde{x}_{i}^{*}$, where
$\tilde{x}_{i}^{\rm(st)}$ is the stationary-state abundance from the solution
to Eq. (1) with the annealed adjacency matrix $\tilde{A}_{ij}$ used. See
Appendix E for different kinds of species abundances used in this paper.

From Eq. (7), the full-coexistence fixed point in Eq. (20) is stable only if
all $\tilde{x}_{i}^{*}$’s are positive. This holds for
$\tilde{m}<\tilde{m}_{e}^{*}(c)\equiv\min\left(\tilde{m}_{e}^{\rm*(P)}(c),\tilde{m}_{e}^{\rm*(A)}(c)\right)$
with

$\displaystyle\tilde{m}_{e}^{\rm*(G)}(c)\equiv$
$\displaystyle\left\\{\begin{array}[]{ll}1&\ {\rm for}\ c<c_{\rm
min}^{\rm*(G)},\\\ {\sqrt{4\left(1-\bar{k}_{\rm
min}^{\rm(G)}\right)+\left(\bar{k}_{\rm
min}^{\rm(G)}\eta^{\rm(G\bar{G})}\right)^{2}}+\bar{k}_{\rm
min}^{\rm(G)}\eta^{\rm(G\bar{G})}\over 2\left(1-\bar{k}_{\rm
min}^{\rm(G)}\right)}&\ {\rm for}\ c>c_{\rm min}^{\rm*(G)}.\end{array}\right.$
(25)

Here $\bar{k}_{\rm min}^{\rm(G)}$ is the rescaled degree of the group-G
species having the smallest degree $k_{\rm min}^{\rm(G)}$, and $c_{\rm
min}^{\rm*(G)}$ is defined as

$c_{\rm min}^{\rm*(G)}\equiv{k_{\rm min}^{\rm(G)}\over N^{\rm(G)}\left(\langle
k\rangle^{\rm(G)}-k_{\rm min}^{\rm(G)}\right)+k_{\rm min}^{\rm(G)}},$ (26)

such that $\bar{k}_{\rm min}^{\rm(G)}<0$ for $c>c_{\rm min}^{\rm(G)}$.

Figure 2: Phase diagram and stability for the selected community Arroyo  _et
al._ (1982). (a) The phase boundaries $\tilde{m}_{e}^{\rm(st)}$ and
$\tilde{m}_{u}$ based on the stationary-state abundances are compared with the
theoretical prediction $\tilde{m}_{e}^{*}$ and $\tilde{m}_{u}^{*}$ from Eq.
(25) and the condition $\tilde{m}^{(+)}|_{\tilde{m}_{u}^{*}}=1$, respectively.
(b) The largest real part $\Lambda^{\rm(st)}$ of the eigenvalues of the
Jacobian matrix at $\vec{x}^{\rm(st)}$ and that at $\vec{x}^{*}$ approximated
as
$\tilde{\Lambda}^{*}=\max_{i\in\mathcal{S}^{(0)},j\in\mathcal{S}^{(+)}}\left(\tilde{x}_{i}^{*\prime},-\tilde{x}_{j}^{*},-1\right)$
are shown as functions of $\tilde{m}$ for $c=0.1$.

The predicted boundary of the full coexistence phase $\tilde{m}_{e}^{*}(c)$ is
shown by a dashed line in Fig. 2 (a), which is fixed at one for $c<c_{\rm
min}^{*}\equiv\min(c_{\rm min}^{\rm*(P)},c_{\rm min}^{\rm*(A)})$ and decreases
with $c$ for $c>c_{\rm min}^{*}$. For $c<c_{\rm min}^{*}$, all species have
positive rescaled degrees, $\bar{k}_{i}>0$, and their fixed-point abundances
$\tilde{x}_{i}^{*}$ increase with $\tilde{m}$ until diverging at
$\tilde{m}=\tilde{m}_{e}^{*}(c)=1$. Considering the fraction of diverging-
abundance species $r_{u}\equiv S^{-1}\sum_{i}\theta(x_{i}^{\rm(st)}-M)$ with
$M=10^{5}$, we call it the unstable phase if $r_{u}>0$. We define the
instability threshold $\tilde{m}_{u}(c)$ such that $r_{u}>0$ for
$\tilde{m}\geq\tilde{m}_{u}(c)$. One finds in Eq. (20) that the theoretical
prediction for the instability threshold is given by $\tilde{m}_{u}^{*}(c)=1$
for $c<c_{\rm min}^{*}$, which agrees with the boundary of the full
coexistence phase $\tilde{m}_{e}^{*}(c)=1$ in Eq. (25). Notice that the
unstable phase meets the full coexistence phase at $\tilde{m}=1$ for $c<c_{\rm
min}^{*}$ [Fig. 2 (a)].

When $c$ is larger than $c_{\rm min}^{*}$ and $\tilde{m}$ is larger than
$\tilde{m}_{e}^{*}(c)$, there appear some species $i$ with
$\tilde{x}_{i}^{*}<0$ according to Eq. (20), implying that they should go
extinct, having zero abundance in the true stable fixed point. Computing the
fraction of extinct species in terms of their stationary-state abundances as

$r_{e}^{\rm(st)}\equiv S^{-1}\sum_{i=1}^{S}\theta(\epsilon-x_{i}^{\rm(st)}),$
(27)

with $\epsilon=10^{-5}$ introduced under finite precision of numerics and the
Heaviside step function $\theta(x)=1$ for $x>0$ and $0$ otherwise, we find
that $r_{e}^{\rm(st)}$ becomes non-zero as $\tilde{m}$ exceeds the extinction
threshold $\tilde{m}_{e}^{\rm(st)}(c)$ for $c>c_{\rm min}^{*}$, which is well
approximated by the predicted threshold $\tilde{m}_{e}^{*}(c)$ in Eq. (25)
[Fig. 2 (a)]. Let us call it selective extinction if there exist extinct
species ($r_{e}^{\rm(st)}>0$) but no abundance-diverging species ($r_{u}=0$).
Our analysis shows that the full coexistence phase meets the selective
extinction phase at $\tilde{m}_{e}^{\rm(st)}\simeq\tilde{m}_{e}^{*}$ for
$c>c_{\rm min}^{*}$. While we showed that the full-coexistence phase ends at
$\tilde{m}_{e}^{\rm(st)}$, it remains to be explored which species go extinct
and what happens for the remaining surviving species. It will be addressed in
the next subsection in details.

In Fig. 2 (b), the largest real part $\Lambda^{\rm(st)}$ of the eigenvalues of
the Jacobian $\mathbb{H}$ computed at the stationary-state abundance
$\vec{x}^{\rm(st)}=(x^{\rm(st)}_{i})$ is shown to be negative in the full
coexistence phase, demonstrating the stability of $\vec{x}^{\rm(st)}$. The
analytically-obtained fixed point $\vec{\tilde{x}}^{*}=(\tilde{x}_{i}^{*})$ is
stable as well; The largest eigenvalue $\tilde{\Lambda}^{*}$ of the Jacobian
$\tilde{H}_{ij}=\delta_{ij}\left[1+\sum_{\ell=1}^{S}\tilde{B}_{i\ell}\tilde{x}_{\ell}^{*}\right]+\tilde{x}_{i}^{*}\tilde{B}_{ij}$
is evaluated as $\tilde{\Lambda}^{*}=\max_{i}(-\tilde{x}_{i}^{*},-1)$, from
Eq. (7) and the property that $\tilde{H}_{ij}$ has the eigenvalue $-1$
[Appendix B], and remains negative in the full-coexistence phase. Moreover, we
see a good agreement between $\Lambda$ and $\tilde{\Lambda}^{*}$.

###  V.2 Selective extinction phase

If the set $\mathcal{S}^{(0)}$ of extinct species is known, one can apply Eq.
(20) to the subcommunity of the surviving species and obtain their abundances
analytically. Removing the rows and columns of the species belonging to
$\mathcal{S}^{(0)}$ in the full matrix $\mathbb{\widetilde{B}}$ and also in
the adjacency matrix $\mathbb{\widetilde{A}}$, one can obtain the effective
interaction matrix $\mathbb{\widetilde{B}}^{(+)}$ and the effective adjacency
matrix $\mathbb{\widetilde{A}}^{(+)}$ for the surviving-species community,
from which we can compute the effective quantities such as $x_{0}^{\rm(G,+)}$,
$\bar{k}_{i}^{(+)}$, $\tilde{m}^{(+)}$, and $\eta^{\rm(PA,+)}$ as described in
Appendix F. Using them in Eq. (20), one can obtain the approximate stable
fixed point

$\tilde{x}^{*}_{i}=\begin{cases}0&\ {\rm for}\ i\in\mathcal{S}^{(0)},\\\
x_{0}^{({\rm
G}_{i},+)}\left(1+\bar{k}_{i}^{(+)}\tilde{m}^{(+)}{\tilde{m}^{(+)}+\eta^{({\rm
G}_{i}\bar{\rm G}_{i},+)}\over 1-(\tilde{m}^{(+)})^{2}}\right)&\ {\rm for}\
i\in\mathcal{S}^{(+)}.\end{cases}$ (28)

The sets of extinct and surviving species, $\mathcal{S}^{(0)}_{\rm stable}$
and $\mathcal{S}^{(+)}_{\rm stable}$ for the stable fixed point are not given
a priori. Examining all possible sets $\mathcal{S}^{(0)}$ and identifying the
one yielding all negative eigenvalues as given in Eq. (7) could be done but
takes a very long time.

Our analytic results, Eqs. (7) and (28), give a clue to proceed. Suppose that
we have a pair of disjoint sets $\mathcal{S}^{(0)}$ and $\mathcal{S}^{(+)}$
with $\mathcal{S}\equiv\mathcal{S}^{(0)}\cup\mathcal{S}^{(+)}$ the set of all
species. If a species in $\mathcal{S}^{(+)}$ with a small (large) effective
degree has a negative (positive) value of $\tilde{x}_{i}^{*}$ from Eq. (28),
then it will be likely to be in the right set $\mathcal{S}^{(0)}_{\rm stable}$
($\mathcal{S}^{(+)}_{\rm stable}$) for the stable fixed point according to Eq.
(7). Using this reasoning, we can update iteratively $\mathcal{S}^{(0)}$ and
$\mathcal{S}^{(+)}$ towards obtaining $\mathcal{S}^{(0)}_{\rm stable}$ and
$\mathcal{S}^{(+)}_{\rm stable}$ as follows.

Initially we begin with $\mathcal{S}^{(0)}=\emptyset$,
$\mathcal{S}^{(+)}=\mathcal{S}$,
$\mathbb{\widetilde{B}}^{(+)}=\mathbb{\widetilde{B}}$, and $\tilde{x}_{i}^{*}$
evaluated by Eq. (28). Then, the following procedures are repeated until
stopping at the step (iii):
(i) We label as new extinct species all the plant (animal) species $p_{e}$’s
($a_{e}$’s) with their fixed-point abundances $\tilde{x}_{p_{e}}^{\rm*}$
($\tilde{x}_{a_{e}}^{\rm*}$) having different sign from that of the hub plant
species $\tilde{x}_{p_{\rm hub}}^{*}$ (animal species $\tilde{x}_{a_{\rm
hub}}^{*}$), the one having the largest effective degree.
(ii) Go to step iv) if there are such new extinct species, or
(iii) Stop if there are none.
(iv) We remove those new extinct species from $\mathcal{S}^{(+)}$ and add them
to $\mathcal{S}^{(0)}$ and update $\mathbb{\widetilde{B}}^{(+)}$ by
eliminating their rows and columns, and
(v) Set $\tilde{x}_{p_{e}}^{*}=\tilde{x}_{a_{e}}^{*}=0$, and evaluate
$\tilde{x}_{i}^{*}$’s for the remaining surviving species $i$ by Eq. (28) with
using the new $\mathbb{\widetilde{B}}^{(+)}$.

Figure 3: Selective extinction phase. (a) Stationary-state abundances
$x_{i}^{\rm(st)}$ are compared with the stable fixed-point ones $x_{i}^{*}$’s
and $\tilde{x}_{i}^{*}$’s for $c=0.1$ and $\tilde{m}=0.8$. The species index
$i$ is arranged in the descending order of degree among plants and among
animals. Inset: Abundance vs. degree. (b) Fraction $r_{e}^{\rm(st)}$ of
extinct species and $r_{u}$ of the abundance-diverging species based on the
stationary-state abundance for $c=0.1$. They are compared with
$\tilde{r}_{e}^{*}$ and $\tilde{r}_{u}^{*}$ based on the stable fixed point
$\tilde{x}_{i}^{*}$’s and also with $r_{e}^{*}$ based on $x_{i}^{*}$’s. The
critical points are also marked, $\tilde{m}_{e}^{\rm(st)}\simeq 0.53$ and
$\tilde{m}_{e}^{*}\simeq 0.61$. Inset: The effective rescaled mutualism
strength $\tilde{m}^{(+)}$ reaches $1$ at $\tilde{m}=\tilde{m}_{u}^{*}\simeq
1.3$ for $c=0.1$. (c) The collapse of the non-zero minimum abundances $x_{\rm
min}^{\rm(st)}$ rescaled as in Eq. (29) in 46 real-world communities as
functions of $\tilde{m}^{(+)}$.

Note that the effective rescaled mutualism strength $\tilde{m}^{\rm(+)}$ can
be larger than one, making the hub abundances negative according to Eq. (28)
in the middle of the above procedures, which is why we compare the sign of the
abundance of a species with that of hub species to identify extinct species.
At the end of these procedures we are given $\mathcal{S}^{(0)}_{\rm stable}$
and $\mathcal{S}^{(+)}_{\rm stable}$, which we use to obtain the stable fixed
point $\tilde{x}_{i}^{*}$’s from Eq. (28) git .

As shown in Fig. 3 (a), the predicted abundance $\tilde{x}_{i}^{*}$
approximates reasonably the stationary-state abundance $x_{i}^{\rm(st)}$. It
grows with degree $k_{i}$ and takes a zero value for the species with the
smallest degrees, demonstrating the crucial role of degree on extinction. Its
origin can be understood by examining the rescaled degree, the imbalance of
the mutualistic benefit and the competition cost, the former of which is
proportional to the raw degree. The stable fixed point $\tilde{x}_{i}^{*}$
predicts whether a species survives or goes extinct correctly for 80% of
species across parameters; The predicted fraction $\tilde{r}_{e}^{*}\equiv
S^{-1}\sum_{i}\theta(\epsilon-\tilde{x}_{i}^{*})$ of extinct species is in
good agreement with the true value $r_{e}^{\rm(st)}$ in Eq. (27), which plays
the role of the order parameter distinguishing the full-coexistence phase
($r_{e}=0$) and the selective extinction phase ($r_{e}>0$) [Fig. 3 (b)].

Deviations stem from the annealed approximation; The stationary-state
abundance $\tilde{x}_{i}^{\rm(st)}$ under the annealed adjacency matrix agrees
perfectly with $\tilde{x}_{i}^{*}$ (Appendix G), and the predicted fraction of
extinct species $\tilde{r}_{e}^{*}$ explains very well the true value
$\tilde{r}_{e}$ under the annealed adjacency matrix (Appendix H). Instead of
using $\mathbb{\widetilde{B}}$ and Eq. (28), one can also use $\mathbb{B}$ and
Eq. (4) in the above procedures to identify the set of extinct and surviving
species and obtain the stable fixed point $x_{i}^{*}$ from Eq. (4), which
agrees very well with $x_{i}^{\rm(st)}$ as shown in Fig. 3 (a).

One might have expected instability to arise at $\tilde{m}=1$ from Eq. (20).
However the extinction of the small-degree species effectively reduces
$\tilde{m}$ to $\tilde{m}^{\rm(+)}$ [Fig. 3 (b)] and the abundances of the
surviving species are evaluated by Eq. (28). The effective rescaled mutualism
strength $\tilde{m}^{(+)}$ is kept smaller than $1$, preventing the onset of
instability, up to $\tilde{m}_{u}(c)$ for $c>c_{\rm min}^{*}$. As the
smallest-degree species go extinct, we find that the degree heterogeneity is
reduced in the interaction network of the surviving species, which drives the
reduction of the effective rescaled mutualism strength (Appendix F). Like in
the full coexistence phase, the largest real part of the eigenvalues $\Lambda$
at $\vec{x}^{\rm(st)}$ and
$\tilde{\Lambda}^{*}=\max_{i\in\mathcal{S}^{(0)}_{\rm
stable},j\in\mathcal{S}^{(+)}_{\rm
stable}}\left(\tilde{x}_{i}^{*\prime},-\tilde{x}_{j}^{*},-1\right)$ at the
stable fixed point $\tilde{x}_{i}^{*}$ remain negative in the selective
extinction phase, demonstrating stability [Fig. 2 (b)].

The selective extinction phase meets the unstable phase at $\tilde{m}_{u}(c)$
where the fraction of abundance-diverging species $r_{u}$ becomes non-zero.
The instability threshold $\tilde{m}_{u}(c)$ coincides with its theoretical
prediction $\tilde{m}_{u}^{*}$ at which $\tilde{m}^{(+)}=1$ and
$\tilde{x}_{i}^{*}$ in Eq. (28) diverges [Fig. 3 (b)]. For
$\tilde{m}>\tilde{m}_{u}^{*}$, the non-zero components of the stable fixed
point becomes negative, featuring the non-zero fraction
$\tilde{r}_{u}^{*}\equiv S^{-1}\sum_{i}\theta(-\epsilon-\tilde{x}_{i}^{*})$ of
the species having negative $\tilde{x}_{i}^{*}$.

Lastly, to demonstrate the general applicability of our analytic results, we
study the minimum stationary-state abundance of the surviving species $x_{\rm
min}^{\rm(st)}$ in 46 large empirical mutualistic networks with
$N^{\rm(P)},N^{\rm(A)}\geq 20$ web . As in Sec. II, we assume the uniform
intra-group competition with strength $c$, and we use the data-sets to
construct the mutualism adjacency matrices with the mutualism strength $m$.

From the theoretical framework developed in the previous sections, $x_{\rm
min}^{\rm(st)}$ can be approximated by the non-zero minimum component of the
stable fixed point, $\tilde{x}_{\rm min}^{*}=\min_{i\in\mathcal{S}^{(+)}_{\rm
stable}}\tilde{x}_{i}^{*}=\tilde{x}_{i_{\rm min}}^{*}$ representing the
predicted abundance of species $i_{\rm min}$. From Eq. (28), we find that
$\tilde{x}_{\rm min}^{*}$ behaves as a function of $\tilde{m}^{(+)}$ as

${{\tilde{x}_{\rm min}^{*}\over x_{0,{\rm min}}^{(+)}}-1\over\bar{k}_{\rm
min}^{(+)}\eta_{\rm min}^{(+)}}\simeq{\tilde{m}^{(+)}\over
1-\tilde{m}^{(+)}},$ (29)

where $\bar{k}_{\rm min}^{(+)}\equiv\bar{k}_{i_{\rm min}}$, $x_{0,{\rm
min}}^{(+)}\equiv x_{0}^{({\rm G}_{i_{\rm min}},+)}$, and $\eta_{\rm
min}^{(+)}\equiv\eta^{({\rm G}_{i_{\rm min}}\bar{\rm G}_{i_{\rm min}},+)}$,
and we approximate $1+\tilde{m}^{(+)}/\eta^{(+)}_{\rm min}$ by
$1+\tilde{m}^{(+)}$ in the right-hand-side. In Fig. 3 (c), the plots of the
rescaled minimum abundance given by the left-hand-side of Eq. (29) with
$x_{\rm min}^{\rm(st)}$ in place of $\tilde{x}_{\rm min}^{*}$ versus
$\tilde{m}^{(+)}$ for $c=0.1$ in 46 empirical mutualistic networks collapse
reasonably onto ${\tilde{m}^{\rm(+)}\over 1-\tilde{m}^{(+)}}$ in agreement
with Eq. (29). There are outliers, though. About 6% of the data points have
their rescaled minimum abundances negative and they are thus neglected in Fig.
3 (c). Some of the outliers seen in Fig. 3 (c) are attributed to the annealed
approximation, which return close to the theoretical curve in the annealed
network [see Fig. 6 (e)]. Outliers seen for $\tilde{m}^{(+)}\simeq 1$ in Fig.
6 (e) suggest the possible relevance of network characteristics beyond the
degree sequence, which needs further investigation.

##  VI Summary and discussion

While various theoretical approaches have been established for studying the
stability and biodiversity of random unstructured communities, the relation
between the structured interaction, universal in the real world, and the
emergent phenomena of the community has been little understood, partly due to
the lack of an analytically tractable model. Here we considered a model
community of two groups of species - plants and pollinators - under uniform
intra-group competition and empirical heterogeneous inter-group mutualism, and
we investigated analytically and numerically the abundance and extinction of
individual species in that community.

Deriving the stability condition and the stable fixed point of the LV
equation, we quantified the influences of the structural heterogeneity. The
strength of mutualism is rescaled by the degree heterogeneity. The species
with few mutualistic partners are driven to extinction by their little benefit
of mutualism compared with the high cost of competition. As the mutualism
strength increases, more species find benefit falling short of cost, resulting
in the increase of the number of extinct species. The effective rescaled
mutualism among the surviving species is reduced with respect to the original
one, which enables the community of the surviving species to be stable for a
wide range of parameters, delaying the onset of instability. The number of
extinct species and the fraction of the abundance-diverging species play the
roles of the order parameter in the phase diagram.

Going beyond the annealed approximation to identify further network
characteristics than degree may provide rich concepts and methods to
characterize the structure-function relationship of ecological communities.
Contrary to the unstructured communities where the interspecific interaction
patterns and strengths are random and distributed identically across species,
the number of mutualistic partners is different from species to species in the
structured community that we study in the present work. We showed how the
abundance and the likelihood of extinction depend on the degree of a species.
Real-world communities should exhibit both non-uniform interaction strengths
and heterogeneous connection patterns. If our analytic framework can be
generalized to handle not only heterogeneity but also the randomness of the
interaction matrix, it will help us to better understand how structural
heterogeneity and randomness together govern the stability and species
extinction of real-world communities.

## Acknowledgment

We thank Matthieu Barbier, Hye Jin Park, Yongjoo Baek, and Sang Hoon Lee for
valuable comments. This work was supported by the National Research Foundation
of Korea (NRF) grants funded by the Korean Government (No. 2019R1A2C1003486
(D.-S.L) and No. 2020R1A2C1005334 (J.W.L)) and a KIAS Individual Grant (No.
CG079901) at Korea Institute for Advanced Study (D.-S.L).

##  Appendix A Properties of the matrices of $1$’s

The matrix of $1$’s denoted by ${\bf J}$ having all elements equal to $1$,
$J_{ij}=1$ for all $i$ and $j$, is used in the present work to represent the
all-to-all uniform competition among plants and among animals via ${\bf
J}^{\rm(PP)}$ and ${\bf J}^{\rm(AA)}$ of dimension $N^{\rm(P)}\times
N^{\rm(P)}$ and $N^{\rm(A)}\times N^{\rm(A)}$, respectively, and also to
represent the uniform coupling between plants and animals, appearing in the
annealed adjacency matrix, via ${\bf J}^{\rm(PA)}$ and ${\bf J}^{\rm(AP)}$ of
dimension $N^{\rm(P)}\times N^{\rm(A)}$ and $N^{\rm(A)}\times N^{\rm(P)}$,
respectively. We also consider its integrated versions in block-matrix form

$\displaystyle\mathbb{J}^{(0)}$ $\displaystyle\equiv\begin{pmatrix}{\bf
J}^{\rm(PP)}&{\bf 0}\\\ {\bf 0}&{\bf J}^{\rm(AA)}\end{pmatrix}={\bf
J}^{\rm(PP)}\oplus{\bf J}^{\rm(AA)},$ $\displaystyle\mathbb{J}^{(1)}$
$\displaystyle\equiv\begin{pmatrix}{\bf 0}&{\bf J}^{\rm(PA)}\\\ {\bf
J}^{\rm(AP)}&{\bf 0}\end{pmatrix}={\bf J}^{\rm(PA)}\oplus{\bf J}^{\rm(AP)},$
(30)

with $\oplus$ representing the direct sums of two matrices of $1$’s defined on
different groups of nodes, and the rescaled matrices given by

$\displaystyle\mathbb{\widetilde{J}}^{(0)}$
$\displaystyle\equiv\begin{pmatrix}{{\bf J}^{\rm(PP)}\over N^{\rm(P)}}&{\bf
0}\\\ {\bf 0}&{{\bf J}^{\rm(AA)}\over N^{\rm(A)}}\end{pmatrix}={{\bf
J}^{\rm(PP)}\over N^{\rm(P)}}\oplus{{\bf J}^{\rm(AA)}\over N^{\rm(A)}},$
$\displaystyle\mathbb{\widetilde{J}}^{(1)}$
$\displaystyle\equiv\begin{pmatrix}{\bf 0}&{{\bf
J}^{\rm(PA)}\over\sqrt{N^{\rm(P)}N^{\rm(A)}}}\\\ {{\bf
J}^{\rm(AP)}\over\sqrt{N^{\rm(P)}N^{\rm(A)}}}&{\bf 0}\end{pmatrix}={{\bf
J}^{\rm(PA)}\oplus{\bf J}^{\rm(AP)}\over\sqrt{N^{\rm(P)}N^{\rm(A)}}}.$ (31)

In this appendix, we present their useful properties, which are used to derive
the analytic results presented in the main text.

Let us denote the matrices of $1$’s of dimension $N_{1}\times N_{2}$ by ${\bf
J}^{(N_{1}\times N_{2})}$. If one multiplies ${\bf J}^{(N_{1}\times N_{2})}$
and ${\bf J}^{(N_{2}\times N_{3})}$, then she obtains ${\bf J}^{(N_{1}\times
N_{2})}{\bf J}^{(N_{2}\times N_{3})}=N_{2}{\bf J}^{\rm(N_{1}\times N_{3})},$
since $\sum_{j=1}^{N_{2}}J_{ij}^{\rm(N_{1}\times
N_{2})}J_{jk}^{\rm(N_{2}\times N_{3})}=\sum_{j=1}^{N_{2}}1=N_{2}$ for all
$1\leq i\leq N_{1}$ and $1\leq k\leq N_{3}$. Therefore we have

${\bf J}^{\rm(G_{1}G_{2})}{\bf J}^{\rm(G_{2}G_{3})}=N^{\rm(G_{2})}{\bf
J}^{\rm(G_{1}G_{3})},$ (32)

where ${\rm G_{1},G_{2},G_{3}\in\\{P,A\\}}$. Using this result, one can see
that the rescaled block matrices of $1$’s satisfy

$\displaystyle\mathbb{\widetilde{J}}^{(0)}\mathbb{\widetilde{J}}^{(0)}$
$\displaystyle=\mathbb{\widetilde{J}}^{\rm(0)},$
$\displaystyle\mathbb{\widetilde{J}}^{(1)}\mathbb{\widetilde{J}}^{(1)}$
$\displaystyle=\mathbb{\widetilde{J}}^{\rm(0)},$
$\displaystyle\mathbb{\widetilde{J}}^{(0)}\mathbb{\widetilde{J}}^{(1)}$
$\displaystyle=\mathbb{\widetilde{J}}^{(1)}\mathbb{\widetilde{J}}^{(0)}=\mathbb{\widetilde{J}}^{\rm(1)}.$
(33)

The multiplication of ${\bf J}$ with the adjacency matrix ${\bf A}^{\rm(PA)}$
or ${\bf A}^{\rm(AP)}={{\bf A}^{\rm(PA)}}^{\intercal}$ is evaluated as

$\displaystyle{\bf A}^{\rm(G_{1}G_{2})}{\bf J}^{\rm(G_{2}G_{3})}={\bf
K}^{\rm(G_{1})}{\bf J}^{\rm(G_{1}G_{3})},$ $\displaystyle{\bf
J}^{\rm(G_{1}G_{2})}{\bf A}^{\rm(G_{2}G_{3})}={\bf J}^{\rm(G_{1}G_{3})}{\bf
K}^{\rm(G_{3})},$ (34)

where we use for instance that
$\sum_{a^{\prime}}A_{pa^{\prime}}J_{a^{\prime}a}=\sum_{a^{\prime}}A_{pa^{\prime}}=k_{p}$
and $\sum_{p^{\prime}}J_{pp^{\prime}}A_{p^{\prime}a}=k_{a}$. Note that
$K_{pp^{\prime}}=k_{p}\delta_{pp^{\prime}}$ and
$K_{aa^{\prime}}=k_{a}\delta_{aa^{\prime}}$. The block adjacency matrix
defined as

$\mathbb{A}\equiv\begin{pmatrix}{\bf 0}&{\bf A}^{\rm(PA)}\\\ {\bf
A}^{\rm(AP)}&{\bf 0}\end{pmatrix}={\bf A}^{\rm(PA)}\oplus{\bf A}^{\rm(AP)},$
(35)

the block rescaled matrices of $1$’s $\mathbb{\widetilde{J}}^{(0)}$ and
$\mathbb{\widetilde{J}}^{(1)}$, and the block degree matrix defined as

$\mathbb{K}\equiv\begin{pmatrix}{\bf K}^{\rm(PP)}&{\bf 0}\\\ {\bf 0}&{\bf
K}^{\rm(AA)}\end{pmatrix}={\bf K}^{\rm(PP)}\oplus{\bf K}^{\rm(AA)}$ (36)

satisfy the following equalities

$\displaystyle\mathbb{A}\mathbb{\widetilde{J}}^{(0)}$
$\displaystyle=\mathbb{K}\mathbb{\widetilde{N}}\mathbb{\widetilde{J}}^{(1)},\
\mathbb{\widetilde{J}}^{(0)}\mathbb{A}=\mathbb{\widetilde{J}}^{(1)}\mathbb{\widetilde{N}}\mathbb{K},$
$\displaystyle\mathbb{A}\mathbb{\widetilde{J}}^{(1)}$
$\displaystyle=\mathbb{K}\mathbb{\widetilde{N}}\mathbb{\widetilde{J}}^{(0)},\
\mathbb{\widetilde{J}}^{(1)}\mathbb{A}=\mathbb{\widetilde{J}}^{(0)}\mathbb{\widetilde{N}}\mathbb{K}$
(37)

with $\mathbb{\widetilde{N}}\equiv\begin{pmatrix}\sqrt{N^{\rm(P)}\over
N^{\rm(A)}}{\bf I}^{\rm(PP)}&{\bf 0}\\\ {\bf 0}&\sqrt{N^{\rm(A)}\over
N^{\rm(P)}}{\bf I}^{\rm(AA)}\end{pmatrix}=\sqrt{N^{\rm(P)}\over
N^{\rm(A)}}{\bf I}^{\rm(PP)}\oplus\sqrt{N^{\rm(A)}\over N^{\rm(P)}}{\bf
I}^{\rm(AA)}$.

Multiplying the degree matrices ${\bf K}^{\rm(PP)}$ and ${\bf K}^{\rm(AA)}$ by
${\bf J}$ matrices yields

$\displaystyle{\bf J}^{\rm(G_{1}G_{2})}{\bf K}^{\rm(G_{2})}{\bf
J}^{\rm(G_{2}G_{3})}=L{\bf J}^{\rm(G_{1}G_{3})},$ $\displaystyle{\bf
J}^{\rm(G_{1}G_{2})}\left({\bf K}^{\rm(G_{2})}\right)^{2}{\bf
J}^{\rm(G_{2}G_{3})}=N^{\rm(G_{2})}\langle k^{2}\rangle^{\rm(G_{2})}{\bf
J}^{\rm(G_{1}G_{3})},$ (38)

where we used
$\sum_{a_{1}a_{2}}J_{aa_{1}}K_{a_{1}a_{2}}J_{a_{2}p}=\sum_{a_{1}}k_{a_{1}}=L$
and
$\sum_{p_{1},p_{2},p_{3}}J_{pp_{1}}K_{p_{1}p_{2}}K_{p_{2}p_{3}}J_{p_{3}a}=\sum_{p_{1}}k_{p_{1}}^{2}=N^{\rm(P)}\langle
k^{2}\rangle^{\rm(P)}$, and $\langle
k^{2}\rangle^{\rm(G)}=\sum_{i\in\mathcal{S}^{\rm(G)}}k_{i}^{2}/N^{\rm(G)}$ is
the mean of the square of the degree of species of group G with
$\mathcal{S}^{\rm(G)}$ the set of group-G species. The block matrices satisfy

$\displaystyle\mathbb{\widetilde{J}}^{(0)}\mathbb{K}\mathbb{\widetilde{J}}^{(0)}$
$\displaystyle=\mathbb{\widetilde{J}}^{(0)}\langle\mathbb{K}\rangle,$
$\displaystyle\mathbb{\widetilde{J}}^{(0)}\mathbb{K}\mathbb{\widetilde{J}}^{(1)}$
$\displaystyle=\langle\mathbb{K}\rangle\mathbb{\widetilde{J}}^{(1)},$
$\displaystyle\mathbb{\widetilde{J}}^{(1)}\mathbb{K}\mathbb{\widetilde{J}}^{(0)}$
$\displaystyle=\mathbb{\widetilde{J}}^{(1)}\langle\mathbb{K}\rangle,$
$\displaystyle\mathbb{\widetilde{J}}^{(1)}\mathbb{K}\mathbb{\widetilde{J}}^{(1)}$
$\displaystyle=\langle k\rangle^{\rm(P)}\langle
k\rangle^{\rm(A)}\langle\mathbb{K}\rangle^{-1}\mathbb{\widetilde{J}}^{(0)},$
$\displaystyle\mathbb{\widetilde{J}}^{(1)}\mathbb{K}^{2}\mathbb{\widetilde{J}}^{(1)}$
$\displaystyle=\langle k^{2}\rangle^{\rm(P)}\langle
k^{2}\rangle^{\rm(A)}\langle\mathbb{K}^{2}\rangle^{-1}\mathbb{\widetilde{J}}^{(0)},$
(39)

with $\langle\mathbb{K}\rangle\equiv\langle k\rangle^{\rm(P)}{\bf
I}^{\rm(PP)}+\langle k\rangle^{\rm(A)}{\bf I}^{\rm(AA)}$ being the sum of the
identity matrices multiplied by the group averages. These relations are valid
also for a function $f(\mathbb{K})$ of $\mathbb{K}$ as

$\displaystyle\mathbb{\widetilde{J}}^{(0)}f(\mathbb{K})\mathbb{\widetilde{J}}^{(0)}$
$\displaystyle=\mathbb{\widetilde{J}}^{(0)}\langle f(\mathbb{K})\rangle,$
$\displaystyle\mathbb{\widetilde{J}}^{(0)}f(\mathbb{K})\mathbb{\widetilde{J}}^{(1)}$
$\displaystyle=\langle f(\mathbb{K})\rangle\mathbb{\widetilde{J}}^{(1)},$
$\displaystyle\mathbb{\widetilde{J}}^{(1)}f(\mathbb{K})\mathbb{\widetilde{J}}^{(0)}$
$\displaystyle=\mathbb{\widetilde{J}}^{(1)}\langle f(\mathbb{K})\rangle,$
$\displaystyle\mathbb{\widetilde{J}}^{(1)}f(\mathbb{K})\mathbb{\widetilde{J}}^{(1)}$
$\displaystyle=\langle f(\mathbb{K})\rangle^{\rm(P)}\langle
f(\mathbb{K})\rangle^{\rm(A)}\langle
f(\mathbb{K})\rangle^{-1}\mathbb{\widetilde{J}}^{(0)},$ (40)

where

$\displaystyle\langle f(\mathbb{K})\rangle$ $\displaystyle\equiv{{\rm
Tr}\,f({\bf K}^{\rm(PP)})\over N^{\rm(P)}}{\bf I}^{\rm(PP)}\oplus{{\rm
Tr}\,f({\bf K}^{\rm(AA)})\over N^{\rm(A)}}{\bf I}^{\rm(AA)}$
$\displaystyle=\langle f(\mathbb{K})\rangle^{\rm(P)}{\bf
I}^{\rm(PP)}\oplus\langle f(\mathbb{K})\rangle^{\rm(A)}{\bf I}^{\rm(AA)},$
(41)

is the sum of the group averages of $f(\mathbb{K})$, and its inverse means
$\langle f(\mathbb{K})\rangle^{-1}={1\over\langle
f(\mathbb{K})\rangle^{\rm(P)}}{\bf I}^{\rm(PP)}\oplus{1\over\langle
f(\mathbb{K})\rangle^{\rm(A)}}{\bf I}^{\rm(AA)}$. For general $z$, $\langle
f(\mathbb{K})\rangle^{z}=\left(\langle f\rangle^{\rm(P)}\right)^{z}{\bf
I}^{\rm(PP)}\oplus\left(\langle f\rangle^{\rm(A)}\right)^{z}{\bf
I}^{\rm(AA)}$. The multiplication of $\langle f(\mathbb{K})\rangle$ and
$\mathbb{\widetilde{J}}^{(1)}$ is not commutative:

$\displaystyle\langle f(\mathbb{K})\rangle\mathbb{\widetilde{J}}^{(1)}$
$\displaystyle=\langle f(\mathbb{K})\rangle^{\rm(P)}\langle
f(\mathbb{K})\rangle^{\rm(A)}\mathbb{\widetilde{J}}^{(1)}\langle
f(\mathbb{K})\rangle^{-1},$ $\displaystyle\mathbb{\widetilde{J}}^{(1)}\langle
f(\mathbb{K})\rangle$ $\displaystyle=\langle
f(\mathbb{K})\rangle^{\rm(P)}\langle f(\mathbb{K})\rangle^{\rm(A)}\langle
f(\mathbb{K})\rangle^{-1}\mathbb{\widetilde{J}}^{(1)},$ (42)

which can be seen by considering for instance the P-block of $\langle
f(\mathbb{K})\rangle\mathbb{\widetilde{J}}^{(1)}\mathbb{X}$ as $\langle
f\rangle^{\rm(P)}{\bf J}^{\rm(PA)}{\bf X}^{\rm(A)}=\langle
f\rangle^{\rm(P)}\langle f\rangle^{\rm(A)}{\bf J}^{\rm(PA)}{1\over\langle
f\rangle^{\rm(A)}}{\bf X}^{\rm(A)}$ with $\mathbb{X}={\bf
X}^{\rm(P)}\oplus{\bf X}^{\rm(A)}$.

##  Appendix B Eigenvalues of the Jacobian matrix

To help understand Eq. (7), here we construct the Jacobian matrix for a small
community as a concrete example and present the eigenvalue perturbation theory
that we used in the main text.

###  B.1 Jacobian matrix for a small community

Let us consider a community consisting of $S=4$ species and a fixed point
$\vec{x}^{*}=(x_{1}^{*}\neq 0,0,0,x_{4}^{*}\neq 0)$ of the LV equation for the
community, which corresponds to $\mathcal{S}^{(0)}=\\{2,3\\}$ and
$\mathcal{S}^{(+)}=\\{1,4\\}$. As shown in Eq. (4), the non-zero components
(abundances) of species $1$ and $4$ in $\mathcal{S}^{(+)}$ satisfy

$\mathbb{B}^{(+)}\begin{pmatrix}x_{1}^{*}\\\
x_{4}^{*}\end{pmatrix}=\begin{pmatrix}B_{11}&B_{14}\\\
B_{41}&B_{44}\end{pmatrix}\begin{pmatrix}x_{1}^{*}\\\
x_{4}^{*}\end{pmatrix}=-\begin{pmatrix}1\\\ 1\end{pmatrix},$

where we introduced the effective interaction matrix $\mathbb{B}^{(+)}$ by
eliminating the rows and columns of the species, $2$ and $3$, of
$\mathcal{S}^{(0)}$ in the original $4\times 4$ interaction matrix
$\mathbb{B}$ while keeping the original indices of rows and columns. Using Eq.
(5) and rearranging the order of indices as $(2\ 3\ 1\ 4)$, we find the
Jacobian matrix at the considered fixed point given by

$\displaystyle\mathbb{H}=$
$\displaystyle\begin{pmatrix}1+B_{21}x_{1}^{*}+B_{24}x_{4}^{*}&0&0&0\\\
0&1+B_{31}x_{1}^{*}+B_{34}x_{4}^{*}&0&0\\\
x_{1}^{*}B_{12}&x_{1}^{*}B_{13}&x_{1}^{*}B_{11}&x_{1}^{*}B_{14}\\\
x_{4}^{*}B_{42}&x_{4}^{*}B_{43}&x_{4}^{*}B_{41}&x_{4}^{*}B_{44}\end{pmatrix},$

which is represented as in Eq. (6) in terms of the block matrices

$\displaystyle{\bf H}^{(00)}$ $\displaystyle=\begin{pmatrix}H_{22}&0\\\
0&H_{33}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}1+B_{21}x_{1}^{*}+B_{24}x_{4}^{*}&0\\\
0&1+B_{31}x_{1}^{*}+B_{34}x_{4}^{*}\end{pmatrix}$

and

${\bf H}^{(++)}=\begin{pmatrix}H_{11}&H_{14}\\\
H_{41}&H_{44}\end{pmatrix}=\begin{pmatrix}x_{1}^{*}B_{11}&x_{1}^{*}B_{14}\\\
x_{4}^{*}B_{41}&x_{4}^{*}B_{44}\end{pmatrix}.$

Note that $B^{(+)}_{ij}=B_{ij}$ if $i,j$ = $1$ or $4$ in our example.

The diagonal elements $H_{22}$ and $H_{33}$ are the eigenvalues of ${\bf
H}^{(00)}$, as ${\bf H}^{(00)}$ is a diagonal matrix. The diagonal elements
can be further simplified. Let us first consider
$H_{22}=1+B_{21}x_{1}^{*}+B_{24}x_{4}^{*}$. A clue is obtained by considering
a neighboring fixed point
$\vec{x}^{*\prime}=(x_{1}^{*\prime},x_{2}^{*\prime},0,x_{4}^{*\prime})$ with
$\mathcal{S}^{(0)\prime}=\\{3\\}$ and $\mathcal{S}^{(+)\prime}=\\{1,2,4\\}$,
where species $2$ has a non-zero component, as well as species $1$ and $4$.
Their abundances satisfy

$\begin{pmatrix}B_{11}&B_{12}&B_{14}\\\ B_{21}&B_{22}&B_{24}\\\
B_{41}&B_{42}&B_{44}\end{pmatrix}\begin{pmatrix}x_{1}^{*\prime}\\\
x_{2}^{*\prime}\\\ x_{4}^{*\prime}\end{pmatrix}=-\begin{pmatrix}1\\\ 1\\\
1\end{pmatrix}.$

Among the three equalities from this equation is
$B_{21}x_{1}^{*\prime}+B_{22}x_{2}^{*\prime}+B_{24}x_{4}^{*\prime}=-1$.
Rearranging the terms and recalling $B_{22}=-1$, we find that
$x_{2}^{*\prime}=1+B_{21}x_{1}^{*\prime}+B_{24}x_{4}^{*\prime}$. If we assume
that $x_{1}^{*\prime}\simeq x_{1}^{*}$ and $x_{4}^{*\prime}\simeq x_{4}^{*}$,
i.e., that the non-zero components of species $1$ and $4$ are similar between
the two fixed points $\vec{x}^{*}$ and $\vec{x}^{*\prime}$, then we find that

$H_{22}=1+B_{21}x_{1}^{*}+B_{24}x_{4}^{*}\simeq
1+B_{21}x_{1}^{*\prime}+B_{24}x_{4}^{*\prime}=x_{2}^{*\prime}.$

The assumption is expected to be valid if $S$ is large and the number of
species having non-zero components is sufficiently large at the fixed point,
for which allowing one more species to have a non-zero component would not
change much the non-zero components of other species. Similarly, one can
simplify $H_{33}$ as $H_{33}\simeq x_{3}^{*\prime}$ by using another
neighboring fixed point corresponding to $\mathcal{S}^{(0)\prime}=\\{2\\}$ and
$\mathcal{S}^{(+)\prime}=\\{1,3,4\\}$. In general, one can obtain the
approximate expression $H_{ii}\simeq x_{i}^{*\prime}$ for every
$i\in\mathcal{S}^{(0)}$ by considering the neighboring fixed point
corresponding to $\mathcal{S}^{(0)\prime}=\mathcal{S}^{(0)}-\\{i\\}$ and
$\mathcal{S}^{(+)\prime}=\mathcal{S}^{(+)}\cup\\{i\\}$. Therefore, given a
fixed point corresponding to $\mathcal{S}^{(0)}$ and $\mathcal{S}^{(+)}$,
$x_{i}^{*\prime}$ for $i\in\mathcal{S}^{(0)}$ means the abundance that the
species $i$ would have if it had a non-zero component like the species of
$\mathcal{S}^{(+)}$.

###  B.2 Eigenvalues of $H_{ij}^{(++)}$ and its largest one

To obtain the eigenvalues $\lambda_{i}^{(+)}$’s of the $S^{(+)}\times S^{(+)}$
matrix $H_{ij}^{(++)}\equiv x_{i}^{*}B^{(+)}_{ij}$ at $x_{i}^{*}$, we
decompose $H^{(++)}_{ij}$ as $H^{(++)}_{ij}=-x_{i}^{*}\delta_{ij}+V_{ij}$ with
$V_{ij}\equiv-cx_{i}^{*}(1-\delta_{ij})+mx_{i}^{*}A_{ij}$. Considering the
eigenvalue expansion $\lambda_{i}^{(+)}\simeq-x_{i}^{*}+V_{ii}+\sum_{j\neq
i}{x_{j}^{*}V_{ji}V_{ij}x_{i}^{*}\over x_{i}^{*}-x_{j}^{*}}\simeq-x_{i}^{*}$
and noting that $V_{ii}=0$, one finds that the eigenvalues can be approximated
by the zeroth-order term as $\lambda_{i}^{(+)}\simeq-x_{i}^{*}$ when $V_{ij}$
is sufficiently small.

It should be also noted that $H_{ij}^{(++)}$ has $\vec{x}^{*}=(x_{i}^{*})$ as
an eigenvector with eigenvalue $-1$;
$\sum_{j}H^{(++)}_{ij}x_{j}^{*}=\sum_{j}x_{i}^{*}B^{(+)}_{ij}\left(-\sum_{\ell}((B^{(+)})^{-1})_{j\ell}\right)=-x_{i}^{*}$.
Therefore the largest real part of the eigenvalues $\lambda^{(+)}_{i}$’s is
approximated by $\max\left(-x_{i}^{*},-1\right)$.

##  Appendix C $\mathbb{B}^{-1}$ to the first order in $m$

For $\mathbb{B}=\mathbb{B}_{0}+m\mathbb{A}$, one can expand its inverse
$\mathbb{B}^{-1}$ in terms of the mutualism strength $m$ as

$\displaystyle\mathbb{B}^{-1}$
$\displaystyle=\mathbb{B}_{0}^{-1}\sum_{n=0}^{\infty}\left(-m\mathbb{A}\mathbb{B}_{0}^{-1}\right)^{n}$
$\displaystyle=\mathbb{B}_{0}^{-1}-m\mathbb{B}_{0}^{-1}\mathbb{A}\mathbb{B}_{0}^{-1}+O(m^{2}).$
(43)

Using Eqs. (37), (39), and (10), one can evaluate $\mathbb{B}^{-1}$ up to the
first order of $m$ as

$\displaystyle\mathbb{B}^{-1}\simeq\mathbb{B}_{0}^{-1}-{m\over(1-c)^{2}}\left[-\mathbb{I}+\mathbbm{\tilde{c}}\mathbb{\widetilde{J}}^{(0)}\right]\mathbb{A}\left[-\mathbb{I}+\mathbbm{\tilde{c}}\mathbb{\widetilde{J}}^{(0)}\right]$
$\displaystyle=\mathbb{B}_{0}^{-1}-{m\over(1-c)^{2}}\left[\mathbb{A}-\mathbb{A}\mathbbm{\tilde{c}}\mathbb{\widetilde{J}}^{(0)}-\mathbbm{\tilde{c}}\mathbb{\widetilde{J}}^{(0)}\mathbb{A}+\mathbbm{\tilde{c}}\mathbb{\widetilde{J}}^{(0)}\mathbb{A}\mathbbm{\tilde{c}}\mathbb{\widetilde{J}}^{(0)}\right]$
$\displaystyle=\mathbb{B}_{0}^{-1}-{m\over(1-c)^{2}}\left[\mathbb{A}-\mathbb{K}\mathbb{\widetilde{N}}\mathbb{\widetilde{J}}^{(1)}\mathbbm{\tilde{c}}-\mathbbm{\tilde{c}}\mathbb{\widetilde{J}}^{(1)}\mathbb{\widetilde{N}}\left(\mathbb{K}-\langle\mathbb{K}\rangle\mathbbm{\tilde{c}}\right)\right].$
(44)

Then the abundance of a plant species $p$ is given by

$\displaystyle x_{p}^{*(1)}=-\sum_{j}(B^{-1})_{pj}=x_{0}^{\rm(P)}$
$\displaystyle+{m\over(1-c)^{2}}\left[k_{p}-k_{p}\tilde{c}^{\rm(A)}-\tilde{c}^{\rm(P)}{L\over
N^{\rm(P)}}+\tilde{c}^{\rm(P)}{L\over N^{\rm(P)}}\tilde{c}^{\rm(A)}\right]$
$\displaystyle=x_{0}^{\rm(P)}+{m\over(1-c)^{2}}(1-\tilde{c}^{\rm(A)})\left(k_{p}-\tilde{c}^{\rm(P)}\langle
k\rangle^{\rm(P)}\right)$ $\displaystyle=x_{0}^{\rm(P)}\left[1+{m\over
1-c}{1-\tilde{c}^{\rm(A)}\over
1-\tilde{c}^{\rm(P)}}\left(k_{p}-\tilde{c}^{\rm(P)}\langle
k\rangle^{\rm(P)}\right)\right],$ (45)

and that of an animal species $a$ is

$\displaystyle x_{a}^{*(1)}=x_{0}^{\rm(A)}\left[1+{m\over
1-c}{1-\tilde{c}^{\rm(P)}\over
1-\tilde{c}^{\rm(A)}}\left(k_{a}-\tilde{c}^{\rm(A)}\langle
k\rangle^{\rm(A)}\right)\right].$ (46)

##  Appendix D $\mathbb{B}^{-1}$ under the annealed approximation for the
adjacency matrix

In the annealed approximation, the adjacency matrix element $A_{ij}$ is
approximated by the probability that the two nodes $i$ and $j$ are connected
by a link in the ensemble of networks preserving the given degree sequence
$\\{k_{i}\\}$ as

$\tilde{A}_{ij}={k_{i}k_{j}\over L}$ (47)

with $L$ the total number of links. Equivalently, the block adjacency matrix
takes the form

$\mathbb{\widetilde{A}}={1\over
L}\mathbb{K}\mathbb{J}^{(1)}\mathbb{K}={1\over\sqrt{\langle
k\rangle^{\rm(P)}\langle
k\rangle^{\rm(A)}}}\mathbb{K}\mathbb{\widetilde{J}}^{(1)}\mathbb{K}.$ (48)

Using this, one finds that the terms
$\mathbb{B}_{0}^{-1}(-m\mathbb{\widetilde{A}}\mathbb{B}_{0}^{-1})^{n}$
appearing in Eq. (12) are simplified. Evaluating the first three terms, one
can find the expression for general $n$ by induction. Let us first consider
the term with $n=1$, which is evaluated as

$\displaystyle\mathbb{B}_{0}^{-1}(-m\mathbb{\widetilde{A}})\mathbb{B}_{0}^{-1}$
$\displaystyle={-m\over(1-c)^{2}}{\left(-\mathbb{I}+\mathbbm{\tilde{c}}\mathbb{\widetilde{J}}^{(0)}\right)\mathbb{K}\mathbb{\widetilde{J}}^{(1)}\mathbb{K}\left(-\mathbb{I}+\mathbbm{\tilde{c}}\mathbb{\widetilde{J}}^{(0)}\right)\over\sqrt{\langle
k\rangle^{\rm(P)}\langle k\rangle^{\rm(A)}}}$
$\displaystyle={-m\over(1-c)^{2}}\times$
$\displaystyle{\mathbb{K}\mathbb{\widetilde{J}}^{(1)}\mathbb{K}-\mathbbm{\tilde{c}}\langle\mathbb{K}\rangle\mathbb{\widetilde{J}}^{(1)}\mathbb{K}-\mathbb{K}\mathbb{\widetilde{J}}^{(1)}\mathbbm{\tilde{c}}\langle\mathbb{K}\rangle+\mathbbm{\tilde{c}}\langle\mathbb{K}\rangle\mathbb{\widetilde{J}}^{(1)}\langle\mathbb{K}\rangle\mathbbm{\tilde{c}}\over\sqrt{\langle
k\rangle^{\rm(P)}\langle k\rangle^{\rm(A)}}}$
$\displaystyle={-m\over(1-c)^{2}}{\left(\mathbb{K}-\mathbbm{\tilde{c}}\langle\mathbb{K}\rangle\right)\mathbb{\widetilde{J}}^{(1)}\left(\mathbb{K}-\mathbbm{\tilde{c}}\langle\mathbb{K}\rangle\right)\over\sqrt{\langle
k\rangle^{\rm(P)}\langle k\rangle^{\rm(A)}}}$ $\displaystyle=-{\tilde{m}\over
1-c}\mathbb{\widetilde{K}}\mathbb{\widetilde{J}}^{(1)}\mathbb{\widetilde{K}},$
(49)

where Eq. (40) is used and the rescaled degree matrix is introduced,

$\mathbb{\widetilde{K}}\equiv{\mathbb{K}-\mathbbm{\tilde{c}}\langle\mathbb{K}\rangle\over\sqrt{\langle\mathbb{K}(\mathbb{K}-\mathbbm{\tilde{c}}\langle\mathbb{K}\rangle)\rangle}}={\mathbb{K}-\mathbbm{\tilde{c}}\langle\mathbb{K}\rangle\over\langle\mathbb{\widetilde{K}}\mathbb{K}\rangle},$
(50)

and the rescaled mutualism strength $\tilde{m}$, defined in Eq. (17) in the
main text, is here also represented as

$\displaystyle\tilde{m}$ $\displaystyle={m\over
1-c}{\langle\mathbb{K}\mathbb{\widetilde{K}}\rangle^{\rm(P)}\langle\mathbb{K}\mathbb{\widetilde{K}}\rangle^{\rm(A)}\over\sqrt{\langle\mathbb{K}\rangle^{\rm(P)}\langle\mathbb{K}\rangle^{\rm(A)}}}.$
(51)

Notice that $\langle k\rangle^{\rm(P)}=\langle\mathbb{K}\rangle^{\rm(P)}$ and
$\langle\mathbb{K}\mathbb{\widetilde{K}}\rangle^{\rm(P)}=\langle
k^{2}\rangle^{\rm(P)}-\tilde{c}^{\rm(P)}(\langle
k\rangle^{\rm(P)})^{2}=(\langle
k\rangle^{\rm(P)})^{2}\left(\xi^{\rm(P)}-\tilde{c}^{\rm(P)}\right)$ from Eq.
(41).

The term with $n=2$ is evaluated as

$\displaystyle\mathbb{B}_{0}^{-1}(-m\mathbb{\widetilde{A}})\mathbb{B}_{0}^{-1}(-m\mathbb{\widetilde{A}})\mathbb{B}_{0}^{-1}$
$\displaystyle={-\tilde{m}\over
1-c}\mathbb{\widetilde{K}}\mathbb{\widetilde{J}}^{(1)}\mathbb{\widetilde{K}}{-m\over\sqrt{\langle
k\rangle^{\rm(P)}\langle
k\rangle^{\rm(A)}}}\mathbb{K}\mathbb{\widetilde{J}}^{(1)}\mathbb{K}{-\mathbb{I}+\mathbbm{\tilde{c}}\mathbb{\widetilde{J}}^{(0)}\over
1-c}$ $\displaystyle={\tilde{m}\over
1-c}\mathbb{\widetilde{K}}{m\langle\mathbb{\widetilde{K}}\mathbb{K}\rangle^{\rm(P)}\langle\mathbb{\widetilde{K}}\mathbb{K}\rangle^{\rm(A)}\langle\mathbb{\widetilde{K}}\mathbb{K}\rangle^{-1}\mathbb{\widetilde{J}}^{(0)}\over\sqrt{\langle
k\rangle^{\rm(P)}\langle
k\rangle^{\rm(A)}}}{-\mathbb{K}+\mathbbm{\tilde{c}}\langle\mathbb{K}\rangle\over
1-c}$ $\displaystyle=-{\tilde{m}^{2}\over
1-c}\mathbb{\widetilde{K}}\mathbb{\widetilde{J}}^{(0)}\mathbb{\widetilde{K}}.$
(52)

The term with $n=3$ is

$\displaystyle\mathbb{B}_{0}^{-1}(-m\mathbb{\widetilde{A}})\mathbb{B}_{0}^{-1}(-m\mathbb{\widetilde{A}})\mathbb{B}_{0}^{-1}(-m\mathbb{\widetilde{A}})\mathbb{B}_{0}^{-1}$
$\displaystyle={-\tilde{m}^{2}\over
1-c}\mathbb{\widetilde{K}}\mathbb{\widetilde{J}}^{(0)}\mathbb{\widetilde{K}}{-m\over\sqrt{\langle
k\rangle^{\rm(P)}\langle
k\rangle^{\rm(A)}}}\mathbb{K}\mathbb{\widetilde{J}}^{(1)}\mathbb{K}{-\mathbb{I}+\mathbbm{\tilde{c}}\mathbb{\widetilde{J}}^{(0)}\over
1-c}$ $\displaystyle={\tilde{m}^{2}\over
1-c}\mathbb{\widetilde{K}}{m\langle\mathbb{\widetilde{K}}\mathbb{K}\rangle\mathbb{\widetilde{J}}^{(1)}\over\sqrt{\langle
k\rangle^{\rm(P)}\langle
k\rangle^{\rm(A)}}}{-\mathbb{K}+\mathbbm{\tilde{c}}\langle\mathbb{K}\rangle\over
1-c}$ $\displaystyle=-{\tilde{m}^{2}\over
1-c}\mathbb{\widetilde{K}}{m\langle\mathbb{\widetilde{K}}\mathbb{K}\rangle^{\rm(P)}\langle\mathbb{\widetilde{K}}\mathbb{K}\rangle^{\rm(A)}\over(1-c)\sqrt{\langle
k\rangle^{\rm(P)}\langle
k\rangle^{\rm(A)}}}\mathbb{\widetilde{J}}^{(1)}\langle\mathbb{K}\mathbb{\widetilde{K}}\rangle^{-1}\left(\mathbb{K}-\mathbbm{\tilde{c}}\langle\mathbb{K}\rangle\right)$
$\displaystyle=-{\tilde{m}^{3}\over
1-c}\mathbb{\widetilde{K}}\mathbb{\widetilde{J}}^{(1)}\mathbb{\widetilde{K}}.$
(53)

From Eqs. (49), (52), and (53), one can obtain by induction Eq. (16) and we
find that the inverse of the interaction matrix is evaluated as

$\displaystyle\mathbb{\widetilde{B}}^{-1}=\mathbb{B}_{0}^{-1}-$
$\displaystyle{1\over
1-c}\mathbb{\widetilde{K}}\left(\mathbb{\widetilde{J}}^{(1)}\sum_{n=1,3,5,\ldots}\tilde{m}^{n}+\mathbb{\widetilde{J}}^{(0)}\sum_{n=2,4,6,\ldots}\tilde{m}^{n}\right)\mathbb{\widetilde{K}}$
$\displaystyle=\mathbb{B}_{0}^{-1}-{1\over 1-c}{\tilde{m}\over
1-\tilde{m}^{2}}\mathbb{\widetilde{K}}\left(\tilde{m}\mathbb{\widetilde{J}}^{(0)}+\mathbb{\widetilde{J}}^{(1)}\right)\mathbb{\widetilde{K}},$
(54)

which is given also in Eq. (19).

##  Appendix E Different measures of the species abundance in the long-time
limit

Figure 4: Species abundance in the full coexistence phase. Stationary-state
abundances $x_{i}^{\rm(st)}$ of individual species $i$ are compared with the
stable fixed-point ones $x_{i}^{*}$’s and $\tilde{x}_{i}^{*}$’s for $c=0.1$
and $\tilde{m}=0.4$. The species index $i$ is arranged in the descending order
of degree among plants and among animals. Inset: Abundance vs. degree.

In the present study appear a couple of different measures of the species
abundance, being numerical or analytical solutions to the LV equations. We
summarize their notations here to help distinguish and understand them.

  * •

$x_{i}^{\rm(st)}$ : It is defined in Eq. (3) and represents the numerical
solution $x_{i}(T)$ at the final step to the LV equation in Eq. (1) with using
the interaction matrix $\mathbb{B}$ constructed from the original adjacency
matrix $\mathbb{A}$.

  * •

$\tilde{x}_{i}^{\rm(st)}$ : It represents the numerical solution
$\tilde{x}_{i}(T)$ at the final step to the LV equation in Eq. (1) with using
the interaction matrix $\mathbb{\widetilde{B}}$ constructed from the
factorized adjacency matrix $\mathbb{\widetilde{A}}$ in Eq. (15) under the
annealed approximation.

  * •

$x_{i}^{*}$ : It is defined in Eq. (4) and represents the stable fixed point
to the LV equation in Eq. (1) with using $\mathbb{B}$. The sets of surviving
and extinct species, $\mathcal{S}^{(+)}$ and $\mathcal{S}^{(0)}$ can be
obtained by updating iteratively $\mathbb{B}^{(+)}$ and $x_{i}^{*}$, using Eq.
(4), as described in Sec. V.2.

  * •

$\tilde{x}_{i}^{*}$ : It represents the stable fixed point to the LV equation
in Eq. (1) with using $\mathbb{\widetilde{B}}$ under the annealed
approximation. While it can be evaluated numerically by Eq. (4) with using
$\mathbb{\widetilde{B}}$, its analytic expression is available in Eq. (28).
The sets of surviving and extinct species, $\mathcal{S}^{(+)}$ and
$\mathcal{S}^{(0)}$, can be obtained by updating iteratively
$\mathbb{\widetilde{B}}^{(+)}$ and $\tilde{x}_{i}^{*}$, using Eq. (28), as
described in Sec. V.2.

In Fig. 4, we present $x_{i}^{\rm(st)}$, $x_{i}^{*}$ and $\tilde{x}_{i}^{*}$
for $c=0.1$ and $\tilde{m}=0.4$ in the full coexistence phase.

##  Appendix F Effective quantities in Eq. (28)

Figure 5: Relative variations of the network structural properties with the
rescaled mutualism strength $\tilde{m}$ for $c=0.1$. Shown are the relative
variations of the mean degree $\langle k\rangle$ and the degree heterogeneity
minus the rescaled competition $\xi-\tilde{c}$ for plants and animals. The
relative variation is evaluated e.g., as ${\langle k\rangle^{\rm(P,+)}-\langle
k\rangle^{\rm(P)}\over\langle k\rangle^{\rm(P)}}$.

The effective interaction matrix $\mathbb{\widetilde{B}}^{(+)}$ for the
surviving species is obtained by removing the rows and columns corresponding
to the species belonging to $\mathcal{S}^{(0)}$ in the full interaction matrix
$\mathbb{\widetilde{B}}$. For later use, we introduce $\mathcal{S}^{\rm(P,+)}$
and $\mathcal{S}^{\rm(A,+)}$ to denote the set of plant and animal species,
respectively, in $\mathcal{S}^{(+)}$ to be assigned non-zero components, and
$\mathcal{S}^{\rm(P,0)}$ and $\mathcal{S}^{\rm(A,0)}$ to denote the set of
plant and animal species, respectively, in $\mathcal{S}^{(0)}$ to be assigned
zero components. Their sizes are $N^{\rm(P,+)}$, $N^{\rm(A,+)}$,
$N^{\rm(P,0)}$, and $N^{\rm(A,0)}$. The effective interaction matrix
$\mathbb{\widetilde{B}}^{(+)}$ is of size $S^{(+)}\times S^{(+)}$ with
$S^{(+)}=N^{\rm(P,+)}+N^{\rm(A,+)}$ and takes the form

$\mathbb{\widetilde{B}}^{(+)}=-\mathbb{I}^{(+)}+c(\mathbb{J}^{\rm(0,+)}-\mathbb{I}^{(+)})+m\mathbb{\widetilde{A}}^{(+)},$
(55)

where $\mathbb{\widetilde{A}}^{(+)}$ is the effective adjacency matrix for the
plant and animal species in $\mathcal{S}^{(+)}$.

The effective adjacency matrix
$\mathbb{\widetilde{A}}^{(+)}={\bf\tilde{A}}^{\rm(PA,+)}\oplus{\bf\tilde{A}}^{\rm(AP,+)}$
is obtained by removing in $\mathbb{\widetilde{A}}$ the rows and columns of
the species in $\mathcal{S}^{(0)}$. Therefore it holds that
$\tilde{A}^{(+)}_{pa}=\tilde{A}_{pa}=k_{p}k_{a}/L$ if $p$ and $a$ are in
$\mathcal{S}^{(+)}$. Once $\mathbb{\widetilde{A}}^{(+)}$ is given, the
effective network quantities such as the effective degree can be derived from
$\mathbb{\widetilde{A}}^{(+)}$. Moreover, $\mathbb{\widetilde{A}}^{(+)}$
maintains the factorized form in terms of the effective degrees and the
effective number of links. Therefore the effective quantities can be inserted
into Eq. (20), developed originally with the factorized adjacency matrix, to
yield Eq. (28). Below we present how to evaluate them specifically.

(i) The effective rescaled competition strength is
$\tilde{c}^{\rm(G,+)}\equiv{cN^{\rm(G,+)}\over cN^{\rm(G,+)}+1-c}$ with G
being P or A.

(ii) The effective zeroth-order abundance is
$x_{0}^{\rm(G,+)}\equiv{1-\tilde{c}^{\rm(G,+)}\over 1-c}$.

(iii) The total number of links is
$L^{(+)}\equiv\sum_{p\in\mathcal{S}^{\rm(P,+)},a\in\mathcal{S}^{\rm(A,+)}}\tilde{A}^{(+)}_{pa}={1\over
L}\sum_{p\in\mathcal{S}^{\rm(P,+)}}k_{p}\sum_{a\in\mathcal{S}^{\rm(A,+)}}k_{a}=L\,\ell^{\rm(P,+)}\ell^{\rm(A,+)}$
with $\ell^{\rm(P,+)}\equiv{1\over L}\sum_{p\in\mathcal{S}^{\rm(P,+)}}k_{p}$
and $\ell^{\rm(A,+)}\equiv{1\over L}\sum_{a\in\mathcal{S}^{\rm(A,+)}}k_{a}$
denoting the ratio of the links incident on the plant and animal species of
$\mathcal{S}^{(+)}$, respectively, to the original number of links
$L=\sum_{p}k_{p}=\sum_{a}k_{a}$.

(iv) The effective degree of a plant or animal species in $\mathcal{S}^{(+)}$
is evaluated as
$k_{p}^{(+)}\equiv\sum_{a\in\mathcal{S}^{\rm(A,+)}}\tilde{A}^{(+)}_{pa}=k_{p}\sum_{a\in\mathcal{S}^{\rm(A,+)}}{k_{a}\over
L}=k_{p}\ell^{\rm(A,+)}$ and
$k_{a}^{(+)}\equiv\sum_{p\in\mathcal{S}^{\rm(P,+)}}\tilde{A}^{(+)}_{pa}=k_{a}\ell^{\rm(P,+)}$,
satisfying
$L^{(+)}=\sum_{p\in\mathcal{S}^{\rm(P,+)}}k_{p}^{(+)}=\sum_{a\in\mathcal{S}^{\rm(A,+)}}k_{a}^{(+)}$.

(v) The effective adjacency matrix maintains its factorized form
$\tilde{A}^{(+)}_{pa}=\tilde{A}_{pa}={k_{p}k_{a}\over
L}={k_{p}^{(+)}k_{a}^{(+)}\over L^{(+)}}$ in terms of the effective degrees
and the effective numbers of links defined above.

(vi) The effective degree heterogeneity is evaluated as
$\xi^{\rm(P,+)}={\langle k^{2}\rangle^{\rm(P,+)}\over\langle
k\rangle^{\rm(P,+)}}$ with the moments given by $\langle
k^{n}\rangle^{\rm(P,+)}={1\over
N^{\rm(P,+)}}\sum_{p\in\mathcal{S}^{\rm(P,+)}}\left(k_{p}^{(+)}\right)^{n}$.
$\xi^{\rm(A,+)}$ and $\langle k^{n}\rangle^{\rm(A,+)}$ are evaluated in the
same manner.

(vii) The effective rescaled degree is evaluated by
$\bar{k}_{i}^{(+)}\equiv{{k_{i}^{(+)}\over\langle k\rangle^{({\rm
G}_{i},+)}}-\tilde{c}^{({\rm G}_{i},+)}\over\xi^{({\rm
G}_{i},+)}-\tilde{c}^{({\rm G}_{i},+)}}$.

(viii) The effective asymmetry factor is evaluated by
$\eta^{\rm(PA,+)}\equiv{1-\tilde{c}^{\rm(A,+)}\over
1-\tilde{c}^{\rm(P,+)}}\sqrt{\langle
k\rangle^{\rm(P,+)}(\xi^{\rm(P,+)}-\tilde{c}^{\rm(P,+)})\over\langle
k\rangle^{\rm(A,+)}(\xi^{\rm(A,+)}-\tilde{c}^{\rm(A,+)})}={1\over\eta^{\rm(AP,+)}}$.

(ix) The effective rescaled mutualism strength $\tilde{m}^{(+)}$ is evaluated
by

$\displaystyle\tilde{m}^{(+)}\equiv$ $\displaystyle\frac{m}{1-c}$
$\displaystyle\times$ $\displaystyle\sqrt{\langle k\rangle^{(P,+)}\langle
k\rangle^{(A,+)}(\xi^{(P,+)}-\tilde{c}^{(P,+)})(\xi^{(A,+)}-\tilde{c}^{(A,+)})}$
(56)

exactly in the same manner as Eq. (17) with using the effective quantities.

The relative variations of the effective quantities with respect to their
original values are shown in Fig. 5. While the extinction of small-degree
species may make the effective mean degree $\langle k\rangle^{(+)}$ larger
than the original one $\langle k\rangle$, the hub plants and animals lose
their significant portions of partners, resulting in the reduction of the
effective degree heterogeneity. The rescaled competition $\tilde{c}^{(+)}$
decreases as more species go extinct with increasing $\tilde{m}$. These
quantities together determine the effective rescaled mutualism
$\tilde{m}^{(+)}$, which turns out to be smaller than $\tilde{m}$ as shown in
Fig. 3 (b).

##  Appendix G Species abundance under the annealed adjacency matrix

Figure 6: Species abundance under the annealed adjacency matrix. (a) Time-
evolution of the abundances of individual species (different lines) obtained
by numerically integrating Eq. (1) with the annealed adjacency matrix
$\tilde{A}_{ij}$ used for $c=0.1$ and $m=0.2$. (b) The stationary-state
abundance vs. degree for animal species with $c=0.1$. (c) Stationary-state
abundances $\tilde{x}_{i}^{\rm(st)}$ of individual species $i$ are compared
with the stable fixed point components $\tilde{x}_{i}^{*}$’s for $c=0.1$ and
$\tilde{m}=0.8$ belonging to the selective extinction phase. (d) Fraction
$\tilde{r}_{e}^{\rm(st)}$ of extinct species and $\tilde{r}_{u}$ of the
abundance-diverging species based on the stationary-state abundance
$\tilde{x}_{i}^{\rm(st)}$. They are compared with the theoretical predictions
$\tilde{r}_{e}^{*}$ and $\tilde{r}_{u}^{*}$ based on the stable fixed point
$\tilde{x}_{i}^{*}$’s. (e) The collapse of the minimum abundances
$\tilde{x}_{\rm min}^{\rm(st)}$ rescaled as in Eq. (29) in 46 real-world
communities as functions of $\tilde{m}^{(+)}$.

Here we present the plots of the species abundances, the fraction of extinct
and abundance-diverging species, and the rescaled minimum abundance in case of
the annealed interaction matrix $\mathbb{\widetilde{B}}$ and the annealed
adjacency matrix $\mathbb{\widetilde{A}}$ in Fig. 6.

##  Appendix H Accuracy in the prediction of the extinction of individual
species

Figure 7: Accuracy of the analytic formula of the stable fixed-point
abundances in predicting the extinction or survival of individual species.
Shown in the $(c,\tilde{m})$ plane is the fraction of the species whose
survival or extinction is predicted identically by (a) both the stationary-
state abundance $x_{i}^{\rm(st)}$ under the original interaction matrix
$\mathbb{B}$, and the stable fixed point $\tilde{x}_{i}^{*}$, and (b) both
$\tilde{x}_{i}^{\rm(st)}$ under the annealed interaction matrix
$\mathbb{\widetilde{B}}$, and the stable fixed point $\tilde{x}_{i}^{*}$.

We consider a species extinct if its abundance is smaller than
$\epsilon=10^{-5}$ and surviving otherwise. The criterion is used to assess
the stationary-state abundance $x_{i}^{\rm(st)}$ and $\tilde{x}_{i}^{\rm(st)}$
and discriminate the fate of $i$ evolving under the original and the annealed
interaction matrix, respectively. To illuminate the predictive power of the
stable fixed point abundance $\tilde{x}_{i}^{*}$ for the fate - survival or
extinction - of individual species, we compute the fraction of the species
that are correctly predicted, i.e., found to be surviving in both abundances,
$x_{i}^{\rm(st)}\geq\epsilon$ and $\tilde{x}_{i}^{*}\geq\epsilon$ or found to
be extinct in both, $x_{i}^{\rm(st)}<\epsilon$ and
$\tilde{x}_{i}^{*}<\epsilon$, which we can consider as the accuracy of the
stable fixed point abundances in the prediction of species extinction and
present in Fig. 7 (a). We also do the same analysis with
$\tilde{x}_{i}^{\rm(st)}$ and $\tilde{x}_{i}^{*}$ and show the result in Fig.
7 (b). On the average across parameters in the selective extinction phase, the
accuracy of $\tilde{x}_{i}^{*}$ in predicting extinction/survival amounts to
79.3% for $x_{i}^{\rm(st)}$ under the original interaction matrix and 99.2%
for $\tilde{x}_{i}^{\rm(st)}$ under the annealed interaction matrix.

## References

  * May (1972) R. M. May, Nature 238, 413 (1972).
  * Goh (1979) B. S. Goh, Am. Nat. 113, 261 (1979).
  * Allesina and Tang (2015) S. Allesina and S. Tang, Popul. Ecol. 57, 63 (2015).
  * Stone (2016) L. Stone, Nat. Commun. 7, 1 (2016).
  * Grilli  _et al._ (2017) J. Grilli, M. Adorisio, S. Suweis, G. Barabás, J. R. Banavar, S. Allesina, and A. Maritan, Nat. Commun. 8, 14389 (2017).
  * Bunin (2017) G. Bunin, Phys. Rev. E 95, 042414 (2017).
  * Cui  _et al._ (2020) W. Cui, R. Marsland, and P. Mehta, Phys. Rev. Lett. 125, 048101 (2020).
  * Pettersson  _et al._ (2020a) S. Pettersson, V. M. Savage, and M. N. Jacobi, J. R. Soc. Interface 17, 20190391 (2020a).
  * Pettersson  _et al._ (2020b) S. Pettersson, V. M. Savage, and M. N. Jacobi, Phys. Rev. E 102, 062405 (2020b).
  * Montoya and Solé (2002) J. M. Montoya and R. V. Solé, J. Theor. Biol. 214, 405 (2002).
  * Dunne  _et al._ (2002) J. A. Dunne, R. J. Williams, and N. D. Martinez, Proc. Natl. Acad. Sci. U.S.A. 99, 12917 (2002).
  * Bascompte  _et al._ (2003) J. Bascompte, P. Jordano, C. J. Melián, and J. M. Olesen, Proc. Natl. Acad. Sci. U.S.A. 100, 9383 (2003).
  * Jordano  _et al._ (2003) P. Jordano, J. Bascompte, and J. M. Olesen, Ecol. Lett. 6, 69 (2003).
  * Montoya  _et al._ (2006) J. M. Montoya, S. L. Pimm, and R. V. Solé, Nature 442, 259 (2006).
  * Bascompte and Jordano (2007) J. Bascompte and P. Jordano, Annu. Rev. Ecol. Evol. Syst. 38, 567 (2007).
  * Guimarães  _et al._ (2007) P. R. Guimarães, G. Machado, M. A. de Aguiar, P. Jordano, J. Bascompte, A. Pinheiro, and S. F. dos Reis, J. Theor. Biol. 249, 181 (2007).
  * Olesen  _et al._ (2007) J. M. Olesen, J. Bascompte, Y. L. Dupont, and P. Jordano, Proc. Natl. Acad. Sci. U.S.A. 104, 19891 (2007).
  * Thébault and Fontaine (2010) E. Thébault and C. Fontaine, Science 329, 853 (2010).
  * Maeng and Lee (2011) S. E. Maeng and J. W. Lee, J. Korean Phys. Soc. 58, 851 (2011).
  * Bastolla  _et al._ (2009) U. Bastolla, M. A. Fortuna, A. Pascual-García, A. Ferrera, B. Luque, and J. Bascompte, Nature 458, 1018 (2009).
  * Suweis  _et al._ (2013) S. Suweis, F. Simini, J. R. Banavar, and A. Maritan, Nature 500, 449 (2013).
  * Saavedra  _et al._ (2016) S. Saavedra, R. P. Rohr, J. M. Olesen, and J. Bascompte, Ecol. Evol. 6, 997 (2016).
  * Yan  _et al._ (2017) G. Yan, N. D. Martinez, and Y. Y. Liu, J. Roy. Soc. Interface 14, 20170189 (2017).
  * Maeng  _et al._ (2012) S. E. Maeng, J. W. Lee, and D.-S. Lee, Phys. Rev. Lett. 108, 108701 (2012).
  * Pascual-García and Bastolla (2017) A. Pascual-García and U. Bastolla, Nat. Commun. 8, 14326 (2017).
  * Barbier  _et al._ (2018) M. Barbier, J.-F. Arnoldi, G. Bunin, and M. Loreau, Proc. Natl. Acad. Sci. U.S.A. 115, 2156 (2018).
  * Gracia-Lázaro  _et al._ (2018) C. Gracia-Lázaro, L. Hernández, J. Borge-Holthoefer, and Y. Moreno, Sci. Rep. 8, 9253 (2018).
  * Wang  _et al._ (2021) X. Wang, T. Peron, J. L. A. Dubbeldam, S. Kèfi, and Y. Moreno, arXiv:2102.02259 (2021).
  * Maeng  _et al._ (2019) S. E. Maeng, J. W. Lee, and D.-S. Lee, J. Stat. Mech. 2019, 033502 (2019).
  * Cai  _et al._ (2020) W. Cai, J. Snyder, A. Hastings, and R. M. D’Souza, Nat. Commun. 11, 5470 (2020).
  * James  _et al._ (2012) A. James, J. W. Pitchford, and M. J. Plank, Nature 487, 227 (2012).
  * Saavedra and Stouffer (2013) S. Saavedra and D. B. Stouffer, Nature 500, E1 (2013).
  * Allesina and Tang (2012) S. Allesina and S. Tang, Nature 483, 205 (2012).
  * Stone (2020) L. Stone, Nat. Commun. 11, 2648 (2020).
  * (35) Http://www.web-of-life.es/.
  * Arroyo  _et al._ (1982) M. T. K. Arroyo, R. Primack, and J. Armesto, Am. J. Bot. 69, 82 (1982).
  * Lee  _et al._ (2009) S. H. Lee, M. Ha, H. Jeong, J. D. Noh, and H. Park, Phys. Rev. E 80, 051127 (2009).
  * (38) The matlab code is available at https://github.com/deoksunlee/Stability-and-selective-extinction-in-complex-mutualistic-networks.