# Evaluating The Impact Of Species Specialisation On Ecological Network
Robustness Using Analytic Methods

Chris Jones  School of Mathematics, University of Bristol, Bristol, UK
0000-0003-4698-5573 Damaris Zurell  Institute of Biochemistry and Biology,
University of Potsdam, Potsdam, Germany  0000-0002-4628-3558 Karoline Wiesner
Institute of Physics and Astronomy, University of Potsdam, Potsdam, Germany
0000-0003-2944-1988

###### Abstract

Ecological networks describe the interactions between different species,
informing us of how they rely on one another for food, pollination and
survival. If a species in an ecosystem is under threat of extinction, it can
affect other species in the system and possibly result in their secondary
extinction as well. Understanding how (primary) extinctions cause secondary
extinctions on ecological networks has been considered previously using
computational methods. However, these methods do not provide an explanation
for the properties which make ecological networks robust, and can be
computationally expensive. We develop a new analytic model for predicting
secondary extinctions which requires no non-deterministic computational
simulation. Our model can predict secondary extinctions when primary
extinctions occur at random or due to some targeting based on the number of
links per species or risk of extinction, and can be applied to an ecological
network of any number of layers. Using our model, we consider how false
positives and negatives in network data affect predictions for network
robustness. We have also extended the model to predict scenarios in which
secondary extinctions occur once species lose a certain percentage of
interaction strength, and to model the loss of interactions as opposed to just
species extinction. From our model, it is possible to derive new analytic
results such as how ecological networks are most robust when secondary species
degree variance is minimised. Additionally, we show that both specialisation
and generalisation in distribution of interaction strength can be advantageous
for network robustness, depending upon the extinction scenario being
considered.

##  1 Introduction

No species exists in isolation, depending upon interactions with other species
to feed, reproduce or maintain a stable population [1, 2]. Modelling the
interactions between species is therefore of great importance in ecology, and
one approach to this problem is to model interactions as an ecological network
[3, 4]. Ecosystems are increasingly threatened by the effects of climate
change [5, 6] which can cause sudden and widespread extinction events.
Therefore, it is useful to model extinctions on ecological networks in order
to understand the possible knock-on effects of species extinctions, as this
may help to identify methods for conserving or reinforcing ecosystems in the
future [7, 8, 9].

Species extinctions on ecological networks have been extensively studied in
the past 20 years, with simplistic topological models providing predictions
for the impact of extinctions under scenarios including: extinctions which
occur at random or with some ordering [10], extinctions on networks made up of
numerous trophic levels [11], and extinctions which occur due to a loss of
interaction strength over a certain threshold [12].

The models used are not the only possible approach to understanding the
robustness of ecological networks. Other models consider the size of the
largest component in the interaction network [13, 14], or take a more
dynamical approach as is the case with Bayesian network models [15, 16]. Here
we restrict ourselves to what we refer to as simplistic topological network
models, which originate from Memmott et al. [10]. In these models, we are
concerned with the point at which a given species goes extinct due to losing
either a certain number of neighbours or a certain amount of interaction
strength.

Previous work on simplistic topological network models have been largely
computational, where extinctions are simulated in order to obtain predictions.
Limited analytic work has been done to predict the robustness of ecological
networks which are either maximally or minimally nested [17], but there is no
existing analytic framework which can predict the robustness of any given
simple ecological network. In the following, we develop such a model, which
improves upon computational methods by providing an insight into the
properties that make ecological networks robust, and by cutting computational
cost.

We start by considering the same scenario put forward by Memmott et al. [10],
where a bipartite mutualistic network (such as a plant pollinator network)
undergoes extinctions on one trophic level, with species on the other trophic
level experiencing secondary extinctions if they lose some or all of their
neighbours. Secondary extinctions may be predicted for random or targeted
primary extinctions, and secondary extinctions are predictable on networks
with more than two trophic layers [11]. Our model may also be used to predict
the effects of errors in network data, where interactions are erroneously
included or excluded. The model is then developed further, taking into account
the variable interaction strengths of neighbouring species, where a species
will go extinct if it loses a certain amount of interaction strength, as
considered by Schleuning et al. [12]. We also consider the scenario in which
species go extinct gradually, modelled by the loss of interaction strength as
opposed to entire species.

Having developed an analytic model for these scenarios, we can determine the
topological properties which make ecological networks robust. Previously, the
roles of nestedness and specialisation have been debated as possible sources
of robustness [10, 17, 18]. We use our model to demonstrate that, when
interaction strength is irrelevant, a network is most robust against random
extinctions when the variance of its secondary species degree distribution is
minimised. When interaction strength is included, we show that if secondary
species’ interaction strength is maximally specialised then network robustness
is constant regardless of network degree distribution or extinction
sensitivity. If interaction strength is maximally generalised, networks with
high degree secondary species have robustness that is solely dependent upon
extinction sensitivity. As a result, high specialisation makes a network more
robust if it is highly sensitive to interaction loss, and high generalisation
is better for robustness if interaction loss sensitivity is low.

##  2 Introducing the Analytic Framework

In the model of Memmott et al. [10], species in one trophic level (e.g.
pollinators) undergo extinctions, and this impacts species in an adjacent,
secondary trophic level (e.g. plants). If species in the secondary level lose
all of their neighbours, they suffer a secondary extinction. Primary
extinctions may occur at random or according to some ordering, such as highest
to lowest degree, where the degree of a species is the number of
interactions/links/edges it has. We can plot the proportion of secondary
species which survive against the proportion of primary extinctions in order
to visualise how robust a given ecological network is against extinction, and
an example of such a “robustness curve” is shown in Figure 1.

Figure 1: (a) An example plant pollinator network and (b) its associated
robustness curve. For this network, pollinators are treated as primary
species, and plants as secondary species. Random primary extinctions are
simulated repeatedly and the proportion of surviving secondary species is
recorded in order to generate the robustness curve.

The area under the robustness curve in Figure 1(b) may be calculated in order
to give a single metric for ecological network robustness, and this is given
by [17]

$R=\frac{1}{N_{p}}\sum_{\varphi=0}^{N_{p}}Pr(\text{survive}|\varphi),$ (1)

where $N_{p}$ is the total number of primary species, $\varphi$ is the number
of primary species which have gone extinct at a given point, and
$Pr(\text{survive}|\varphi)$ is the average probability of a randomly chosen
secondary species surviving after some $\varphi$ primary species have been
removed. Previously, calculations of the robustness curve and the robustness
value $R$ have been done computationally, with some analytic results being
derived for extreme cases [17]. As we show in the following, it is in fact
possible to analytically predict the robustness curve of any given simple
ecological network for a variety of extinction scenarios.

Let us consider some species $A$ in the secondary trophic level, which
initially has degree $k_{A}$ and therefore $k_{A}$ unique neighbours in the
primary level. If species in the primary level go extinct at random, we want
to know the probability that species $A$ has degree $k_{A}-j$ (i.e. $j$
extinct neighbours) after some $\varphi$ number of primary species have gone
extinct. If there are $N_{p}$ primary species, then there are
${N_{p}}\choose{\varphi}$ different possible combinations of primary species
extinctions. We then need to find how many of those combinations include $j$
neighbours of $A$. There are $k_{A}\choose j$ possible combinations for
removing $j$ neighbours of $A$, and therefore there are
$N_{p}-k_{A}\choose\varphi-j$ possible combinations for removing $\varphi-j$
species which are not neighbours of $A$, so long as $\varphi\geq j$.
Multiplying $k_{A}\choose j$ by $N_{p}-k_{A}\choose\varphi-j$ gives us the
total number of combinations of length $\varphi$ which include $j$ neighbours
of $A$, and so we may write the probability of $A$ having degree $k_{A}-j$
once $\varphi$ primary species are extinct as

$Pr(k^{\prime}_{A}=k_{A}-j|\varphi)=\begin{cases}\frac{{k_{A}\choose
j}{N_{p}-k_{A}\choose\varphi-j}}{{N_{p}\choose\varphi}}&\text{if }\varphi\geq
j,\\\ 0&\text{otherwise},\end{cases}$ (2)

where $k^{\prime}_{A}$ refers to $A$’s actual degree value once $\varphi$
primary species have gone extinct. This is the hypergeometric distribution,
which describes a process of sampling without replacement where each sample
may pass (a neighbour of $A$ is removed) or fail (a non-neighbouring primary
species is removed). If we specify that species $A$ goes extinct once its
degree is $k^{\prime}_{A}=k_{A}-i_{k}$ or below (i.e. it has lost at least
$i_{k}$ neighbours), then the disconnection probability for secondary species
$A$ once some $\varphi$ primary species are extinct is

$Pr(A\text{
extinct}|\varphi)=\sum_{j=i_{k}}^{k_{A}}Pr(k^{\prime}_{A}=k_{A}-j|\varphi),$
(3)

Since extinction probability is only dependent upon the initial degree of a
given secondary species species, the total number of primary species in the
network and the number of primary species removed, we may extend this to all
secondary species of initial degree $k$. Consequently, the average secondary
extinction probability over the entire network is

$Pr(\text{extinct}|\varphi)=\sum_{k=0}p(k)\sum_{j=i_{k}}^{k}Pr(k^{\prime}=k-j|\varphi),$
(4)

where $p(k)$ is the probability of some randomly chosen secondary species
having an initial degree of $k$. Given that $Pr(\text{survive}|\varphi)$ is
simply $1-Pr(\text{extinct}|\varphi)$, we can rewrite the expression for
robustness $R$ from Equation 1 as

$R=1-\frac{1}{N_{p}}\sum_{\varphi=0}^{N_{p}}Pr(\text{extinct}|\varphi).$ (5)

Here we note that this analytic model is considerably computationally cheaper
than brute force simulation. With an efficient implementation, calculating
$Pr(\text{extinct}|\varphi)$ analytically takes $O(p)$ time, where $p$ is the
number of unique non-zero entries in the secondary species degree
distribution. By contrast, estimating $Pr(\text{extinct}|\varphi)$
computationally once takes $O(N_{s})$ time, where $N_{s}$ is the number of
secondary species and $N_{s}\geq p$. In practice, it is often necessary to run
several thousand simulations in order to produce an accurate estimate of
$Pr(\text{extinct}|\varphi)$, and so our analytic approach is substantially
computationally cheaper than the brute force method.

In Figure 2(a), we demonstrate the results of our model by comparing the
analytically predicted robustness curve for an ecological network against the
average curve obtained computationally when all neighbours must be removed for
extinction to occur (i.e. $i_{k}=k$). The ecological data used is from a study
of plant pollinator networks in Japan by Kato [19]. We can see that the
computationally obtained curve converges to our predicted curve as the number
of simulations increases, indicating that our method accurately predicts the
average robustness curve.

Figure 2: (a) Analytically predicted and computationally simulated robustness
curves of a real world network from a study by Kato [19] and (b) the absolute
curve divergence between the analytic and simulated curve as computational
simulations increase, plotted on a log-log scale.

We also compare the absolute curve divergence between predicted and simulated
curves for an increasing number of simulations in Figure 2(b). The absolute
curve divergence $D$ is given by

$D=\frac{\sum_{\varphi=0}^{N_{p}}|Pr(\text{survive}|\varphi)_{predict}-Pr(\text{survive}|\varphi)_{sim}|}{N_{p}},$
(6)

where $Pr(\text{survive}|\varphi)_{predict}$ and
$Pr(\text{survive}|\varphi)_{sim}$ are the predicted and average simulated
secondary species survival probabilities respectively. In Figure 2 we can see
that the computationally generated result converges towards our prediction in
the limit of a large number of simulations.

##  3 Minimum Variance Maximises Robustness Against Random Extinctions

Using this analytic model, we can prove that an ecological network is more
robust when the secondary species degree distribution’s variance is minimised,
under the condition that secondary species average degree is held constant. In
other words, the network is most robust if all secondary species have the same
number of links to primary species as one another (or are as close to equal as
possible). Let us consider the scenario in which a secondary species only goes
extinct if it loses all of its primary neighbours. Now let us take two
secondary species $A$ and $B$, where $k_{A}<k_{B}-1$, so if we were to remove
an edge from $B$ and add it to $A$ they would be closer together in degree and
$A$ would not have more edges than $B$. How would rewiring an edge like this
affect their extinction probabilities, and by extension, the network’s
robustness? For a given number of extinct primary neighbours $\varphi$ we can
write the change in average extinction probability for the network as

$\displaystyle\Delta P$
$\displaystyle=Pr(\text{extinct}|\varphi)_{\text{rewired}}-Pr(\text{extinct}|\varphi)_{\text{initial}}$
$\displaystyle=Pr(k^{\prime}_{A}+1=0|\varphi)-Pr(k^{\prime}_{A}=0|\varphi)+Pr(k^{\prime}_{B}-1=0|\varphi)-Pr(k^{\prime}_{B}=0|\varphi)$
$\displaystyle=\frac{1}{{N_{p}\choose\varphi}}\Bigg{[}{N_{p}-(k_{A}+1)\choose\varphi-(k_{A}+1)}-{N_{p}-k_{A}\choose\varphi-
k_{A}}+{N_{p}-(k_{B}-1)\choose\varphi-(k_{B}-1)}-{N_{p}-k_{B}\choose\varphi-
k_{B}}\Bigg{]}.$ (7)

This may be simplified using the identity ${x\choose y}={x-1\choose
y-1}+{x-1\choose y}$ to give

$\Delta
P=\frac{1}{{N_{p}\choose\varphi}}\Bigg{[}{N_{p}-k_{B}\choose\varphi-(k_{B}-1)}-{N_{p}-(k_{A}+1)\choose\varphi-
k_{A}}\Bigg{]},$ (8)

and if we expand the terms for ${N_{p}-k_{B}\choose\varphi-(k_{B}-1)}$ and
${N_{p}-(k_{A}+1)\choose\varphi-k_{A}}$ into their factorial forms, we can
rearrange to give

$\Delta P=\frac{(N_{p}-k_{B})!}{{N\choose\varphi}(\varphi-
k_{A})!(N-\varphi-1)!}\Bigg{[}\prod_{k=k_{A}+1}^{k_{B}-1}(\varphi-(k-1))-\prod_{k=k_{A}+1}^{k_{B}-1}(N_{p}-k)\Bigg{]}.$
(9)

Since
$\prod_{k=k_{A}+1}^{k_{B}-1}(\varphi-(k-1))\leq\prod_{k=k_{A}+1}^{k_{B}-1}(N_{p}-k)$
when $\varphi<N$ and $k_{A}<k_{B}-1$, we know that the change in disconnection
probability $\Delta P\leq 0$ under the same conditions. Therefore, the
disconnection probability must decrease or remain the same under the rewiring
procedure, and so robustness must increase or remain the same. Additionally,
we can prove that this rewiring will always reduce the variance of the
secondary species degree distribution. The change in variance
$\Delta\text{Var}$ is given by

$\displaystyle\Delta\text{Var}$
$\displaystyle=E(k^{2})_{\text{rewired}}-E(k^{2})_{\text{initial}}$
$\displaystyle=\frac{1}{N_{s}}\Big{[}(k_{A}+1)^{2}-k_{A}^{2}+(k_{B}-1)^{2}+k_{B}^{2}\Big{]}$
$\displaystyle=\frac{2}{N_{s}}(k_{A}-k_{B}+1),$ (10)

where $N_{s}$ is the number of secondary species in the network. From this, we
know $\Delta\text{Var}<0$ when $k_{A}<k_{B}-1$, so variance always decreases
for rewiring an edge from species $B$ to species $A$ given that initially
$k_{A}<k_{B}-1$. This proves that, for a secondary species degree distribution
of fixed degree, lower degree distribution variance entails higher robustness,
and vice versa.

This result tells us exactly the structural properties that make secondary
species in ecological networks robust against random primary species
extinctions, namely equally distributed interactions. Previous research has
indicated this before [17], however, this was not conclusively proven, nor was
the relationship between robustness and secondary species degree variance
established. In the work of Burgos et al. [17], robustness is related to
nestedness, which is the propensity of primary species to interact with
secondary species which other, higher degree primary species also interact
with. Nestedness affects the degree distributions of both primary and
secondary species, however, we know from our model that for robustness against
random primary extinctions, the primary species degree distribution is
irrelevant. Therefore, for random primary extinctions, nestedness is not
necessarily an indicator of secondary species robustness.

To illustrate the relation between degree variance and robustness, we provide
a series of example networks in Figure 3, each with equal numbers of
interactions and primary and secondary species, but different robustness and
secondary degree distribution variance. We also show the correlation between
robustness and secondary degree distribution variance for a network undergoing
edge rewiring. The edge rewiring procedure starts on a network where a single
secondary species is connected to all primary species, with all other
secondary species having one primary neighbour, and one by one edges are
swapped from the highest to lowest degree secondary species until secondary
species degree variance is minimised.

Figure 3: (a) Example networks with different second neighbour degree
variances and (b) network robustness plotted against second neighbour degree
variance. Network (i) exhibits the highest variance and lowest robustness,
networks (ii) and (iii) have the same variance and robustness even though they
have different primary species degree distributions, and network (iv) has the
lowest variance and highest robustness.

We can see clearly that the highest variance network (i) has lowest
robustness, and the lowest variance network (iv) has the highest robustness.
Additionally, networks (ii) and (iii) have the same variance and robustness
values as one another, as the only difference between them is the degree
distribution of primary species. The difference in primary species degree
distribution means that they are not considered to have the same nestedness as
one another, but they are equally robust, demonstrating the fact that
nestedness and robustness are not necessarily related.

##  4 Targeted Species Extinctions

The model demonstrated in the preceding sections only predicts robustness when
primary species are removed at random, but it is also possible to adjust the
model to predict robustness when primary species are removed in descending or
ascending degree order. This means that species with many links (descending
order) or few links (ascending order) go extinct first. In these scenarios,
primary species are effectively sorted into some $n$ different groups based on
degree value. All species within a group are removed in a random order before
moving onto the next group which is higher or lower in degree value, depending
on the scenario. If we consider some secondary species, it will have some
$k_{l}$ neighbours in a given primary species group where all primary species
have degree $l$. As before, we set some threshold number $i_{k}$ of
neighbouring species which must be lost before a given secondary species of
degree $k$ goes extinct. We can then say that the secondary species will go
extinct as we remove primary species from some group $d$ if it satisfies the
conditions $\sum_{l=d}^{n}k_{l}\geq i_{k}$ and $\sum_{l=d+1}^{n}k_{l}<i_{k}$
for descending degree order removal, and $\sum_{l=0}^{d}k_{l}\geq i_{k}$ and
$\sum_{l=0}^{d-1}k_{l}<i_{k}$ for ascending degree order removal. We can
therefore write the probability of a some secondary species $A$ going extinct
when primary species are removed in descending degree order as

$Pr(\text{A extinct}|\varphi)=\begin{cases}0&\text{if
}\sum_{l=d}^{n}k_{l}<i_{k}\text{ or }\varphi_{d}<j_{d},\\\ \frac{{k_{d}\choose
j_{d}}{N_{d}-k_{d}\choose\varphi_{d}-j_{d}}}{{N_{d}\choose\varphi_{d}}}&\text{if
}\sum_{l=d}^{n}k_{l}\geq i_{k}\text{ and }\sum_{l=d+1}^{n}k_{l}<i_{k}\text{
and }\varphi_{d}\geq j_{d},\\\ 1&\text{if }\sum_{l=d+1}^{n}k_{l}\geq
i_{k},\end{cases}$ (11)

where $N_{d}$ is the number of primary species in group $d$, $\varphi_{d}$ is
the number of primary species removed from group $d$ and
$j_{d}=i_{k}-\sum_{l=d+1}^{n}k_{l}$. Similarly, for primary species removal in
ascending degree order we have

$Pr(\text{A extinct}|\varphi)=\begin{cases}0&\text{if
}\sum_{l=0}^{d}k_{l}<i_{k}\text{ or }\varphi_{d}<j_{d},\\\ \frac{{k_{d}\choose
j_{d}}{N_{d}-k_{d}\choose\varphi_{d}-j_{d}}}{{N_{d}\choose\varphi_{d}}}&\text{if
}\sum_{l=0}^{d}k_{l}\geq i_{k}\text{ and }\sum_{l=0}^{d-1}k_{l}<i_{k}\text{
and }\varphi_{d}\geq j_{d},\\\ 1&\text{if }\sum_{l=0}^{d-1}k_{l}\geq
i_{k},\end{cases}$ (12)

where $j_{k}=i_{k}-\sum_{l=0}^{d-1}k_{l}$. The extinction probabilities
$Pr(\text{A extinct}|\varphi)$ from Equations 11 and 12 can be averaged over
the degree distribution of the network in a manner similar to Equation 4,
which allows us to predict the average extinction probability for any number
of removed primary species. Therefore, we can predict the robustness curves
for descending and ascending degree removal of primary species, and an example
is given in Figure 4.

Figure 4: Analytically predicted robustness curves for targeted primary
species removal. Predictions are made on the plant pollinator network from a
study by Kato [19]. The blue curve is for removal of primary species in
ascending degree order (lowest degree first) and the orange curve is for
removal of primary species in descending degree order (highest degree first).

As before, our analytic model can successfully predict the average extinction
probabilities for secondary species as primary species are removed. The
scenarios in which primary species are removed in descending and ascending
degree order have been referred to as the “worst” and “best” case scenarios
respectively. However, our analytic model suggests that from the perspective
of robustness, this may not exactly be the case. Under descending degree order
removal, a secondary species’ extinction probability only depends upon how
many lowest degree neighbours it has, and for ascending degree order removal
it depends upon the number of highest degree neighbours. In Figure 5, we
provide an example ecosystem for which descending degree order removal gives
higher network robustness than ascending degree removal.

Figure 5: (a) Example ecological network and (b) targeted removal robustness
curves. Removing low degree primary species (pollinators) first gives a lower
robustness than removing high degree primary species.

While this is a specifically constructed example, it demonstrates that finding
the true worst or best case scenario for secondary extinctions is not
necessarily a case of removing primary species in descending or ascending
degree order respectively. As such, a possible future line of enquiry is to
try and establish the true worst and best case scenarios for species
extinction on any given network.

Thusfar in this section we have considered species extinctions which are
targeted based on degree value, but this is only one possible extinction
ordering. Recent research has examined extinction scenarios in which species
are lost according to their extinction risk as assessed by the IUCN
(International Union for Conservation of Nature) Red List [20]. Species are
ranked in the Red List from Critically Endangered to Least Concern, and in
work by Lamperty and Brosi [20] frugivore species in a seed dispersal network
are removed from highest to lowest extinction risk. Since species only belong
to one of these extinction risk categories, it is necessary to simulate
extinctions from each risk category in descending risk order, with the order
of extinctions within each group randomised. This framework fits well with our
targeted species extinction model.

Using data from Bello et al. [21], we can replicate the results of Lamperty
and Brosi [20], predicting the survival of plant species in a seed dispersal
network as frugivore species are lost in descending extinction risk order. We
calculate extinction probability in the same way as Equation 11, except
instead of our primary species groups being organised by degree value, they
are now organised by extinction risk. In Figure 6, we show simulated and
predicted plant species survival probabilities as frugivore species are
removed in descending extinction risk order, with plant species going extinct
once they have lost all of their frugivore neighbours.

Figure 6: Analytically predicted and computationally simulated robustness
curves of a real world seed dispersal networks, where primary frugivore
species go extinct according to their IUCN extinction risk. This replicates
results from Lamperty and Brosi [20] using data from Bello et al. [21].

This demonstrates the fact that our framework for targeted species extinctions
may be extended to any ordering of primary species loss where primary species
are sorted into groups which go extinct in some order, but extinction within
groups occurs at random.

##  5 Multi-layer Ecosystems

So far, we have only demonstrated our model for bipartite systems, i.e. those
which include only two groups that interact with one another. However, real
world ecosystems can exist on several distinct layers, for example, predators
may feed on pollinators which in turn pollinate plants. Another extension for
our model is to predict species extinction in a group of species not directly
adjacent to the group undergoing primary extinction. This is predictable
analytically, but only for the scenario in which a species must lose all of
its neighbours in order to go extinct.

Scenarios such as this have previously been considered by Pocock et al. [11],
using computational methods. Here we construct an example network to
demonstrate robustness predictions on multi-layer networks. Let us consider a
system of plants, pollinators and predators, where predators feed on
pollinators, who in turn feed on plants. We want to know the probability of a
predator going extinct after a certain number of plant extinctions. For some
predator species A, the species will go extinct if all of the pollinator
species it is connected to go extinct, which only occurs once all of their
plant species neighbours go extinct. Therefore, the extinction probability of
predator species A is simply dependent upon the $u_{A}$ initial number of
unique plant species to which it is connected via its pollinator species
neighbours. Therefore, we can treat the predator species in this system as our
secondary species and the plant species as primary species.

Similar to Equation 2, the extinction probability for a secondary species A
once some $\varphi$ number of primary species have been removed is

$Pr(u^{\prime}_{A}=0|\varphi)=\begin{cases}\frac{{N_{p}-u_{A}\choose\varphi-
u_{A}}}{{N_{p}\choose\varphi}}&\text{if }\varphi\geq u_{A},\\\
0&\text{otherwise},\end{cases}$ (13)

As before, we can average Equation 13 over the distribution of secondary
species connected to $u$ unique primary species to predict secondary
extinctions as primary species are removed. The analytically predicted and
computationally simulated robustness curves for this are given in Figure 7.

Figure 7: (a) Example three layer network of plants (primary species),
pollinators and predators (secondary species), and (b) the associated
analytically predicted and computationally simulated robustness curves.

While this expands the scope of our analytic model beyond simply two layer
ecosystems, it is important to note that thusfar we can only model extinctions
on multi-layer systems if species go extinct after losing all of their
neighbours. Therefore, we cannot consider as many different extinction
scenarios on multi-layer networks than we can for bipartite networks.

##  6 False Positives and Negatives in Network Data

Beyond simply predicting robustness, we may also be interested in how
predictions of robustness are affected by errors in network data, as
ecological data can be error prone [22, 23]. For example, networks may vary
across environmental gradients [24] or may constitute metawebs inferred from
proxies [25, 26]. One may be interested in how the robustness of networks
change as false edges are added in (false positives) or true edges are removed
(false negatives). In the simplest case, let us consider the random addition
and removal of edges. Since robustness against random primary species removal
only depends upon the degree distribution of secondary species, we can
analytically predict how robustness will change as edges are randomly added or
removed by modelling the changes to the secondary degree distribution.

For random edge addition and removal, we can define recursive formulae which
describe how the secondary degree distribution will change. For random edge
addition, the recursive formula which describes the probability of randomly
choosing a secondary species with degree $k$ after some $t$ edges have been
added is

$\displaystyle p(k)_{t}=p(k)_{t-1}$
$\displaystyle-p(k)_{t-1}\frac{N_{p}-k}{N_{p}N_{s}-(E+t-1)}$ (14)
$\displaystyle+p(k-1)_{t-1}\frac{N_{p}-(k-1)}{N_{p}N_{s}-(E+t-1)},$

where $E$ is the total number of edges in the network before any additional
edges have been added. Note that $N_{p}N_{s}-(E+t-1)$ is the total possible
number of edges which could be added to the network once $t-1$ edges have been
added.

For random edge removal, the recursive formula for the probability of choosing
a secondary species with degree $k$ after some $v$ edges have been removed is

$p(k)_{v}=p(k)_{v-1}-p(k)_{v-1}\frac{k}{E-(v-1)}+p(k+1)_{v-1}\frac{k+1}{E-(v-1)}.$
(15)

In order to simulate the presence of false positives or negatives in network
data, we can apply these recursive formulae a certain number of times in order
to adjust the degree distribution, and then assess the impact of false
positives and negatives by predicting robustness for either scenario.

In Figure 8, we show the robustness curve for a plant pollinator network
undergoing random extinctions where $i_{k}=k$, comparing the original
predicted curve against the predicted curves for including false positives and
negatives.

Figure 8: (a) Analytically predicted robustness curves for a real world
ecological network from a study by Kato [19] for the original network, the
network with 200 edges randomly added and the network with 200 edges randomly
removed. The original network has 1125 edges. (b) Shows how network robustness
changes as edges are added or removed, and (c) shows the difference between
the robustness of networks with errors and the original networks as edges are
added or removed.

In terms of robustness, false positives have a more significant impact than
false negatives in small quantities. The network we use to generate the data
shown in Figure 8 has 1125 edges, and we see that up to a change in edges
$\Delta E=500$ (roughly $44\%$ the original number of edges) false positives
increase robustness more than false negatives decrease it. For false
positives, as $R\rightarrow 1$, additional edges contribute less and less to
robustness and for false negatives, as $\Delta E\rightarrow E$, $R\rightarrow
0$, so we see a larger change in robustness due to false negatives for large
numbers of errors. We verify similar results on a dataset of several networks,
and details of this data from 18 real world plant pollinator networks is given
in the Supplementary Materials. For this dataset, we find that if we measure
robustness where $20\%$ of the original number of edges are added/removed for
false positives/negatives, then the net change in robustness is always
positive, i.e. false positives always increase robustness more than false
negatives decrease it.

If we assume that ecological network data gathering in the real world is
reasonably accurate, i.e. unlikely to over/under record interactions by more
than $20\%$, then we would expect false positives to introduce more error into
calculations of robustness than false negatives. Particular care is needed for
robustness analyses based on metawebs of potential trophic interactions for
which the false positive and false negative rates are difficult to ascertain
[25, 26].

While this result indicates that false positives have more impact than false
negatives, it only provides one perspective for how these errors may be
introduced into network data. One future avenue of enquiry is to establish the
likely sources of errors and model those, as opposed to modelling errors
randomly.

##  7 Species Specialisation and Generalisation

In previous sections, we have only considered networks in which interactions
are weighted equally, i.e. each one of a secondary species’ interactions is as
important for its survival. However, on real ecological networks, a secondary
species may interact more with one primary species than another, and this has
an impact on a secondary species survivability [27]. We can specify a certain
percentage of total interaction strength that a species must lose before it
goes extinct, an approach used before before by Schleuning et al. [12]. We can
update our extinction probability for some secondary species A to

$Pr(A\text{
extinct}|\varphi)=\sum_{j=1}^{k_{A}}Pr(k^{\prime}_{A}=k_{A}-j|\varphi)Pr(\sum_{0}^{j}W\geq
i_{A}),$ (16)

where $W$ is a random variable representing some randomly chosen weight
corresponding to the interaction strength with a neighbour of $A$, $i_{A}$ is
the interaction strength threshold for $A$ that must be removed before $A$
goes extinct, and $Pr(\sum_{0}^{j}W\geq i_{A})$ is the probability of choosing
some $j$ weighted interactions which exceed the threshold $i_{A}$. We can get
the value of $i_{A}$ from specifying some ratio of interaction strength that
must be lost for secondary extinction to occur as the sensitivity threshold
$T$, and calculating $i_{A}=\lceil Tk_{A}\rceil$. The notation of $\lceil
x\rceil$ refers to $x$ being rounded up to the nearest integer. The extinction
probability $Pr(A\text{ extinct}|\varphi)$ may be averaged over all species
then over values of $\varphi$ to obtain a robustness value for the network.

There is no closed form solution for $Pr(\sum_{0}^{j}W\geq i_{A})$ since the
weights $w_{z}$ do not necessarily follow a particular distribution. It is
instead necessary to estimate $Pr(\sum_{0}^{j}W\geq i_{A})$ in some way. The
brute force method is to randomly sample $j$ weights using a Monte Carlo
method, however, this must be repeated many times in order to give an accurate
estimate, and is subject to statistical fluctuations.

Instead, we have developed a deterministic sampling method, where a species’
weights $w_{z}$ are arranged in size order and assigned a variable $y_{z}$,
which takes values of 0 or 1. For $j$ removals there will be $j$ values of
$y_{z}=1$, with the rest equal to 0. We can express the sequence of weights as
a sequence of 0 or 1 $y_{z}$ values, giving us a binary number. If weights
were ordered as powers of 2, i.e. $w_{z}=2^{k_{A}-z}$, then we could find the
$n^{th}$ binary sequence of $y_{z}$ values with $j$ values equal to 1 above
which all $\sum_{z=0}^{k_{A}}w_{z}y_{z}\geq i_{A}$ and below which all
$\sum_{z=0}^{k_{A}}w_{z}y_{z}<i_{A}$, allowing us to calculate
$Pr(\sum_{0}^{j}W\geq i_{A})$ exactly. However, weights are not typically
ordered as powers of two, so finding an exact result this way is rarely
possible. Instead we can order binary sequences of $y_{z}$ values and sample
from these orderings at some specified “depth” in order to estimate
$Pr(\sum_{0}^{j}W\geq i_{A})$, where depth effectively determines how many
samples are taken. This sampling method is fully deterministic, so for a given
sequence of weights and a specified depth, we always return the same estimate
for $Pr(\sum_{0}^{j}W\geq i_{A})$. Further details about how this algorithm
operates are given in the Supplementary Materials.

Using our deterministic sampling method, we can provide quasi-analytic
predictions for secondary species survival on networks undergoing random
primary extinctions, where interaction strength is weighted unevenly and
secondary extinctions occur after the loss of a certain percentage of
interaction strength. Example predictions are given in Figure 9, alongside
results showing how our deterministic estimate becomes increasingly accurate
with greater depth, and a comparison between the time taken to estimate
$Pr(\sum_{0}^{j}W\geq i_{A})$ and prediction accuracy for our deterministic
sampling method and for a brute force Monte Carlo method.

Figure 9: (a) Analytically predicted and computationally simulated robustness
curves for a real world network from a study by Kato [19] where unevenly
weighted interaction strength is taken into account and extinctions occur over
a specified threshold $T$ of interaction strength loss. Analytic predictions
are given for threshold values of $70\%,50\%$ and $30\%$, and computationally
simulated curve is given for $50\%$. (b) Comparison between depth of the
estimation for $Pr(\sum_{0}^{j}W\geq i_{A})$ and divergence between
analytically predicted robustness curve and the simulated curve averaged over
5000 iterations. (c) Divergence between the predicted robustness curve and the
5000 iteration simulation curve compared against the time taken, with data for
both the deterministic estimation method and the Monte Carlo method.

From this, we can see that as the depth of the deterministic estimation
increases, we get diminishing returns in terms of prediction accuracy, and
that it is more computationally efficient to use the deterministic estimation
method as opposed to Monte Carlo simulation in order to get the same level of
prediction accuracy.

Having developed an analytic framework for secondary species extinctions when
interaction strength is weighted unevenly, we can examine some extreme
scenarios of interaction strength weighting. One property of interest in
ecological networks is specialisation [28], where specialist species tend to
interact with a small number of species very strongly, and generalist species
tend to interact with many species evenly. Given an ecological network with a
set number of primary and secondary species, and a set distribution of
interactions, we can examine the most specialist interaction weighting and the
most generalist interaction weighting.

In the most specialist case, each secondary species weights one of its
interactions at close to $100\%$ of its interaction strength, and all others
close to $0\%$. Therefore, a given secondary species $A$ only goes extinct
when it loses the neighbour with which it shares almost all interaction
strength. If a neighbour of $A$ goes extinct, the probability of losing the
heavily weighted neighbour is simply $\frac{1}{k_{A}}$, so
$Pr(\sum_{0}^{j}W\geq i_{A})=\frac{j}{k_{A}}$. This give an extinction
probability for some species $A$ of

$\displaystyle Pr(A\text{ extinct}|\varphi)$
$\displaystyle=\sum_{j=1}^{k_{A}}Pr(k^{\prime}_{A}=k_{A}-j|\varphi)\frac{j}{k_{A}},$
$\displaystyle=\frac{\mathbb{E}(k^{\prime}_{A})}{k_{A}}$
$\displaystyle=\frac{\varphi}{N_{p}},$ (17)

which we derive from the fact that
$\mathbb{E}(k^{\prime}_{A})=k_{A}\frac{\varphi}{N_{p}}$ since
$Pr(k^{\prime}_{A}=k_{A}-j|\varphi)$ describes the hypergeometric
distribution. This results in a robustness value of $R=0.5$, regardless of the
secondary degree distribution, number of primary species or threshold.

For the most generalist case, each secondary species weights all of its
interactions evenly, which means $Pr(\sum_{0}^{j}W\geq i_{A})=0$ when
$j<i_{k}$, and $Pr(\sum_{0}^{j}W\geq i_{A})=1$ when $j\geq i_{k}$. Therefore,
$Pr(A\text{
extinct}|\varphi)=\sum_{j=i_{k}}^{k_{A}}Pr(k^{\prime}_{A}=k_{A}-j|\varphi)$,
the same as Equation 3.

Given these results, when is it more advantageous for a network to be highly
specialist or highly generalist in terms of robustness? Let us consider some
secondary species $A$ which is connected to all primary species in its
network, i.e. $k_{A}=N_{p}$. Therefore, the extinction probability for $A$ is
given by

$\displaystyle Pr(A\text{ extinct}|\varphi)$
$\displaystyle=\sum_{j=i_{k}}^{k_{A}}\begin{cases}\frac{{N_{p}\choose
j}{0\choose\varphi-j}}{{N_{p}\choose\varphi}}&\text{if }\varphi\geq j,\\\
0&\text{otherwise},\end{cases}$ $\displaystyle=\begin{cases}1&\text{if
}\varphi\geq i_{k},\\\ 0&\text{otherwise},\end{cases}$ (18)

If all secondary species in a network have $k=N_{p}$, then when they are
maximally generalist, the network robustness is $R=\frac{i_{k}}{N_{p}}$.
Therefore, such a network is more robust when secondary species are maximally
generalist if more than $50\%$ of interaction strength must be lost before a
secondary species goes extinct, i.e. when the sensitivity threshold $T>0.5$.
Conversely, the network is more robust when secondary species are maximally
specialist if less than $50\%$ of interaction strength must be lost to make
secondary species go extinct, so $T<0.5$. To illustrate this, we provide
robustness curves in Figure 10 of (a) a single secondary species with
$k_{A}=N_{p}$, and of (b) an entire real world network.

Figure 10: (a) Shows species survival for maximal generalists and maximal
specialists at various sensitivity thresholds $T$ when $k_{A}=N_{p}$. (b)
Gives robustness curves for maximal generalisation and maximal specialisation
on a real world network from a study by Kato [19] at various sensitivity
thresholds $T$.

From these results, we know that either extreme of species specialisation can
be advantageous from the perspective of maximising network robustness,
depending upon the sensitivity of the network, i.e. the proportion of
interaction strength that must be lost for secondary species to go extinct.
However, we see in Figure 10 (b) that it is not strictly the case on real
networks that maximum generalisation is always better for robustness than
maximum specialisation when $T>0.5$, as the maximum generalist curve when
$T=0.5$ gives $R=0.471$. This is due to the fact that secondary species
typically have $k<N_{p}$ on real networks. Additionally, we note that the
robustness values from the maximum generalist and maximum specialist
interaction weightings do not necessarily give the maximum and minimum
robustness values for a given threshold. Nevertheless, these results are still
indicative of the fact that species generalisation and specialisation can both
improve network robustness in different contexts, and so we might expect that
in the real world, a network that has developed to be highly generalist is
less sensitive to interaction loss than a network which has developed in order
to be highly specialist.

##  8 Interaction Loss

The models we have considered are only concerned with the loss of primary
species as a whole, where a primary species is removed at each “step” in the
extinction process. Considering the loss of entire species at a time can skew
our understanding of network robustness, for example a plant animal network
with more animal than plants will appear more robust against primary
extinctions of animal than against primary extinctions of plants [12].
However, is this a realistic understanding of how species go extinct? There
may be more animal species than plant species, but what if there is a very
large population of each plant species and a small population of each animal
species? Extinctions may be experienced more gradually, where a species’
population dies off over time rather than all at once [29]. This process can
be modelled by examining the loss of interactions as opposed to the loss of
species, which in network terms entails considering edge removal as opposed to
node removal.

If interaction strength between species is represented as integer values, then
we can treat each unit of interaction strength as an edge, so secondary
species have degree values equal to the sum of their interaction strength with
other species. We then have $E$ “edges” (i.e. total interaction strength on
the network), and we remove some $\varphi$ units of interaction strength. For
a given secondary species $A$, after removing some $\varphi$ interaction
strength the probability that it has lost some interaction strength $j$ is

$Pr(k^{\prime}_{A}=k_{A}-j|\varphi)=\begin{cases}\frac{{k_{A}\choose
j}{E-k_{A}\choose\varphi-j}}{{E\choose\varphi}}&\text{if }\varphi\geq j,\\\
0&\text{otherwise},\end{cases}$ (19)

where, as before in the case of Equation 2, $k_{A}$ is the initial
degree/total interaction strength of $A$, and $k_{A}^{\prime}$ is the
degree/total interaction strength of $A$ after removing $\varphi$ interaction
strength. From this, it is straightforward to predict secondary species’
survival probability and network robustness using a similar logic as in
Section 2. Predictions for interaction loss on an ecological network are given
in Figure 11.

Figure 11: Robustness curves for interaction loss on a real world network from
a study by Kato [19] at varying sensitivity thresholds $T$.

These predictions of interaction loss give different robustness values than
predictions of species extinctions. For example, for species extinctions on a
real network [19] when true interaction strength values are used (as shown in
Figure 9 (a)) and $T=0.7$, we have $R_{species}=0.619$. By contrast, for
interaction loss on the same network when $T=0.7$, we have $R_{inter}=0.653$.
Therefore, modelling secondary species extinctions as an outcome of
interaction loss as opposed to primary species extinctions gives a different
perspective on network robustness, allowing one to identify networks which are
fragile against primary species loss but robust against interaction loss, and
vice versa. A similar logic to that presented in Section 2 may be followed in
order to show that a network with a set number of primary species, secondary
species and total interaction strength is maximally robust against interaction
loss when the variance in interaction strength per secondary species is
minimised.

##  9 Discussion

In conclusion, we have successfully extended the robustness framework of
Memmott et al. [10] such that we may make predictions of ecological network
robustness analytically. For random extinctions, we have shown that networks
with low second species degree variance are highly robust. We are also able to
predict secondary extinctions as primary species go extinct according to some
degree or extinction risk based targeting, and we can predict secondary
extinctions on ecological networks with more than two layers. Additionally, we
can model the influence of random false positives and negatives in network
data on robustness, finding that in small quantities false positives have a
greater impact than false negatives on network robustness.

Our model is also capable of predicting the robustness of networks where
interaction strength is weighted unevenly between different secondary species’
neighbours, and species go extinct once a certain proportion of interaction
strength has been lost. We have given results for robustness when interaction
strength is equally distributed (maximally generalist), and when interaction
strength is shared solely with one neighbouring species (maximally
specialist). From this, we know that maximal generalisation and maximal
specialisation can both produce a more robust network, depending on the
proportion of interaction strength that must be lost before secondary species
extinction. Finally, we have demonstrated the fact that it is also possible to
model interaction strength loss as opposed to simply species extinction,
representing a more “gradual” extinction scenario.

These results represent a substantial advancement in analytic understanding of
ecological network robustness. However, there are still many open questions.
We can predict the average secondary species extinction probability for a
given number of primary extinctions, but we may also want to analytically
predict the possible error in robustness curves by finding the standard
deviation in secondary species extinction probability for a given number of
primary extinctions. Additionally, we may want to establish the true worst and
best case scenarios for secondary extinctions, as these have not been
definitively identified. For errors (i.e. false positives and negatives) in
network data, our current results examine errors which occur at random, but
this may not be the case in the real world. Errors may occur due to some
specific reason or dynamic, and identifying what this is may allow us to
better mathematically model data errors and their influence.

Beyond these possible improvements, it is also important to acknowledge that
in recent years, ecologists have considered properties of ecological networks
which affect robustness and go beyond simpler models of species extinction.
For example, ontogenetic niche shifts, where species change their diets when
they undergo changes such as growing from a larvae to an adult, can affect the
structure and robustness of interaction networks [30]. Another consideration
is how interactions can be “rewired” after species extinctions [12, 31], which
to predict analytically would likely require combinatoric methods for sampling
with fuzzy replacement [32]. These examples highlight the fact that there is
still considerable room for analytic models of ecological network robustness
to develop, and there are ongoing areas of research in both ecological
networks and combinatorics which may complement one another well, so it may be
useful for there to be a greater dialogue between these fields in future.

## Funding

C.J. was supported by Engineering and Physical Sciences Research Council
Doctoral Training Partnership funding (EP/R513179/1). We thank the Research
Focus Data-centric Sciences of the University of Potsdam for financial
support.

## Acknowledgements

We thank Jane Memmott for providing feedback on a summary version of this
paper.

## Conflict of Interest Statement

There authors have no conflict of interest to declare.

## Data Availability

The code used in generating the results for this paper may be found online at
https://github.com/cj14373/eco-analytic.git

## References

  * [1] E. Haeckel, Generelle morphologie der organismen. Allgemeine grundzüge der organischen formen-wissenschaft, mechanisch begründet durch die von Charles Darwin reformirte descendenztheorie.  Berlin, G. Reimer, 1866, 1866.
  * [2] J. E. Cohen, Food Webs and Niche Space. (MPB-11), Volume 11.  Princeton University Press, 1978.
  * [3] L.-F. Bersier, A history of the study of ecological networks, pp. 365–421.  World Scientific, 12 2007.
  * [4] T. C. Ings and J. E. Hawes, The History of Ecological Networks, pp. 15–28.  Cham: Springer International Publishing, 2018.
  * [5] T. P. Dawson, S. T. Jackson, J. I. House, I. C. Prentice, and G. M. Mace, “Beyond predictions: Biodiversity conservation in a changing climate,” Science, vol. 332, no. 6025, pp. 53–58, 2011.
  * [6] C. Bellard, C. Bertelsmeier, P. Leadley, W. Thuiller, and F. Courchamp, “Impacts of climate change on the future of biodiversity,” Ecology Letters, vol. 15, no. 4, pp. 365–377, 2012.
  * [7] R. H. Jongman, “Nature conservation planning in europe: developing ecological networks,” Landscape and urban planning, vol. 32, no. 3, pp. 169–183, 1995\.
  * [8] M. L. Forup, K. S. Henson, P. G. Craze, and J. Memmott, “The restoration of ecological interactions: plant–pollinator networks on ancient and restored heathlands,” Journal of Applied Ecology, vol. 45, no. 3, pp. 742–752, 2008\.
  * [9] J. M. Tylianakis, E. Laliberté, A. Nielsen, and J. Bascompte, “Conservation of species interaction networks,” Biological conservation, vol. 143, no. 10, pp. 2270–2279, 2010.
  * [10] J. Memmott, N. M. Waser, and M. V. Price, “Tolerance of pollination networks to species extinctions,” Proceedings: Biological Sciences, vol. 271, no. 1557, pp. 2605–2611, 2004.
  * [11] M. J. O. Pocock, D. M. Evans, and J. Memmott, “The robustness and restoration of a network of ecological networks,” Science, vol. 335, no. 6071, pp. 973–977, 2012.
  * [12] M. Schleuning, J. Fründ, O. Schweiger, E. Welk, J. Albrecht, M. Albrecht, M. Beil, G. Benadi, N. Blüthgen, H. Bruelheide, K. Böhning-Gaese, M. Dehling, C. Dormann, N. Exeler, N. Farwig, A. Harpke, T. Hickler, A. Kratochwil, M. Kuhlmann, and C. Hof, “Ecological networks are more sensitive to plant than to animal extinction under climate change,” Nature Communications, vol. 7, 12 2016.
  * [13] R. V. Solé and J. M. Montoya, “Complexity and fragility in ecological networks,” Proceedings: Biological Sciences, vol. 268, no. 1480, pp. 2039–2045, 2001.
  * [14] J. Montoya, S. Pimm, and R. Sole, “Ecological networks and their fragility,” Nature, vol. 442, pp. 259–64, 08 2006.
  * [15] P. A. Aguilera, A. Fernández, R. Fernández, R. Rumí, and A. Salmerón, “Bayesian networks in environmental modelling,” Environmental Modelling & Software, vol. 26, no. 12, pp. 1376–1388, 2011.
  * [16] P. Ramazi, M. Kunegel-Lion, R. Greiner, and M. A. Lewis, “Exploiting the full potential of bayesian networks in predictive ecology,” Methods in Ecology and Evolution, vol. 12, no. 1, pp. 135–149, 2021.
  * [17] E. Burgos, H. Ceva, R. P. Perazzo, M. Devoto, D. Medan, M. Zimmermann, and A. María Delbue, “Why nestedness in mutualistic networks?,” Journal of Theoretical Biology, vol. 249, no. 2, pp. 307–313, 2007.
  * [18] A. Nielsen and J. Bascompte, “Ecological networks, nestedness and sampling effort,” Journal of Ecology, vol. 95, no. 5, pp. 1134–1141, 2007.
  * [19] M. Kato, “Anthophilous insect community and plant-pollinator interactions on amami islands in the ryukyu archipelago, japan,” Contributions from the Biological Laboratory, Kyoto University, vol. 29, no. 2, pp. 157–254, 2000.
  * [20] T. Lamperty and B. J. Brosi, “Loss of endangered frugivores from seed dispersal networks generates severe mutualism disruption,” Proceedings of the Royal Society B, vol. 289, no. 1984, p. 20220887, 2022.
  * [21] C. Bello, M. Galetti, D. Montan, M. Pizo, T. Mariguela, L. Culot, F. Bufalo, F. Labecca, F. Pedrosa, R. Constantini, C. Emer, W. Silva, F. Da Silva, O. Ovaskainen, and P. Jordano, “Atlantic frugivory: A plant-frugivore interaction dataset for the atlantic forest,” Ecology, vol. 98, 03 2017\.
  * [22] A. Kangas, T. Packalen, K. Korhonen, and J. Vauhkonen, “Sources and types of uncertainties in the information on forest-related ecosystem services,” Forest Ecology and Management, vol. 427, 05 2018.
  * [23] M. A. de Aguiar, E. A. Newman, M. M. Pires, J. D. Yeakel, C. Boettiger, L. A. Burkle, D. Gravel, P. R. Guimarães Jr, J. L. O’Donnell, T. Poisot, et al., “Revealing biases in the sampling of ecological interaction networks,” PeerJ, vol. 7, p. e7566, 2019.
  * [24] L. Pellissier, C. Albouy, J. Bascompte, N. Farwig, C. Graham, M. Loreau, M. A. Maglianesi, C. J. Melián, C. Pitteloud, T. Roslin, et al., “Comparing species interaction networks along environmental gradients,” Biological Reviews, vol. 93, no. 2, pp. 785–800, 2018.
  * [25] I. Morales-Castilla, M. G. Matias, D. Gravel, and M. B. Araújo, “Inferring biotic interactions from proxies,” Trends in ecology & evolution, vol. 30, no. 6, pp. 347–356, 2015.
  * [26] L. Maiorano, A. Montemaggiori, G. F. Ficetola, L. O’connor, and W. Thuiller, “Tetra-eu 1.0: a species-level trophic metaweb of european tetrapods,” Global Ecology and Biogeography, vol. 29, no. 9, pp. 1452–1457, 2020.
  * [27] E. L. Berlow, S. A. Navarrete, C. J. Briggs, M. E. Power, and B. A. Menge, “Quantifying variation in the strengths of species interactions,” Ecology, vol. 80, no. 7, pp. 2206–2224, 1999.
  * [28] N. Blüthgen, F. Menzel, and N. Blüthgen, “Measuring specialization in species interaction networks,” BMC ecology, vol. 6, p. 9, 02 2006.
  * [29] A. Valiente-Banuet, M. A. Aizen, J. M. Alcántara, J. Arroyo, A. Cocucci, M. Galetti, M. B. García, D. García, J. M. Gómez, P. Jordano, R. Medel, L. Navarro, J. R. Obeso, R. Oviedo, N. Ramírez, P. J. Rey, A. Traveset, M. Verdú, and R. Zamora, “Beyond species loss: the extinction of ecological interactions in a changing world,” Functional Ecology, vol. 29, no. 3, pp. 299–307, 2015.
  * [30] T. Nakazawa, “Ontogenetic niche shifts matter in community ecology: a review and future perspectives,” Population Ecology, vol. 57, no. 2, pp. 347–354, 2015.
  * [31] K. C. Baldock, M. A. Goddard, D. M. Hicks, W. E. Kunin, N. Mitschunas, H. Morse, L. M. Osgathorpe, S. G. Potts, K. M. Robertson, A. V. Scott, et al., “A systems approach reveals urban pollinator hotspots and conservation opportunities,” Nature ecology & evolution, vol. 3, no. 3, pp. 363–373, 2019.
  * [32] O. Kesemen, B. Tiryaki, Ö. Tezel, E. Özkul, and E. Naz, “Random sampling with fuzzy replacement,” Expert Systems with Applications, vol. 185, p. 115602, 07 2021.