Saint- Clair Chabert-Liddell
email: saint-clair.chabert-liddell@inrae.fr

Pierre Barbillon

Sophie Donnet

Impact of the mesoscale structure of a bipartite ecological interaction network on its robustness through a probabilistic modeling *

Keywords: Ecological network, Bipartite Stochastic Block Model, Partially Observed Network

The robustness of an ecological network quantifies the resilience of the ecosystem it represents to species loss. It corresponds to the proportion of species that are disconnected from the rest of the network when extinctions occur sequentially. Classically, the robustness is calculated for a given network, from the simulation of a large number of extinction sequences. The link between network structure and robustness remains an open question.

Setting a joint probabilistic model on the network and the extinction sequences allows analysis of this relation. Bipartite stochastic block models have proven their ability to model bipartite networks e.g. plant-pollinator networks: species are divided into blocks and interaction probabilities are determined by the blocks of membership. Analytical expressions of the expectation and variance of robustness are obtained under this model, for different distributions of primary extinction sequences. The impact of the network structure on the robustness is analyzed through a set of properties and numerical illustrations. The analysis of a collection of bipartite ecological networks allows us to compare the empirical approach to our probabilistic approach, and illustrates the relevance of the latter when it comes to computing the robustness of a partially observed or incompletely sampled network.

Introduction

In response to the rapid evolution of ecosystems due to climate change (habitat losses and species extinctions) on the one hand, and the increasing number of available data sets of ecological interaction networks on the other hand, the study of the robustness of ecological networks to species loss has become an active area of research in the ecological scientific community (see Landi et al., 2018, for a review). Given a primary extinction sequence, the influence of these extinctions on the other species is studied by monitoring the proportion of species that get disconnected from the rest of the network [START_REF] Dunne | Network structure and biodiversity loss in food webs: robustness increases with connectance[END_REF].

Bipartite networks are used to represent interactions between two separated sets of species that do not interact within their own set. These interactions may be either mutualistic when the species of both groups benefit from the interactions (such as pollination for plant-pollinator networks or seed-dispersal for plant-frugivore bird networks) or antagonistic when one group of species benefits from the interactions at the expense of the other group of species (e.g. plantherbivore or host-parasite networks). Robustness is a numeric indicator quantifying the impact of the disappearance of one set of (primary) species on the other (secondary species) under given extinction sequences (see [START_REF] Memmott | Tolerance of pollination networks to species extinctions[END_REF][START_REF] Curtsdotter | Robustness to secondary extinctions: comparing traitbased sequential deletions in static and dynamic food webs[END_REF], for examples of primary extinction sequences). In a few words, the object of interest is the proportion of secondary species that remain connected to other species after a given number of primary species have disappeared. Since this proportion depends on the order in which the species disappear, robustness is defined as the average of these proportions over a large number of primary extinction sequences. 1The definition of ecological robustness used in this paper is not the only one, in particular the robustness of ecological networks is also studied through dynamic approaches. Among recent developments in the dynamic approach, [START_REF] Song | A guideline to study the feasibility domain of multi-trophic and changing ecological communities[END_REF] are interested in the feasibility domain of ecological communities, i.e. the set of environmental conditions under which all species have positive abundances. [START_REF] Barabás | Sensitivity analysis of coexistence in ecological communities: theory and application[END_REF] estimate extinction risk to small environmental variation through the use of sensitivity analysis while others study the capacity of coexistence of species through the so-called invasion criterion (see [START_REF] Grainger | The invasion criterion: a common currency for ecological research[END_REF], for a review).

Obviously, the robustness defined above is strongly related to the size of the network consid-ered (i.e. the number of species involved or species richness) and to its density (i.e. the number of interactions observed compared to the total number of possible interactions or connectance).

To go further, the question arises to what extent the topological properties of a network also influence its robustness. Indeed, observed ecological networks present different topological structures depending on the type of interactions: for example, mutualistic networks are known to have a strong nested structure, while this is not necessarily the case for antagonistic networks [START_REF] Fortuna | Nestedness versus modularity in ecological networks: two sides of the same coin[END_REF][START_REF] Bascompte | The nested assembly of plant-animal mutualistic networks[END_REF]. A key question is to identify the relationship between the ability of a network to withstand species extinctions and its characteristics such as its size and its mesoscale structure.

To that purpose, we propose to assume a parametric probabilistic model for the bipartite network. This model will embed the topological properties of the network in a few number of parameters. Then we propose to study the behavior of the robustness in this framework and in particular the variations of the robustness with respect to the model parameters. More precisely, we focus on the expectation and the variance of the robustness under a network probabilistic model. For some particular probabilistic models and some particular extinction sequences, the relation between the parameters and the expected robustness can be provided in a closed-form, making a fine study of the relation possible.

Furthermore, relying on a probabilistic model which can be adjusted to account for an observation process may enable for correction of the sampling effect on the robustness. Indeed, even though ecological networks are often considered to describe all the possible interactions between species, the sampling may be incomplete [START_REF] Blüthgen | What do interaction network metrics tell us about specialization and biological traits[END_REF] which may bias the computed network statistics [START_REF] Rivera-Hutinel | Effects of sampling completeness on the structure of plant-pollinator networks[END_REF][START_REF] De Manincor | How biased is our perception of plant-pollinator networks? a comparison of visit-and pollen-based representations of the same networks[END_REF] such as the robustness.

The stochastic block model (SBM [START_REF] Nowicki | Estimation and prediction for stochastic blockstructures[END_REF] and its extensions for bipartite networks, such as the biSBM [START_REF] Govaert | Clustering with block mixture models[END_REF], also referred to as Latent Block Model) or its degree corrected counterpart (DCbiSBM [START_REF] Larremore | Efficiently inferring community structure in bipartite networks[END_REF] have gained a lot of attention in the statistical and network science research fields during the last decade. While the Erdős-Rényi model [START_REF] Erdős | On the evolution of random graphs[END_REF] assumes that any pair of species has the same probability of interaction, SBMs introduce heterogeneity in the connection behavior. Specifically, in bipartite SBMs, the two sets of nodes are divided into clusters/blocks/groups and the probability that two nodes are connected depends on which blocks the nodes belong to. Depending on the number of blocks, their size, and interaction probabilities (inter-and intra-block), SBMs encompass a wide variety of topologies (such as assortative community or core-periphery) and encode them in a small number of parameters.

Since [START_REF] Allesina | Food web models: a plea for groups[END_REF] has advocated for the use of groups in ecological networks, SBMs (sometimes referred to as group models) have gained in popularity. Some vari-ants adapted to multilayer ecological networks have been proposedfor: multiplex networks [START_REF] Kéfi | How structured is the entangled bank? the surprisingly simple organization of multiplex ecological networks leads to increased persistence and resilience[END_REF], multipartite networks [START_REF] Bar-Hen | Block models for generalized multipartite networks: Applications in ecology and ethnobiology[END_REF] or temporal networks [START_REF] Matias | Statistical clustering of temporal networks through a dynamic stochastic block model[END_REF]. Besides, they have been used to answer specific ecological questions. To name a few, [START_REF] Michalska-Smith | Understanding the role of parasites in food webs using the group model[END_REF] explore the structural role of parasite species in food webs, while [START_REF] Miele | Core-periphery dynamics in a plantpollinator network[END_REF] use biSBM to assess the core-periphery structure of plant-pollinator networks. Furthermore, SBMs provide an ecological interpretation in terms of functional groups in the ecosystem: species in the same block interact in a similar way, which means that the exchangeability of species in a block model is related to the concept of ecological equivalence [START_REF] Sander | What can interaction webs tell us about species roles?[END_REF].

To our knowledge, the behavior of the robustness has never been studied in the SBM framework. Some grouping algorithms are used to derive extinction sequences. [START_REF] Cai | The robustness of ecosystems to the species loss of community[END_REF] optimize an objective function to determine communities in the network and then generates primary extinction based on these communities. More generally, some models are used to generate primary extinction sequences, or to model how species of the other functional group react to those extinctions (such as rewiring or cascading, see [START_REF] Bane | Effects of model choice, network structure, and interaction strengths on knockout extinction models of ecological robustness[END_REF]Vizentin-Bugoni et al., 2020, for examples). In an approach more similar to ours, [START_REF] Burgos | Why nestedness in mutualistic networks?[END_REF] studied the relationship between nestedness and robustness of mutualistic networks by using the self-organizing network model [START_REF] Medan | Analysis and assembling of network structure in mutualistic systems[END_REF] to generate nested networks and derived analytical expression of the robustness under this model.

The robustness given as a decreasing function and a robustness statistic classically used in ecology are described in Section 2. Section 3 is dedicated to the introduction of the biSBM, the DCbiSBM and related models for sequential species extinctions. Section 4 supplies the expression of the expectation and the variance of the robustness under a biSBM for different distributions of primary extinction sequences. The analytical properties of the expected robustness are given in Section 5 together with more general studies to illustrate the impact of the network structure on the robustness. In Section 6, we apply our approach on a dataset composed of both mutualistic and antagonistic bipartite networks and compare our results to the classical approach. Finally, on the same dataset, we show how our approach allows us to calculate the expected robustness of partially observed ecological networks.

Robustness of bipartite ecological networks

Robustness aims at measuring the tolerance of a network to species extinctions by quantifying the proportion of remaining species along a species extinction sequence. Formally, let A ∈ {0, 1} nr×nc be the n r × n c incidence matrix of a bipartite network representing ecological interactions between two groups of species of respective sizes n r and n c (such as plant-pollinators or hosts-parasites). Then:

A ij =    1 if row species i interacts with column species j, 0 otherwise.
Let s be an extinction sequence on the row species: s ∈ S nr where S nr is the symmetric group. The extinctions in row (whose order is given by s) are the primary extinctions. These row extinctions lead to secondary extinctions among the column species if these column species remain isolated after the disappearance of row species. More precisely, a column species j is said to be be extinct after m row primary extinctions if these m primary extinctions caused species j to lose all its connections, or equivalently if after these m primary row extinctions, j has no connections left with the remaining rows, which is equivalent to nr i=m+1 A s(i)j = 0. For a given sequence s and a given number of primary extinctions m, R(A, s, m) is the proportion of remaining column species:

R(A, s, m) = 1 - 1 n c nc j=1 1 { nr i=m+1 A s(i)j =0} . (1)
Note that for all m ∈ {0, . . . , n r }, 0 ≤ R(A, s, m) ≤ 1 with R(A, s, n r ) = 0. The robustness function is defined as the expectation of R(A, S, m) against the primary extinction sequences:

m → R S (A, m) = E S [R(A, S, m)] with S ∼ S . ( 2
)
where S is a probability distribution on S nr . Equation ( 2) is a weighted summation over the n r ! possible permutations of S nr , which may render the computation intractable. In practice, it is generally approximated by a Monte Carlo integration:

R S (A, m) = 1 B B b=1 R(A, S (b) , m) where S (b) ∼ i.i.d S, b = 1, . . . , B. (3)
The robustness statistic of a network A and a primary extinction sequences distribution S is defined as:

R S (A) = 1 n r nr m=0 R S (A, m), (4)
which corresponds to the area under the curve (AUC) of the robustness function where the x-axis has been re-normalized to match with the proportion of removed row species. R S (A) is the Monte Carlo version of R S (A).

Other statistics exist in the literature, such as the median of the robustness function [START_REF] Dunne | Network structure and biodiversity loss in food webs: robustness increases with connectance[END_REF], i.e. the proportion of primary extinctions needed to provoke a secondary extinction on half of the species. However, the statistic (4) is widely used in ecology and is mathematically convenient for our purpose.

About the extinction sequence distribution S Several choices for S corresponding to various ecological scenarios are commonly considered in the literature. The first choice is the uniform distribution over S nr (S = U Sn r ), assuming that the species disappear without specific order. For the sake of simplicity, we use U instead of U Sn r .

Distributions S depending on A are suggested in the literature (see [START_REF] Curtsdotter | Robustness to secondary extinctions: comparing traitbased sequential deletions in static and dynamic food webs[END_REF], for examples). Among these approaches, we will focus on sequences depending on the row degree sequences. For any row species i = 1, . . . , n r , let D i = nc j=1 A ij be its degree, i.e. the number of edges involving i. On the one hand, the worst-case scenario in terms of ecological extinctions assumes that the row species with the highest degrees (most connected row species) disappear first. In this case, the species in rows are ordered in decreasing degrees and the primary extinction sequences follow this order; row species of equal degrees disappear in a uniformly distributed order. On the other hand, the generation of primary extinction sequences that first eliminate species of lower degree mimics a more favorable ecological scenario. In these two cases, the support of S is of cardinal k #{i : D i = k}!. Depending on the sequence (D i ) i=1,...,nr , the corresponding R S (A, m) may become tractable or not. If not, a Monte Carlo approximation is used.

Instead of considering a strict monotone ordering of the row degrees, one might relax the constraint and set the probability for a row species to disappear proportional to a function of its degree, as seen in [START_REF] Liu | Node attacks on two-mode networks[END_REF]. This would correspond to sampling without replacement a sequence s where the weights of each species i ∈ {1, . . . , n r } are, for instance given by:

w i ∝ D α i . (5)
α = 0 coincides to the uniform distribution U while, if α = 1, the primary extinction sequence depends linearly on the degrees. The increasing order is obtained by reversing this sequence of primary extinction.

Remark 1. The R S (A, m)'s computed for the first two S's described above are the most widely used. They are implemented in the R package bipartite [START_REF] Dormann | Introducing the bipartite package: analysing ecological networks[END_REF] available on cran.

The robustness function (2) and corresponding statistic (4) are then computed from a unique observed network A and a probability distribution S conditional to A (through the row degrees):

R S (A, m) = E S|A [R(A, S, m)] with S|A ∼ S . (6)
This robustness computed for a given A will be referred to as the empirical robustness.

In order to understand the variability of the robustness with respect to the network structure we set a probabilistic model on the network A ∼ A θ where the parameters θ embeds the network structure in a small number of parameters and then study the random variable R S (A, m) for various distributions on (S, A). The following section is dedicated to the description of flexible probabilistic models on A and adapted joint distributions on A and S.

3 Bipartite Stochastic Block Model and related sequential extinctions

3.1 Probabilistic model on bipartite ecological networks Bipartite stochastic block models The bipartite SBM, a.k.a. the Latent Block Model [START_REF] Govaert | Clustering with block mixture models[END_REF], is a mixture model on the edges adapted to bipartite networks. It relies on a clustering of nodes in rows and a clustering of nodes in columns. The nodes which belong to the same cluster (or equivalently block) are assumed to share the same connectivity profile in the network and the probability of interaction between two nodes depends on the blocks they belong to. More precisely, each of the n r row species is attributed to a block k ∈ {1, . . . , Q r } independently from the other species. Let Z i be such that Z i = k if row species i belongs to block k. The Z i 's are assumed to be independent and identically distributed (i.i.d.) and:

P(Z i = k) = π k i ∈ {1, . . . , n r }, (7)
with Qr k=1 π k = 1. For j = 1, . . . , n c , let W j be such that W j = q if column species j belongs to block q ∈ {1, . . . , Q c }. The W j 's are assumed to be i.i.d. and:

P(W j = q) = ρ q j ∈ {1, . . . , n c }, (8)
with Qc q=1 ρ q = 1. Then, conditionally to their respective latent blocks, the interactions between two species are distributed independently as:

A ij |{Z i = k, W j = q} ∼ B(δ kq ), ( 9
)
where B is the Bernoulli distribution and δ kq ∈ (0, 1) for (k, q) ∈ {1, . . . , Q r } × {1, . . . , Q c }.

The parameter θ = {ρ q , π k , δ kq } k=1,...,Qr,q=1,...,Qc then encodes the topology of the network. Let n = (n r , n c ) be the numbers of species (richness) in rows and in columns of the network:

biSBM (θ, n) denotes the distribution on the networks defined by equations ( 7), ( 8) and ( 9) for a given value of θ. Note that d = k,q π k δ kq ρ q is the expected connection probability for any pair of species (i, j). d will be referred to as the expected density. . Remark 2. A special case among the biSBM distributions is the one where δ kq = d for any k, q, or equivalently -in terms of models on networks-where Q r = 1 and Q c = 1. In that case, all connections occur independently with the same probability d. This biSBM is the bipartite Erdős-Rényi network.

Degree corrected stochastic block models The degree corrected stochastic block model [START_REF] Karrer | Stochastic blockmodels and community structure in networks[END_REF]) is an extension of the SBM which accounts for the heterogeneity in the degree distribution. This heterogeneity does not only depend on the blocks but also on some parameters related to the nodes themselves. The extension to bipartite network is done in [START_REF] Larremore | Efficiently inferring community structure in bipartite networks[END_REF]. In these models, the distribution on the dyads is a Poisson distribution even if the observed network consists of binary edges. This is a reasonable approximation in the case where the network is large and sparse. However, in ecological interaction networks, this assumption is hardly met. That is why we use a true binary definition of the degree corrected biSBM that we denote by DCbiSBM in the following. On top of the blocks for the nodes in rows and columns given by variables Zs and W s as defined in Equations ( 7) and ( 8), two vector parameters γ r ∈ R nr and γ c ∈ R nc are associated to the nodes and control their degree. The larger γ r,i (resp. γ c,j ) the higher the probability of connection involving species i (resp. j), compared to other species within the same block. Then, Equation ( 9) is replaced with

A ij |{Z i = k, W j = q} ∼ B(1/(1 + exp(-(δ kq + γ r,i + γ c,j ))) , (10)
with γ r,1 = γ c,1 = 1 for identifiability issues. We denote by DCbiSBM (θ, γ, n) where γ = (γ r , γ c ), this probabilistic model with a given set of parameters. If there is only one block, i.e. Q r = Q c = 1, this model only accounts for the heterogeneity of degrees and so corresponds to a bipartite version of the expected degree distribution (EDD) model [START_REF] Chung | Connected components in random graphs with given expected degree sequences[END_REF].

Extinction sequence distributions adapted to bipartite Block Models

As described in Section 2, classically, the extinction sequences are either distributed uniformly on S nr or conditionally to A. Positing a probabilistic distribution on A leads to a joint distribution on A and S. In the case of a uniform distribution on S: S ∼ U, the joint distribution on (A, S) consists of the product of the distributions for A and S. We will then denote by L θ,n,U and L θ,γ,n,U the corresponding joint distributions when

A ∼ biSBM (θ, n) or A ∼ DCbiSBM (θ, γ, n) respectively.
As described in Section 2, the extinction sequences may depend on A through its degree distributions. When A follows a probabilistic distribution, such an extinction sequence is defined conditionally on the realization of A. We propose another dependence between A and S through the row clustering variables Zs. More precisely, we consider that species of row block 1 disappear first, then block 2, etc. Formally, let B be the uniform probability on the row species following a given ordering of the blocks, then S|Z ∼ B if:

P(S = s|Z) = 1 Z s(1) ≤•••≤Z s(nr ) #{s : Z s(1) ≤ • • • ≤ Z s(nr) } . (11)
The support of B is restricted to Qr k=1 n k ! elements, where n k is the cardinal of block k. Equations ( 7), ( 8), ( 9) and ( 11) define a joint probability distribution on (A, S) such that A and S are independent conditionally to Z. We denote L θ,n,B this joint distribution on (A, S) when A ∼ biSBM (θ, n). Under this joint distribution, if θ is such that δ k+ = Qc q=1 ρ q δ kq is a decreasing (resp. increasing) sequence, then B will generate sequences with an expected decreasing (resp. increasing) sequence of degrees, i.e. with decreasing (resp. decreasing) connectivity. In the following, these extinction sequences will be referred to as block-decreasing (resp. block-increasing). The blocks may also be used to generate extinctions by ecological traits that correspond to the blocks.

If A ∼ DCbiSBM (θ, γ, n), the extinction sequences are not related to expected degrees anymore but it may still make sense to generate extinction sequences linked to ecological traits.

Moments of the robustness statistic

Studying the distribution of R(A, S, m) or R(A, S) = 1 nr M m=1 R(A, S, m) under the joint distribution of A and S could be done by simulation, using the fact that the biSBM and DCbiSBM are generative models. However, this comes at a computational cost. In this section, we prove that, when (A, S) ∼ L θ,n,U or (A, S) ∼ L θ,n,B the first moments of the robustness are tractable in a closed-form. The proof relies on the exchangeability of the nodes under a biSBM and on the fact that the considered extinction sequences are adapted to this exchangeability. When A follows a DCbiSBM, the nodes are no longer exchangeable because of the parameters γ associated with each node. Although we are able to obtain a closed-form expression, it is not tractable. Therefore, we will rely on a Monte Carlo approximation for summing on the extinction sequences.

We now exhibit explicit expressions

of R θ,n,S (m) = E A [R S (A, m)] when (A, S) ∼ L θ,n,U or L θ,n,B and of V A [R U (A, m))] when A ∼ biSBM (θ, n).

Expectation

Proposition 1. Let (A, S) ∼ L θ,n,U and set δ +q = Qr k=1 π k δ kq . Then,

• ∀m = 1, . . . , n r : R θ,n,U (m) = 1 - Qc q=1 ρ q (1 -δ +q ) nr-m , (12)
• Consequently, the robustness statistic is:

R θ,n,U = 1 n r nr m=0 R θ,n,U (m) = 1 - 1 n r Qc q=1 ρ q (1 -δ +q ) -(1 -δ +q ) nr+1 δ +q . ( 13
) Proof. Since (A, S) ∼ L θ,n,U , R A,S (m) = E A [E S [ R(A, S, m) ]] = E A [ R U (A, m) ] = 1 n r ! s∈Sn r 1 - 1 n c nc j=1 E A 1 nr i=m+1 A s(i)j =0
.

Using the property of exchangeability of the biSBM, we have that, for any permutation s,

E A 1 nr i=m+1 A s(i)j =0 = E A 1 nr i=m+1 A ij =0
. Moreover, introducing the latent variables Z and W we have:

E A 1 nr i=m+1 A ij =0 = k m+1:nr ∈{1,...,Qr} nr -m Qc q=1 P nr i=m+1 A ij = 0|Z m+1:nr = k m+1:nr , W j = q × P (Z m+1:nr = k m+1:nr , W j = q) = k m+1:nr ∈{1,...,Qr} nr -m Qc q=1 nr i=m+1 (1 -δ k i q ) nr i=m+1 π k i ρ q = Qc q=1 ρ q k m+1:nr ∈{1,...,Qr} nr -m nr i=m+1 (1 -δ k i q )π k i = Qc q=1 ρ q Qr k=1 π k (1 -δ kq ) nr-m .
As a consequence,

R θ,n,U (m) = 1 - Qc q=1 ρ q Qr k=1 π k (1 -δ kq ) nr-m = 1 - Qc q=1 ρ q (1 -δ +q ) nr-m
where δ +q = Qr k=1 π k δ kq . Then, averaging over m leads to:

R θ,n,U = 1 n r nr m=0 R θ,n,U (m) = 1 - 1 n r Qc q=1 ρ q (1 -δ +q ) -(1 -δ +q ) nr+1 δ +q .
Note that, if the network has no specific structure (

δ kq = d or Q r = Q c = 1) then R θ,n,U (m) = 1 -(1 -d) nr-m .. Proposition 2. Let (A, S) ∼ L θ,n,B , then R θ,n,B (m) = 1 - Qc q=1 ρ q n 1 +•••+n Qr =nr n r ! n 1 ! . . . n Qr ! Qr k=1 π n k k (1 -δ kq ) min + (n k , l≤k n l -m) , (14)
where min + is the positive part of the minimum function: min + (x, y) = max(0, min(x, y)).

The proof of Proposition 2 is provided in the supplementary material S-1. The robustness statistic R θ,n,B is the mean of the R θ,n,B (m)'s: no simplified expression has been obtained. Note that, if

Q r = 1 or if δ kq = d, then R θ,n,B (m) = R θ,n,U (m).
In Equation ( 14), the summation over the partitions of the n r row species into Q r blocks may be burdensome if n r or Q r are large. In such cases, a Monte Carlo approximation may be used. Among the ecological networks we consider in Section 6, only a few of them require this approximation.

Variance

We now aim at computing the variance of the robustness V A [R U (A, m)], when A ∼ biSBM (θ, n), thus quantifying the variability of the robustness for a sample of networks sharing the same block models parameters (i.e. the same mesoscale patterns of connectivity).

Proposition 3. Let A ∼ biSBM (θ, n), η q = 1 -δ +q and η qq = Qr k=1 π k (1 -δ kq )(1 -δ kq ).
1. Then:

V A [R U (A, m)] = 1 n c min(2m,nr) l=m m l-m nr-m l-m nr m Qr q=1 ρ q η l q -( Qr q=1 ρ q η m q ) 2 + (n c -1) n c max(2m,nr) l=m m l-m nr-m l-m nr m Qr q,q =1 ρ q ρ q (η q η q ) l-m η 2m-l qq .
2. The variance of the robustness statistic due to the network variability under a given biSBM is:

V A [R U (A)] = 1 n 2 r n 2 c n c nr m,m =0 min(m+m ,nr) l=max(m,m ) m l-m nr-m l-m nr m Qr q=1 ρ q η l q -( 1 n r nr m=0 Qr q=1 ρ q η m q ) 2 + 1 n 2 r n 2 c n c (n c -1) nr m,m =0 min(m+m ,nr) l=max(m,m ) m l-m nr-m l-m nr m Qr q,q =1 ρ q ρ q η l-m q η l-m q η m+m -l qq .
The proof is provided in the supplementary material S-1. The expectations of the robustness under the biSBM given in Proposition 1 and 2 only rely on the model parameters θ and the number of rows n r , whereas the expression of the variance also involves the number of columns n c . Note that for the block sequences S = B, the calculus is not so simple, but the value could still be estimated by simulations if required.

V A (R U (A)) quantifies the variability of the robustness statistic among a population of networks distributed as biSBM. However, from an ecological point of view, it could also be interesting to quantify the variability of the robustness of a network with respect not only to the network but also to the extinction sequence. This is obtained by computing

V A,S 1 n r nr m=1 R(A, S, m) = V A,S R(A, S)
Remark 3. This total variance can be reformulated as:

V A,S R(A, S) = E S (V A|S R(A, S) ) + V S (E A|S (R(A, S) )) = E S (V A R(A, S) ) (15) = E A (V S R(A, S) ) + V A (E S (R(A, S) )) (16)
because S and A are independent under L θ,n,U and E A|S (R(A, S)) does not depend on S by exchangeability. The total variance can be expressed explicitly for S = U as

V A,S R(A, S) = n c n 2 r n 2 c nr m,m =0 Qr q=1 ρ q η nr-min(m,m ) q - 1 n r Qc q=1 ρ q η q -η nr+1 q δ +q 2 + n c (n c -1) n 2 r n 2 c nr m,m =1 Qr q,q =1 ρ q ρ q η max(m,m )-min(m,m ) q η nr-max(m,m ) qq
, through the computation of Equation (15). The second term in Equation ( 16) is provided in Proposition 3. Thus we are able to compute the remaining terms to understand the various sources of variability (due to S or A).

Illustration of the variability of the robustness function

In Figure 1 This corresponds to a so-called core-periphery structure where the first blocks (in rows and columns) are well connected (the core) and the second blocks are less connected (the periphery). We represent the functions m → R θ,n,S (m) with S = U and S = B such that the blocks are ordered by decreasing or increasing connectivity (solid black lines in Figure 1). For the uniform distribution, we also plot the area given by m

→ R θ,n,U (m) ± 2 V A (R U (A, m)) (grey ribbon).
The colored dotted lines are Monte Carlo estimates of the robustness functions R S (A, m) corresponding to 10 simulated networks for the same extinction sequence distributions. The Monte Carlo estimates are computed over 300 realizations of extinction sequences. The inflection points and the dispersion of the dotted lines observed respectively around an extinction rate of 0.25 for decreasing extinction and 0.75 for increasing extinction correspond respectively to the proportion of the core or the periphery blocks. When the rate of extinction exceeds one of these proportions, the extinctions which were only happening in one block then happen in the other block. Also notice that when S = B ordered by increasing connectivity, some networks still have a robustness of 1 even after a large fraction of primary extinctions has occurred. The robustness for this primary extinction distribution is highly dependent on the degree of the most connected primary species.

Figure 1: Robustness function computed from a set of biSBM parameters (plain black) and by Monte Carlo for 10 networks generated from the same biSBM distribution (dotted) for block decreasing, uniform and block increasing primary extinction sequences. The grey ribbon on the uniform facet is twice the standard error given in Proposition 3.

Impact of the Network Structure on the Robustness

From the expressions derived in the previous section when A ∼ biSBM (θ, n), we are now able to study the average behavior of the robustness with respect to the mesoscale structure of the network encoded in θ.

Analytical Properties

5.1.1 First properties on R θ,n,S (m)
We first derive the following straightforward but useful properties of the robustness function and statistic:

Properties 1.

1. Under the joint distribution L θ,n,S where S ∈ {U, B}, the following properties hold:

(a) the function m ∈ {0, . . . , n r } → R θ,n,S (m) is a strictly decreasing function provided that δ k+ > 0 for all k > 0,

(b) R θ,n,S (0) ≤ 1 -(1 -d) nr ≤ 1 where d = k,q π k ρ q δ kq , (c) For θ = (π, ρ, δ) and θ = (π , ρ , δ ) such that π = π , ρ = ρ and ∀(k, q) ∈ {1, . . . , Q r } × {1, . . . , Q c } δ kq ≤ δ kq , we have ∀m ∈ 0, . . . , n r , R θ,n,S (m) ≤ R θ ,n,S (m) and R θ,n,S ≤ R θ,n,S .
2. Under the joint distribution L θ,n,U , if n r ≤ n r then R θ,(nr,nc),U < R θ,(n r ,nc),U .

Proof 12) and ( 14).

For Property 2., first notice that R θ,(nr+1,nc),U (m + 1) = R θ,(nr,nc),U (m). Then,

R θ,(nr+1,nc),U = 1 n r + 1 R θ,(nr+1,nc),U (0) + nr m=0 R θ,(nr+1,nc),U (m + 1) = 1 n r + 1 R θ,(nr+1,nc),U (0) + nr m=0 R θ,(nr,nc),U (m) = 1 n r + 1 R θ,(nr+1,nc),U (0) + n r R θ,(nr,nc),U (Property 1.(a)) > 1 n r + 1 R θ,(nr,nc),U (0) + n r R θ,(nr,nc),U > 1 n r + 1 R θ,(nr,nc),U + n r R θ,(nr,nc),U = R θ,(nr,nc),U
Properties 1.(a) and 1.(b) supply a bound depending on the expected density and the number of row species n r . Property 1.(d) implies that for a given block structure, the robustness increases when each element δ kq increases while Property 2. states that the uniform robustness automatically increases with the size n r of the network. However, it is important to note that these properties do not assess that the robustness is an increasing function of the connectance d.

Upper bound for the robustness under a uniform extinction sequence

We now aim at identifying mesoscale structures maximizing the average robustness under the joint distribution L θ,n,U . In order to remove the effect of the mean number of interactions a.k.a. the density d = k,q π k δ kq ρ q , we propose to compare structures encoded in θ leading to the same density value. We also fix the number of row species n r . Thus, for a given density d, we define the set Θ d = {θ = (π, ρ, δ) : k,q π k δ kq ρ q = d}. The following proposition provides an upper bound for the expectation of the robustness function and of the robustness statistic as a function of the density. It also identifies a condition on θ to achieve this upper bound. Eventually, it shows that the parametrization of an Erdős-Rényi distribution satisfies this condition and reaches the lowest variance of the robustness statistic among the parametrizations satisfying this condition.

Proposition 4. Upper bound of robustness under L θ,n,U .

1. For all m ∈ {0 . . . , n r }: R θ,(nr,nc),U (m

) ≤ 1 -(1 -d) nr-m .
2. For all m ∈ {0 . . . , n r -2} arg max

θ∈Θ d R θ,(nr,nc),U (m) = θ : Qr k=1 π k δ kq = Qr k=1 π k δ kq , ∀(q, q ) ∈ {1, . . . , Q c } 2 := Θ max d,nr ,
Moreover, ∀θ ∈ Θ d and ∀n c :

R θ,(nr,nc),U ≤ 1 - 1 n r (1 -d) -(1 -d) nr+1 d . ( 17
)
3. Assume that E is a network with n r rows and n c columns following an Erdős-Rényi distribution with density parameter d. Then:

min A∼biSBM (θ,n):θ∈Θ max d,nr V A (R U (A)) = V E (R U (E)) .
In other words, among all parameters θ ∈ Θ max d,nr , the one of the Erdős-Rényi network minimizes the variance of robustness statistic R U (A) given in Proposition 3.

Proof.

1. Recall that R θ,(nr,nc),U (m) = 1 -Qc q=1 ρ q (1 -δ +q ) nr-m where δ +q = Qr k=1 π k δ kq for any q ∈ {1, . . . , Q r }. So, R θ,(nr,nc),U (n r ) = 0 and for R θ,(nr,nc),U (n r -1) = 1 -(1 -d) : the bound is true for m = n r and m = n r -1. Now, if 0 ≤ m ≤ n r -2, then x → x nr-m is a strictly convex function and the Jensen. inequality applies:

R θ,(nr,nc),U (m) = 1 - Qc q=1 ρ q (1 -δ +q ) nr-m ≤ 1 - Qr q=1 ρ q (1 -δ +q ) nr-m (18) = 1 -1 - Qr q=1 ρ q δ +q nr-m = 1 -(1 -d) nr-m ,
2. Equality at line (18) holds if and only if the term inside the strictly convex function is constant, ie. for any q, q , δ +q = δ +q .

3. When we compute the variance given in Proposition 3(2) on any biSBM with θ ∈ Θ max d,nr , we notice that η q = η q = 1 -d and that only the term q,q ρ q ρ q η l-m q η l-m q η m+m -l qq varies. We can reformulate this quantity and use the Jensen inequality:

q,q ρ q ρ q η l-m q η l-m q η m+m -l qq = (1 -d) l-m (1 -d) l-m q,q ρ q ρ q η m+m -l qq (Jensen) ≥ (1 -d) 2l-m -m q,q ρ q ρ q k π k (1 -δ kq )(1 -δ kq ) m+m -l (Distributivity) = (1 -d) 2l-m -m k π k (1 - q ρ q δ kq )(1 - q ρ q δ kq ) m+m -l (Jensen) ≥ (1 -d) 2l-m -m k π k (1 - q ρ q δ kq ) 2 m+m -l = (1 -d) m +m (Initial term from ER)
Although the homogeneous distribution on the networks (Erdős-Rényi) leads to the maximum robustness in expectation for a given density, a particular realization of a network according to another distribution may be likely to have a larger robustness than a realization according to the Erdős-Rényi distribution. Indeed, the variance is larger when it corresponds to a distribution that represents a more complex structure. This behavior is illustrated in Figure 2.

Analysis for Typical Structures

We now illustrate numerically how the robustness statistic varies with respect to the network topology2 . For that purpose, we fix n r = n c = 100, Q c = Q r = 2, and π = (.25, .75) and make ρ vary. For j ∈ [1/8, 8], we consider the following connectivity matrices.

Modular: δ = j 1 1 j , each block of row species is strongly connected to a block of column species and lowly connected to the other.

Core-periphery: δ = j j j 1 .

• For j > 1, the structure is nested, the core is strongly connected to the whole network while the periphery is lowly connected with itself.

Figure 2: Density curve of the distribution of the expectation of the robustness statistic under biSBMs of the same size and density with three different sets of parameters. Each network is simulated from a given set of parameters, then the robustness is computed by the Monte-Carlo approximation given in (4). 500 simulations each.

• For j < 1, the core is strongly connected with itself but the rest of the network is lowly connected.

Each connectivity matrix is then normalized such that the density of the network is equal to 0.0156, by applying the following transformation: δkq = 0.0156 δ kq k ,q π k δ k q ρ q ∀k, q ∈ {1, 2}. This value is chosen so the robustness statistic associated with an Erdős-Rényi distribution with density d = 0.0156 and n r = 100 rows is approximately equal to 0.5. We recall this value is an upper bound of the expectation of the robustness statistic under the joint distribution L θ,n,U where θ is such that the network has the same density and same number of rows.

The extinction sequence distributions are S ∈ {U, B ↑ , B ↓ } where B ↑ (resp. B ↓ ) corresponds to the block increasing (resp. decreasing) extinction sequences distribution defined in Equation ( 11). For every S, topology (modular or core-periphery) and value of j and ρ, we compute the expectations of the robustness statistics by using the expressions derived in Section 4. These expectations are displayed in the heat maps of Figure 3.

For modular networks, the plots are symmetric with respect to the central dot which corresponds to the case of an Erdős-Rényi distribution. The impact of the modular structure is slighter for S = U than for S ∈ {B ↑ , B ↓ }. For S = U, the strongest impact is observed when the modular structure is strong and the most connected blocks are slightly larger than the least connected ones (Fig. 3.♦). Keeping the same strong modular structure with the most connected blocks slightly smaller than the least connected ones (Fig. 3. ) leads to a negative Figure 3: Robustness for different biSBM topologies. Primary extinction sequences are displayed in rows and topologies in columns. In abscissa, j is a topology strength parameter. rho is the proportion of column species that belongs to the first column block. • represents a topology with no structure (Erdős-Rényi (ER)). Color gradient varies from blue (less robust than an ER), to white (as robust as an ER), to red (more robust than an ER). The symbols corresponds to the following mesoscale structures:

-large core, strongly connected only to itself, -small core, strongly connected to the whole network, -medium size core strongly connected only to itself, -highly modular with unbalanced blocks' sizes, ♦ -highly modular with slightly smaller highly connected blocks, -highly modular with slightly larger highly connected blocks.

impact on the robustness no matter S.

When the network is highly modular and the most connected block is small (j and ρ both small or both large, Fig. 3. ) the robustness is strongly impacted. This impact is negative for B ↓ and positive with B ↑ .

For a core-periphery structure, there is a clear asymmetry between j < 1 and j > 1. The robustness statistic tends to be smaller when the core is mainly connected to itself, especially when the core is small (j < 1 and large ρ: Fig. 3. ) no matter S. On the other hand, coreperiphery structure with small core that are highly connected to the whole network have the strongest impact on the robustness, negatively when S = B ↓ and positively when S = B ↑ (Fig. 3. ). The effect of this structure ( ) is very slight on uniform extinction sequences whereas the effect tends to get larger when the blocks' sizes are more balanced (Fig. 3. ).

Analysis of a collection of observed bipartite ecological networks

In this section, we analyse the robustness of a collection of 136 plant-pollinator, seed-dispersal or host-parasite networks, issued from the Web of Life dataset (www.web-of-life.es). The selected networks involve at least 10 row species and 10 column species.

When observing an interaction network, one can compute its empirical robustness as defined in Equation ( 3). However, understanding the behavior of the robustness statistic for any network with the same probabilistic distribution (i.e. with the same mesoscale structure) may also be attractive and informative. This can be done by positing a model on the network of interest, estimating the corresponding parameters (summarizing its structure) and deriving the moments of the robustness for these inferred parameters using Section 4 (either using a closedform expression or a Monte Carlo integration, depending on the model). This expected version of the standard robustness may be considered as a new robustness indicator.

In what follows, we put in perspective the empirical robustness and its expected versions for the biSBM and DCbiSBM models (Subsection 6.1) and comment the systematic differences we observe. Then, we demonstrate the interest of the expected version when the network is partially observed (Subsection 6.2). Indeed, when inferring the bipartite SBM, the observational process that generates the observed network with possibly missing data could be taken into account [START_REF] Tabouy | Variational inference for stochastic block models from sampled data[END_REF]. This which allow us to compensate for observational biases in the empirical robustness.

Inferring the biSBM and DCbiSBM The parameters of these two models can be inferred by a variational version of the Expectation-Maximization (EM) algorithm. The number of blocks is chosen according to an Integrated Classification Likelihood (ICL) criterion [START_REF] Daudin | A mixture model for random graphs[END_REF]. This variational EM algorithm is theoretically grounded [START_REF] Bickel | Asymptotic normality of maximum likelihood and its variational approximation for stochastic blockmodels[END_REF] and has proven its pratical efficiency [START_REF] Daudin | A mixture model for random graphs[END_REF][START_REF] Mariadassou | Uncovering latent structure in valued graphs: a variational approach[END_REF]. In practice for the inference of these models, we use the blockmodels R package [START_REF] Leger | blockmodels: Latent and Stochastic Block Model Estimation by a 'V-EM[END_REF] and to handle missing observations the GREMLINS R package [START_REF] Bar-Hen | Block models for generalized multipartite networks: Applications in ecology and ethnobiology[END_REF], both available on CRAN .

Computation of the Robustness for the Web of Life Dataset

For each network A, on the one hand, we compute the empirical robustness R U (A) using 300 Monte Carlo realisations. On the other hand, for each model (biSBM and DCbiSBM) we denote θ (resp. θ, γ) the estimated parameters and compute the expectations of the robustness under each model. For the biSBM, we also supply their variance and define the ratio

Z R (A) =| R U (A) -R θ,n,U | / V θ,n ( E U [ R(A, S)|A ] )
which is the number of standard deviations separating the two computed robustness statistics. This quantity helps assess the goodness of fit of the biSBM to the network A with respect to the robustness statistic. biSBM versus DCbiSBM In Figure 4, we plot the points ( R U (A), R θ,n,U ) for the two models (circles are for biSBM and purple triangles for DCbiSBM). For the biSBM, the color of the point depends on (1-d) nr which is the estimated probability to observe a column with no interaction. The observed robustness indicators range from 0.5 to 1. The points roughly follow the identity line for the biSBM (circles) whereas the expected version seems systematically smaller than the empirical robustness for the DCbiSBM (purple triangles). We comment further on these points hereafter.

• Under a biSBM, the smaller (1 -d) nr , the smaller | R U (A) -R θ,n,U | and Z R (A) (the number of standard deviations). When (1 -d) nr is large (red dots), the expected version of the robustness under a biSBM underestimates the standard robustness (red dots at the bottom left of the plot). This phenomenon may be explained as follows. By construction, ecological networks only involve species that have been seen at least one time in interaction (species with no interaction were removed). As a consequence, since the biSBM does not take into account this phenomenon, the space of networks we integrate over when computing the robustness statistic under a biSBM is too large. This remark is especially true for small networks and highlights the limitation of the biSBM model for small ecological networks.

• On 127 networks (over the 136), the DCbiSBM selects 1 block, encoding solely, as the EDD model, the degree distribution and the size of the networks in their parameters. On these networks, the expectation of the robustness statistic under a DCbiSBM is smaller than the empirical robustness most of the time (purple triangles). We can conclude that the DCbiSBM seems to be unable to capture the additional structure beyond their degree distribution that makes them more robust to uniform random extinction than expected. This first study highlights the limitation of the DCbiSBM model to mimic ecological networks (hence, hereafter, we focus our analysis solely on the biSBM). On the contrary, for not too small networks, the expected version of the robustness with biSBM supplies coherent robustness values with the empirical indicator R U (A). This new version of the robustness has the advantage to arise in a closed-form and a quantification of its variance is available. This comparison has been made for S = U, we now consider more elaborate extinction sequence distributions.

About S.

We compute the robustness statistics for decreasing and increasing primary extinction sequences on the degrees with the two methods described in Section 2: the one which strictly depends on the order of the degrees ( R

D ↑ ord (A) and R D ↓ ord (A)
) and the one where the nodes are weighted by a linear function of the degrees as in Equation ( 5) with α = 1 ( R

D ↑ lin (A) and R D ↓ lin (A))
. We examine the fit with R θ,n,S where S ∈ {B ↓ , B ↑ } mimics the degree increasing and decreasing extinction sequences. The comparison are plotted in Figure 5.

When considering primary extinction sequences by decreasing connectivity, there is a pos-

itive bias of R θ,n,B ↓ on R D ↓ ord (A) (top-left) which is attenuated when considering R D ↓ lin (A) (bottom-left).
For primary extinction sequences by increasing connectivity, the empirical robustness is highly dependent on the degrees of the most connected species, hence the fit of R θ,n,B ↑ on R D ↑ ord (A) is very poor with a negative bias (top-right). While the negative bias still remains, the fit on R D ↑ lin (A) is correct for a much higher fraction of the networks (bottom-right).

As a conclusion, the degree decreasing sequences standardly used in the ecological field can be easily reproduced with our block decreasing connectivity sequences, leading to quite comparable values. Once again, our expected version of the robustness can be calculated in a closed-form while the empirical robustness relies on a computationally expensive Monte Carlo integration. The degree decreasing extinction sequences are very sensible to the most highly connected species leading to no agreement between the two versions of the robustness.

Correction for Partially Observed Networks

Although the ecological networks are often considered to describe all the possible interaction between species, the sampling may be incomplete [START_REF] Blüthgen | What do interaction network metrics tell us about specialization and biological traits[END_REF] which may bias the computed network statistics [START_REF] Rivera-Hutinel | Effects of sampling completeness on the structure of plant-pollinator networks[END_REF] such as the robustness. By relying on a probabilistic model such as the biSBM which can be adjusted to account for the observation process, the sampling effect on the robustness can be corrected. We assume that we could have obtained the true interaction network if the sampling effort has been large enough. Instead of this true network, we have only a partially observed interaction network that corresponds to a subset of the true network. We assume that one of the two following frameworks may have generated missing data:

Partially observed species 25% of both the row species and the column species are removed uniformly at random, resulting in the following networks subset: A obs = {A ij : i and j are observed}. In this framework, we need to assume that some other data or some expert knowledge gives us the true number of species and the density of the true network. By relying on this expert knowledge, we are able to adjust the parameters of the biSBM.

Partially observed interactions

The observation process consists of the observation of interaction on two different transects T 1 and T 2 . Since not all the species are present on both transects, the interactions between the species which were not observed on the same transect are labelled as missing and encoded by NAs. This results in the following modified incidence matrix:

A obs =       j ∈ {T 1 ∩ T 2 } j ∈ {T 1 \ T 2 } j ∈ {T 2 \ T 1 } i ∈ {T 1 ∩ T 2 } A ij A ij A ij i ∈ {T 1 \ T 2 } A ij A ij NA i ∈ {T 2 \ T 1 } A ij NA A ij      
.

More precisely, we consider that on average 50% of the species were observed on both transects while 25% were just observed on one of the two transects and select those species uniformly at random, resulting in 12.5% of missing interactions on average. Note that in this case, we do not need any expert or additional knowledge to have an unbiased estimation of the biSBM parameters. It is sufficient that the missing values are taken into account in the biSBM inference.

We performed a simulation study where we considered 98 networks from the web of life dataset as the true interaction networks. We kept the 98 networks where (1 -d) nr < .1 and the error in terms of standard deviation is smaller than 1. For each of these 98 networks, we simulated the two observation processes described above 30 times, then we computed the standard robustness statistic and its expectation under a biSBM. For Figure 6, we computed the root mean squared error (RMSE) for each network and for different cases (in terms of observation process and extinction sequence) between the robustness computations on the true network (fully observed) and on the partially observed network. When the method for computing the robustness is based on MC computations, we considered the uniform primary extinction sequences (named Monte Carlo in Fig. 6) and the increasing and decreasing extinction sequences that depend strictly (named Ordered Monte Carlo in Fig. 6) or proportionally as in Equation (5) (named Linear Monte Carlo in Fig. 6) on the row degrees. More specifically, the RMSE for the methods based on MC computations are computed for any extinction sequences for each of the 98 fully observed networks A and where the Ãb 's are realizations of a partial observation of A. In order to get the robustness under a biSBM, the biSBM was first inferred by taking into account missing data and the robustness statistics were then computed under the inferred biSBM. We computed the corresponding RMSE for the distributions on extinction sequences S ∈ {U, B ↑ , B ↓ } as:

1 30 30 b=1 (R θ,n,S -R θb ,n,S ) 2
where θ is estimated from the fully observed network and the θb 's are estimated from partial observations of the network. Note that the numbers of species in row and in column are the same for the fully and the partially observed networks even in the case of missing species since we assumed that we had some additional knowledge which provided us with this information.

The errors in the prediction of the robustness statistics computed from partially observed networks are much smaller when computing the robustness under a biSBM than when using MC computations. Indeed, the missing data can be accounted for in the biSBM inference and the underlying structure of the network can still be recovered from partial information whereas the Monte Carlo simulations are more sensitive to perturbations in the network. The extinction sequences which depend strictly on the degrees are the most impacted by a partial observation of the network. This impact is rather strong although only 12.5% of the information is missing. Note that adding some randomness in the primary extinction sequences by using a distribution which depends linearly on the degrees instead of a strict order has a stabilizing effect.

In the framework with partially observed species we assumed that we had additional knowledge giving us access to the true number of species and network density. In practice, this is not always the case. When details about the sampling of networks are available, a few methods exist in the literature to estimate the number of species [START_REF] Jiménez-Valverde | Comparing relative model fit of several species-accumulation functions to local Papilionoidea and Hesperioidea butterfly inventories of Mediterranean habitats[END_REF][START_REF] Gotelli | Estimating species richness[END_REF] and similar methods could be used to estimate the density. One could then plug these estimates in the formula of robustness. If neither the sampling scheme nor the expert knowledge is available, the true number of species remains unknown. Note that in most available ecological interaction networks, the number of species is underestimated as species which have not been seen in interaction with other species are not included in the network. Finally, having access to the number of species, but not to the true number of interactions,

Discussion

We proposed an expected version of the empirical robustness for bipartite ecological interaction networks by considering a joint model on the network and on the primary extinction sequences. In particular, we obtained a closed and tractable form of the robustness when considering a biSBM as the network model for uniform and by blocks primary extinction sequences. We validated our method by showing that the obtained values were consistent with the empirical robustness classically used in the ecological network community. Having analytical forms allows us to better understand the impact on the robustness of the topology of the network in terms of number of species, density, and mesoscale connectivity patterns. Furthermore, we used the difference between the empirical and expected versions to show that the biSBM is better suited than the DCbiSBM for bipartite ecological networks. This could come from the fact that the DCbiSBM tends to favor community structure above core-periphery structure (Newman, 2018, 14.7.3). The core-periphery structure is indeed very common in bipartite ecological networks which has a strong impact on the robustness.

On real networks, this method can serve as an alternative to the traditional empirical approach, especially when networks are partially observed (because of the inference stability of the biSBM) or when the sampling effort is incomplete (because the model parameters can be easily corrected). This step could be improved with more precise information by considering specific sampling schemes in the model [START_REF] Tabouy | Variational inference for stochastic block models from sampled data[END_REF] or by obtaining more details on the sampling in order to better estimate the parameters block by block. Moreover, from the observation of a network, the impact of some hypotheses on the structure might be tested by tuning some parameters of the models and computing the corresponding expected robustness. It may also help study the impact of the structure beyond the effect of the number of species or density when comparing the robustness of several networks: indeed, once the biSBM parameters are estimated, the robustness statistics can be computed by setting the same numbers of species and the same density for all the networks.

Although we developed our method for bipartite networks, this approach can be extended to other types of networks in particular multipartite networks [START_REF] Pocock | The robustness and restoration of a network of ecological networks[END_REF] and food webs. Multilayer networks (multiplex or multipartite networks) are gaining a lot of attention with scholars studying ecological networks [START_REF] Hutchinson | Seeing the forest for the trees: Putting multilayer networks to work for community ecology[END_REF]. Extending this framework to multipartite networks by including cascading effect between layers is an interesting perspective once data will be more readily available. The study of robustness for food webs is quite active but requires some ecological insight to properly model the food web as some species are basal species and thus do not prey on other species (in-degree is 0). So food webs which are usually blockmodeled by a directed stochastic block model, might be better handled in our case using a multipartite block model with basal species as a functional group of its own. In this case, it might be important to incorporate rewiring or cascading mechanisms to the modeling of the extinctions. A direction to look at in order to deal with the extinction of basal species and the incorporation of cascading mechanisms is the approach of Bayesian networks to model species extinction in food webs [START_REF] Eklöf | Secondary extinctions in food webs: a bayesian network approach[END_REF][START_REF] Häussler | A bayesian network approach to trophic metacommunities shows that habitat loss accelerates top species extinctions[END_REF].

Lastly, we believe that other ecological indicators could be estimated through a parametric model based approach. Especially, the EDD model and the DCbiSBM seem particularly well suited to study nestedness (see [START_REF] Mariani | Nestedness in complex networks: observation, emergence, and implications[END_REF], for a review) which is a widely used statistic for ecological networks.

S-1 Proof for Section 4 (Moments of the robustness statis

- tic) Proposition 2. Let (A, S) ∼ L θ,n,B , then R θ,n,B (m) = 1 - Qc q=1 ρ q n 1 +•••+n Qr =nr n r ! n 1 ! . . . n Qr ! Qr k=1 π n k k (1 -δ kq ) min + (n k , l≤k n l -m) , (S-19)
where min + is the positive part of the minimum function: min + (x, y) = max(0, min(x, y)).

Proof. Let us assume without loss of generality that the primary extinction sequences go from block 1 to block Q r . Under this assumption, P(S = s|Z) =

1 Z s(1) ≤•••≤Z s(nr )
#{s:Z s(1) ≤•••≤Z s(nr ) } . Thus, an extinction sequence which does not maintain the block ordering has a null probability. Then, conditioning first by the blocks memberships then the primary extinction sequences:

E (A,S) [R(A, S, m)] = Qc q 1:nc ρ 1:nc Qr k 1:nr π 1:nr E (A,S) [R(A, S, m)|Z 1:nr = k 1:nr , W 1:nc = q 1:nc ] (S-20) = 1 - 1 n c nc j=1 Qc q ρ q Qr k 1:nr π 1:nr E (A,S) nr i=m+1 (1 -A S(i)j )|Z 1:nr = k 1:nr , W j = q = 1 - 1 n c nc j=1 Qc q ρ q Qr k 1:nr π 1:nr s∈Sn r P(S = s | Z 1:nr = k 1:nr ) nr i=m+1 (1 -δ k s(i) q ) .
We have for all sequences s, s in the support of S | Z, that for all i, Z s(i) = Z s (i) . Hence, using the exchangeability properties of the biSBM for species belonging to the same block:

nr i=m+1 1 -δ k s(i) q = nr i=m+1 1 -δ k s (i) q . Furthermore, nr i=m+1 1 -δ k s(i) q = nr i=m+1 Qr k=1 (1 -δ kq ) 1 {k s(i) =k} = Qr k=1 (1 -δ kq ) nr i=m+1 1 {k s(i) =k} , with for all s such that k s(1) ≤ • • • ≤ k s(nr) : nr i=m+1 1 {k s(i) =k} =                0 if l≤k n l ≤ m, l≤k n l -m if l≤k n l -n k < m < l≤k n l , n k if m ≤ l≤k n l -n k .
This leads to: Proposition 3. Let A ∼ biSBM (θ, n), η q = 1 -δ +q and η qq = Qr k=1 π k (1 -δ kq )(1 -δ kq ).

1. Then: ρ q ρ q (η q η q ) l-m η 2m-l qq .

V A [R U (A, m)] =
2. The variance of the robustness statistic due to the network variability under a given biSBM is: ρ q ρ q η l-m q η l-m q η m+m -l qq .

V A [R U (A)] = 1 n 2 r n 2
Proof. The calculus of the following proof relies on the following lemma:

Combinatorial lemma 1. Let m, m the first terms of two permutations of S n , the the proportion of couple of permutation (s, s ) that have exactly l unique terms is: Proof. There are n! permutations of size n. For the first permutation of size n, we look at the first m . In order to create the second permutation we must among the first m terms take:

• l -m terms that are not common with the first permutation (among n -m)

• m + m -l terms that are common with the first permutation (among m). The bound on l is straightforward.

We first prove the result for the robustness statistic. The variance of the robustness function is a straightforward derivation from it.

V A ( E S [ R(A, S))|A ] ) = E A [ E 2 S [ 1 -R(A, S)|A ] ] B -E 2 A [ E S [ 1 -R(A, S)|A ] ] C (S-22)
Using Equation ( 12) and Equation (4), we have: C = 1 nr nr m=0 Qc q=1 ρ q η nr-m q 2

. We divide B into 2 terms, based on the column index: (S-23)

To compute B1, we separate the set of the union of the extinction sequences into 3 disjoint sets in order to distribute the mixture parameter π:

B1 = 1 n r !n r ! s,s ∈Sn r E π,ρ   E δ 1 { nr i=m+1 A s(i)j =0} 1 { nr i=m +1 A s (i)j =0} |Z, W   = 1 n r !n r ! s,s ∈Sn r E π,ρ nr i=m+1 (1 -δ Z s(i) W j ) nr i=m +1
(1 -δ Z s (i) W j ) = 1 n r !n r ! s,s ∈Sn r q,q ρ q ρ q k 1:nr π k 1:nr nr i=m+1

(1 -δ k s(i) q ) nr i=m +1

(1 -δ k s (i) q ) = 1 n r !n r ! s,s ∈Sn r q,q ρ q ρ q k 1:nr π k 1:nr i∈s(m+1:nr)\s (m +1:nr)

(1 -δ k i q ) i∈s (m +1:nr)\s(m+1:nr)

(1 -δ k i q ) i∈s(m+1:nr)∩s (m +1:nr)

(1 -δ k i q )(1 -δ k i q ) (Lem.1) = l m m+m -l n-m l-m n m
q,q ρ q ρ q η l-(nr-m ) q η l-(n-m) q η 2n-m-m -l qq For B2, the column indices are the same, hence some entries of A are taken twice:

B2 = 1 n r !n r ! s,s ∈Sn r q ρ q k 1:n π k 1:n E δ 1 { i∈s(m+1:nr ) A ij =0}
•1 { i ∈s (m +1:nr ) A i j =0} |Z 1:nr = k 1:nr , W j = q = 1 n r !n r ! s,s ∈Sn r q ρ q k 1:n

  , we illustrate the variations of the robustness function for a given structure encoded in θ = δ = .4 .15 .25 .05 , π = (0.25, 0.75), ρ = (0.2, 0.8) .

.

  Property 1.(a) comes from the robustness definition and from the fact that s(m + 1 : n r ) ⊂ s(m : n r ) for any extinction sequence s. Property 1.(b) is a consequence of Proposition 4. Property 1.(c) is true for all m ∈ {0, . . . , n r } from Equations (

Figure 4 :

 4 Figure 4: Expected robustness statistic (R θ,n,U ) under a biSBM (circles) and DCbiSBM (purple triangles) as a function of the standard robustness ( R U (A)) for uniform extinction sequences. The grey to red gradient stands for (1 -d) nr , the probability to have an empty column under an Erdős-Rényi distribution. The size of the point is related to Z R (A).

Figure 5 :

 5 Figure 5: R θ,n,S as a function of the classical robustness R D (A) for various S. Block Decreasing (resp. increasing) = B ↓ (resp B ↑ ), Ordered Decreasing (resp. Increasing) = D ↓ ord (resp. D ↑ ord ). Linear Decreasing (resp. Increasing) = D ↓ lin (resp. D ↑ lin ). The grey to red gradient stands for (1 -d) nr , the probability to have an empty column under an Erdős-Rényi distribution. The size of the point is the deviation between the biSBM and the empirical robustness in terms of the number of standard deviations of the robustness under a biSBM distribution for S = U.

S

  = {U, D ↑ ord , D ↑ lin , D ↓ ord , D ↓ lin } listed above, as:

Figure 6 :

 6 Figure 6: Error (RMSE) in the prediction of the robustness of 98 fully observed networks computed from partially observed networks. On the left, 12.5% of possible interactions are recoded as NA. and on the right, 25% of species are missing.

  = s |Z 1:nr = k 1:nr ) nr i=m+1 (1 -δ k s(i) q ) = s∈Sn r P(S = s | Z 1:nr = k 1:nr ) Qr k=1 (1 -δ kq ) min + (n k , l≤k n l -m) δ kq ) min + (n k , l≤k n l -m) .(S-21)Inputting Equation (S-21) in Equation (S-20), we obtain the result:E (A,S) [R(A, S, m)] = 1 -Qc q=1 ρ q n 1 +•••+n Qr =nr n r ! n 1 ! . . . n Qr ! × Qr k=1 π n k k (1 -δ kq ) min + (n k , l≤k n l -m) .

  #{(s, s ) : s(1 : m) ∪ s (1 : m ) = l} #(S n , S n ) = max{m, m } ≤ l ≤ min{m + m , n}0 otherwise, where s(1 : m) = {s(1), . . . , s(m)}.

  Those m terms can be reordered into m ! possible arrangements and the n -m resting terms into n -m ! permutations, giving:# {(s, s ) : s(1 : m) ∪ s (1 : m ) = l} = n! m m + m -l n -m l -m m !(n -m )!We then divide by the (n!) 2 couple of permutations possible. #{(s, s ) : s(1 : m) ∪ s (1 : m ) = l} #(S n ,

E

  A 1 -E S [ R(A, S, m)|A] • E S [ 1 -R(A, S, m )|A]

πq

  k 1:n •E δ 1 { i∈s(1:nr -m)∪s (1:nr -m ) A ij =0} |Z 1:nr = k 1:nr , W j = q = 1 n r !n r ! s,s ∈Sn r q ρ q k 1:n π k 1:n i∈s(1:nr-m)∪s (1:nr-m ) (1 -δ k i q )Finally we use the symmetry between m and n r -m and input first B1 and B2 into Equation (S-23) then B and C into Equation (S-22) to obtain the variance of the robustness statistic.To obtain the variance of the robustness function, we fix m, and set m = m in Lemma 1, and Equation (S-23).

Note that the robustness of ecological networks bears a different meaning than the robustness traditionally used in the information or epidemic networks literature. In these fields, the goal is to see how a network stays connected when some nodes are disabled (or on the opposite how epidemics are still able to spread when nodes become immune), the task is tackled by analyzing the size of the different connected components and in particular the existence of a giant connected component (seeNewman,

2018, chapter 15 and reference therein). In ecology, the key issue is to observe how species get isolated from the rest of the ecosystem. This approach is related to the concept of isolated nodes in a graph[START_REF] Erdős | On the evolution of random graphs[END_REF].

All the simulations and estimations in this section are done using the R package robber (Chabert-Liddell, 2021)available on CRAN and documented at https://chabert-liddell.github.io/robber/.

Aknowledgement

The authors would like to thank Sonia Kefi, François Massol and Vincent Miele for their helpful advice. This work was supported by a public grant as part of the Investissement d'avenir project, reference ANR-11-LABX-0056-LMH, LabEx LMH. This work was partially supported by the grant ANR-18-CE02-0010-01 of the French National Research Agency ANR (project EcoNet).


Benchmarking an improved statistical adiabatic channel model for competing inelastic and reactive processes

Maarten Konings, Benjamin Desrousseaux, François Lique, Jérôme Loreau

To cite this version:

Maarten Konings, Benjamin Desrousseaux, François Lique, Jérôme Loreau. Benchmarking an improved statistical adiabatic channel model for competing inelastic and reactive processes. Virtual Winter School on Computational Chemistry 2022, Feb 2022, Leuven, Belgium. hal-03644133

Benchmarking an improved statistical adiabatic channel model for competing inelastic and reactive processes

KU Leuven, Division of Quantum Chemistry and Physical Chemistry, Department of Chemistry, Belgium.

Université de Rennes 1, CNRS, IPR (Institut de Physique de Rennes), France.

E-mail: maarten.konings@kuleuven.be

The computation of state-to-state rate coefficients for complexforming collisional processes using fully quantum mechanical approaches, is challenging. Due to the statistical nature of said collisions, we may resort to much simpler statistical models, such as the statistical adiabatic channel model (SACM). The state-to-state cross sections can be obtained via, where,