Rémi Bonnefoi

Lilian Besson

Christophe Moy

Emilie Kaufmann

Jacques Palicot

Multi-Armed Bandit Learning in IoT Networks: Learning helps even in non-stationary settings

Keywords: Internet of Things, Multi-Armed Bandits, Reinforcement Learning, Cognitive Radio, Non-Stationary Bandits

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Unlicensed bands are more and more used and considered for mobile and LAN communication standards (WiFi, LTE-U), and for Internet of Things (IoT) standards for short-range (ZigBee, Z-Wave, Bluetooth) and long-range (LoRaWAN, SIGFOX, Ingenu, Weightless) communications [START_REF] Centenaro | Long-range communications in unlicensed bands: the rising stars in the IoT and smart city scenarios[END_REF]. This heavy use of unlicensed bands will cause performance drop, and could even compromise IoT promises.

Efficient Medium Access (MAC) policies allow devices to avoid interfering traffic and can significantly reduce the spectrum contention problem in unlicensed bands. As end-devices battery life is a key constraint of IoT networks, this leads to IoT protocols using as low signaling overhead as possible and simple ALOHA-based mechanisms. In this article, we analyze the performance of Multi-Armed Bandits (MAB) algorithms [START_REF] Lai | Asymptotically efficient adaptive allocation rules[END_REF][START_REF] Bubeck | Regret analysis of Stochastic and Non-Stochastic Multi-Armed Bandit Problems[END_REF], used in combination with a timefrequency slotted ALOHA-based protocol. We consider the Upper-Confidence Bound (UCB 1 ) [START_REF] Auer | Finite-time Analysis of the Multi-armed Bandit Problem[END_REF], and the Thompson-Sampling (TS) algorithms [START_REF] Thompson | On the likelihood that one unknown probability exceeds another in view of the evidence of two samples[END_REF][START_REF] Agrawal | Analysis of Thompson sampling for the Multi-Armed Bandit problem[END_REF][START_REF] Kaufmann | Thompson Sampling: an Asymptotically Optimal Finite-Time Analysis[END_REF]. MAB learning has already been proposed in Cognitive Radio (CR) [START_REF] Haykin | Cognitive Radio: Brain-Empowered Wireless Communications[END_REF], and in particular, for sensing-based Dynamic Spectrum Access (DSA) in licensed bands [START_REF] Jouini | Upper Confidence Bound Based Decision Making Strategies and Dynamic Spectrum Access[END_REF]. Recently, TS and UCB 1 algorithms have been used for improving the spectrum access in (unlicensed) WiFi networks [START_REF] Toldov | A Thompson Sampling approach to channel exploration-exploitation problem in multihop cognitive radio networks[END_REF], and the UCB 1 algorithm was used in a unlicensed and frequency-and time-slotted IoT network [START_REF] Bonnefoi | Advanced metering infrastructure backhaul reliability improvement with Cogn.Radio[END_REF]. Many recent works show that MAB algorithms work well for real-world radio signal. However, even with only one dynamic user using the learning algorithm, the background traffic or the traffic of the other devices is never really stationary or i.i.d (independent and identically distributed). In recent works like [START_REF] Bonnefoi | Advanced metering infrastructure backhaul reliability improvement with Cogn.Radio[END_REF], several devices are using bandit algorithms, and the assumptions made by the stochastic bandit algorithms are not satisfied: as several agents learn simultaneously, their behavior is neither stationary nor i.i.d. As far as we know, we provide the first study to confirm robustness of the use of stochastic bandit algorithms for decision making in IoT networks with a large number of intelligent devices in the network, which makes the environment not stationary at all, violating the hypothesis required for mathematical proofs of bandit algorithms convergence and efficiency.

The aim of this article is to assess the potential gain of learning algorithms in IoT scenarios, even when the number of intelligent devices in the network increases, and the stochastic hypothesis is more and more questionable. To do that, we suppose an IoT network made of two types of devices: static devices that use only one channel (fixed in time), and dynamic devices that can choose the channel for each of their transmissions. Static devices form an interfering traffic, which could have been generated by devices using other standards as well. We first evaluate the probability of collision if dynamic devices randomly select channels (naive approach), and if a centralized controller optimally distribute them in channels (ideal approach). Then, these reference scenarios allow to evaluate the performance of UCB 1 and TS algorithms in a decentralized network, in terms of successful communication rate, as it reflects the network efficiency.

The rest of this article is organized as follows. The system model is introduced in Section 2. Reference policies are described in Section 3, and MAB algorithms are introduced in Section 4. Numerical results are presented in Section 5.

System model and notations

As illustrated in Figure 1, we suppose a slotted protocol. All devices share a synchronized time, and know in advance the finite number of available RF channels. In each time slot, devices try to send packets to the unique Base Station, which listens continuously to all channels, following an ALOHA-based communication (no sensing). Each time slot is divided in two parts: first for uplink communications in which data packets are sent by end-devices to the base station. If only one packet is sent in this part of the slot, the base station can decode it and sends an acknowledgement to the device in the second part. If two or more devices send an uplink packet in the same slot, the uplink packets collide (i.e., there is a collision), and the acknowledgement Ack is not transmitted. This way, no collision can occur on the downlink messages, easing the analysis of collisions. There are two types of end-devices in the network:

-Static end-devices have poor RF abilities, and each of them uses only one channel to communicate with the base station. Their choice is assumed to be fixed in time (stationary) and independent (i.i.d.). The traffic generated by these devices is considered as an interfering traffic for other devices. -Dynamic (or smart ) end-devices have richer RF abilities, they can use all the available channels, by quickly reconfiguring their RF transceiver on the fly. They can also store communication successes or failures they experienced in each channel, in order to change channel, possibly at every time slot.

There are N c ≥ 1 channels, D ≥ 0 dynamic end-devices, and S ≥ 0 static devices. Furthermore, in channel i ∈ 1; N c there are 0 ≤ S i ≤ S static devices (so S = Nc i=1 S i ). We focus on dense networks, in which the number of devices S + D is very large compared to N c (about 1000 to 10000, while N c is about 10 to 50). As this problem is only interesting if devices are able to communicate reasonably efficiently with the base station, we assume devices only communicate occasionally, i.e., with a low duty cycle, as it is always considered for IoT.

We suppose that all devices follow the same emission pattern, being fixed in time, and we choose to model it as a simple Bernoulli process: all devices have the same probability to send a packet in any (discrete) temporal slot, and we denote p ∈ (0, 1) this probability1 . The parameter p essentially controls the frequency of communication for each device, once the time scale is fixed (i.e., real time during two messages), and 1/p is proportional to the duty cycle.

The goal is to design a simple sequential algorithm, to be applied identically by each dynamic device, in a fully distributed setting (each device runs its own algorithm, from its observations), in order to minimize collisions and maximize the fraction of successful transmissions of all the dynamic devices.

Before explaining how this goal presents similarity with a multi-armed bandit problem, we present some natural baseline policies (i.e., algorithms).

Three reference policies

This section presents three different policies that will be used to assess the efficiency of the learning algorithms presented in the next section. The first one is naive but can be used in practice, while the two others are very efficient but require full knowledge on the system (i.e., an oracle) and are thus unpractical.

Naive policy: Random Channel Selection

We derive here the probability of having a successful transmission, for a dynamic device, in the case where all the dynamic devices make a purely random channel selection (i.e., uniform on i ∈ 1;

N c = {1, . . . , N c }).
In this case, for one dynamic device, a successful transmission happens if it is the only device to choose channel i, at that time slot. The probability of successful transmission is computed as follows, because the S i static devices in each channel i are assumed to be independent, and static and dynamic devices are assumed to not transmit at each time t with a fixed probability 1 -p :

P(success|sent) = Nc i=1 P(success|sent in channel i) No one else sent in channel i P(sent in channel i) =1/Nc,by uniform choice (1)
All dynamic devices follow the same policy in this case, so the probability of transmitting at that time in channel i for any dynamic device is p/N c , and there are D -1 other dynamic devices. As they are independent, the probability that no other dynamic device sent in i is q = P(

D-1 k=1 device k did not sent in i) = D-1 k=1 P(device k did not sent in i). And P(device k sent in i) = p × 1/N c ,

by uniform choice on channels and the Bernoulli emission hypothesis. So q

= D-1 k=1 (1 -p/N c ) = (1 -p/N c ) D-1 . Thus we can conclude, P(success|sent) = Nc i=1 (1 -p/N c ) D-1
No other dynamic device

× (1 -p) Si No static device × 1 N c = 1 N c 1 - p N c D-1 Nc i=1 (1 -p) Si . (2)
This expression (2) is constant (in time), and easy to compute numerically, but comparing the successful transmission rate of any policy against this naive policy is important, as any efficient learning algorithm should outperform it.

(Unachievable) Optimal oracle policy

We investigate in this section the optimal policy that can be achieved if the dynamic devices have a perfect knowledge of everything, and a fully centralized decision making2 is possible. We want to find the stationary repartition of devices into channels that maximizes the probability of having a successful transmission.

If the oracle draws once uniformly at random a configuration of dynamic devices, with D i devices affected to channel i is fixed (in time, i.e., stationary), then this probability is computed as before:

P(success|sent) = Nc i=1 P(success|sent in channel i) P(sent in channel i) = Nc i=1 (1 -p) Di-1 Di-1 others × (1 -p) Si No static device × D i /D Sent in channel i . (3)
Consequently, the optimal allocation vector (D 1 , . . . , D Nc ) is the solution of the following real-valued constraint optimization problem :

arg max D1,...,DN c Nc i=1 D i (1 -p) Si+Di-1 , (4a) such that Nc i=1 D i = D, (4b)
D i ≥ 0 ∀i ∈ 1; N c . ( 4c
)
Proposition 1. The Lagrange multipliers method [START_REF] Boyd | Convex Optimization[END_REF] can be used to solve the constraint real-valued maximization problem introduced in equation (4).

It gives a closed form expression for the optimal solution D * i (λ), depending on the system parameters, and the unknown Lagrange multiplier λ ∈ R.

D * i (λ) = 1 log(1 -p) W λe (1 -p) Si-1 -1 † . ( 5
)
Proof.

-In a realistic scenario, we can assume that D i ≤ -2 ln(1-p) ≈ 2 p , ∀i ∈ 1; N c . For such values for D i , the objective function f : (D 1 , . . . , D Nc ) → Nc i=1 D i (1 -p) Si+Di-1 is concave as the sum of concave functions3 .

-The Lagrange multipliers method can be applied to the optimization problem (4a), with a concave objective function f , linear equality constraints (4b) and linear inequality constraints (4c). The strong duality condition is satisfied in this case [START_REF] Boyd | Convex Optimization[END_REF], so finding the saddle points will be enough to find the maximizers.

Where (a) † = max(a, 0), and W denotes the W-Lambert function which is the reciprocal bijection of x → xe x on R + = [0, +∞) [START_REF] Corless | On the Lambert W Function[END_REF]. Moreover, condition (4b) implies that the Lagrange multiplier λ is the solution of the constraint

Nc i=1 D * i (λ) = D. (6)
Equation ( 6) can be solved numerically, with simple one-dimensional root finding algorithms. Solving the optimization problem provides the optimal real number value for D * i , which has to be rounded to find the optimal number of devices for channel i :

D i = ⌊D * i ⌋ for 1 ≤ i < N c , and
D Nc = D - Nc-1 i=1 D i .

A greedy approach of the oracle strategy

We propose a sequential approximation of the optimal policy: the third solution is a sub-optimal naive policy, simple to set up, but also unpractical as it also needs an oracle. End-devices are iteratively inserted in the channels with the lowest load (i.e., the index i minimizing S i + D i (τ ) at global time step τ ). Once the number of devices in each channel is computed, the probability of sending successfully a message is also given by equation ( 3). This is the policy that would be used by dynamic devices if they were inserted one after the other, and if they had a perfect knowledge of the channel loads.

Sequential policies based on bandit algorithms

We now present the stochastic Multi-Armed Bandit (MAB) model, and the two stochastic MAB algorithms used in our experiments [START_REF] Bubeck | Regret analysis of Stochastic and Non-Stochastic Multi-Armed Bandit Problems[END_REF]. While the stochastic MAB model has been used to describe some aspects of Cognitive Radio systems, it is in principle not suitable for our IoT model, due to the non-stationarity of the channels occupancy caused by the learning policy used by dynamic objects.

Stochastic Multi-Armed Bandits

A Multi-Armed Bandit problem is defined as follows [START_REF] Thompson | On the likelihood that one unknown probability exceeds another in view of the evidence of two samples[END_REF][START_REF] Robbins | Some aspects of the sequential design of experiments[END_REF][START_REF] Lai | Asymptotically efficient adaptive allocation rules[END_REF]. There is a fixed number N c ≥ 1 of levers, or "arms", and a player has to choose one lever at each discrete time t ≥ 1, t ∈ N, denoted as A(t) = k ∈ {1, . . . , N c }. Selecting arm k at time t yields a (random) reward, r k (t) ∈ R, and the goal of the player is to maximize the sum of his rewards, r 1...T = T t=1 r A(t) (t). A well-studied version of this problem is the so-called "stochastic" MAB, where the sequence of rewards drawn from a given arm k is assumed to be independent and identically distributed (i.i.d ) under some distribution ν k , that has a mean µ k . Several types of reward distributions have been considered, for example distributions that belong to a one-dimensional exponential family (e.g., Gaussian, Exponential, Poisson or Bernoulli distributions). We consider Bernoulli bandit models, in which r k (t) ∼ Bern(µ k ), that is, r k (t) ∈ {0, 1} and

P(r k (t) = 1) = µ k .
The problem parameters µ 1 , . . . , µ K are unknown to the player, so to maximize his cumulated rewards, he must learn the distributions of the channels, to be able to progressively focus on the best arm (i.e., the arm with largest mean). This requires to tackle the so-called exploration-exploitation dilemma: a player has to try all arms a sufficient number of times to get a robust estimate of their qualities, while not selecting the worst arms too many times.

In a Cognitive Radio application, arms model the channels, and players are the dynamic end-devices. For example in the classical OSA setting with sensing [START_REF] Jouini | Upper Confidence Bound Based Decision Making Strategies and Dynamic Spectrum Access[END_REF], a single dynamic device (a player) sequentially tries to access channels (the arms), and collects a reward of 1 if the channel is available and 0 otherwise. So rewards represent the availability of channels, and the parameter µ k represents the mean availability of channel k.

Before discussing the relevance of a multi-armed bandit model for our IoT application, we present two bandit algorithms, UCB1 and Thompson Sampling, which both strongly rely on the assumption that rewards are i.i.d..

The UCB 1 algorithm

A naive approach could be to use an empirical mean estimator of the rewards for each channel, and select the channel with highest estimated mean at each time. This greedy approach is known to fail dramatically [START_REF] Lai | Asymptotically efficient adaptive allocation rules[END_REF]. Indeed, with this policy, the selection of arms is highly dependent on the first draws, if the first transmission in a channel fails, the device will never use it again. Rather than relying on the empirical mean reward, Upper Confidence Bounds algorithms instead use a confidence interval on the unknown mean µ k of each arm, which can be viewed as adding a "bonus" exploration to the empirical mean. They follow the "optimism-in-face-of-uncertainty" principle : at each step, they play according to the best model, as the statistically best possible arm (i.e., the highest upper confidence bound) is selected.

More formally, for one device, let N k (t) = t τ =1 ½(A(τ) = k) be the number of times channel k was selected up-to time t ≥ 1. The empirical mean estimator of channel k is defined as the mean reward obtained by selecting it up to time t, µ k (t) = 1/N k (t)

t τ =1 r k (τ )½(A(τ ) = k).
For UCB 1 , the confidence term is B k (t) = α log(t)/N k (t), giving the upper confidence bound U k (t) = µ k (t) + B k (t), which is used by the device to decide the channel for communicating at time step t + 1: A(t + 1) = arg max 1≤k≤Nc U k (t). UCB 1 is an index policy.

The UCB 1 algorithm uses a parameter α > 0, originally, α set to 2 [START_REF] Auer | Finite-time Analysis of the Multi-armed Bandit Problem[END_REF], but empirically α = 1/2 is known to work better (uniformly across problems), and α > 1/2 is advised by the theory [START_REF] Bubeck | Regret analysis of Stochastic and Non-Stochastic Multi-Armed Bandit Problems[END_REF]. In our model, every dynamic device implements its own UCB 1 algorithm, independently. For one device, the time t is the number of time it accessed the network (following its Bernoulli transmission process,i.e., its duty cycle), not the total number of time slots from the beginning, as rewards are only obtained after a transmission, and IoT objects only transmit sporadicly, due to low transmission duty cycles.

Our IoT application is challenging in that there are multiple players (the dynamic devices) interacting with the same arms (the channels), without any centralized communication (they do not even know the total number of dynamic devices).

Considered alone, each dynamic device implements a learning algorithm to play a bandit game, the device is consequently a smart device. In each time slot, if it has to communicate (which happens with probability p), then it chooses a channel and it receives a reward 1 if the transmission is successful, 0 otherwise. Each device aims at maximizing the sum of the rewards collected during its communication instants, which shall indeed maximize the fraction of successful transmissions. Besides the modified time scale (rewards are no longer collected at every time step), this looks like a bandit problem. However, it cannot be modeled as a stochastic MAB, as the rewards are clearly not i.i.d : they not only depend on the (stationary, i.i.d ) behavior of the static devices, but also on the behavior of other smart devices, that is not stationary (because of learning).

Despite this, we show in the next section that running a stochastic bandit algorithm for each device based on its own rewards is surprisingly successful.

Multi-Player MAB with collision avoidance? Another idea could be to try to use a multi-player MAB model, as proposed by [START_REF] Liu | Distributed learning in Multi-Armed Bandit with multiple players[END_REF], to describe our problem.

In that case, the static and dynamic devices effect is decoupled, and arms only model the availability of the channels in the absence of dynamic devices : they are i.i.d. with mean µ i = 1 -pS i . Moreover, dynamic devices are assumed to be able to sense a channel before sending [START_REF] Liu | Distributed learning in Multi-Armed Bandit with multiple players[END_REF], and so communicate only if no static device is detected on the channel. The smart devices try to learn the arms with highest means, while coordinating to choose different arms, i.e., avoid collisions in their choice, in a decentralized manner. However, in this model it is assumed that the multiple agents can know that they experienced a collision with another agent, which is non-realistic for our problem at stake, as our model of smart device cannot do sensing nor differentiate collisions between smart and non-smart devices.

Adversarial bandit algorithms? Instead of using MAB algorithms assuming a stochastic hypothesis on the system, we could try to use MAB algorithms designed to tackle a more general problem, that makes no hypothesis on the interfering traffic. The adversarial MAB algorithms is a broader family, and a well-known and efficient example is the Exp3 algorithm [START_REF] Bubeck | Regret analysis of Stochastic and Non-Stochastic Multi-Armed Bandit Problems[END_REF]. Empirically, the Exp3 algorithm turned out to perform worse than both UCB 1 and TS in the same experiments. Contrarily to the two stochastic algorithms, the use of Exp3 is correctly justified, even in the non-stationary and non-i.i.d, as its performance guarantee are true in any setting. But it is not so surprising that it performs worse, as the theoretical performance guarantees of adversarial MAB algorithms are an order of magnitude worse than the one for stochastic ones. More is left on this aspect for our future work.

Experiments and numerical results

We suppose a network with S + D = 2000 end-devices, and one IoT base station. Each device sends packets following a Bernoulli process, of probability p = 10 -3 (e.g., this is realistic: one packet sent about every 20 minutes, for time slots of 1s). The RF band is divided in N c = 10 channels4 . Each static device only uses one channel, and their uneven repartition 4 in the 10 channels is: (S 1 , • • • , S Nc ) = S × (0.3, 0.2, 0.1, 0.1, 0.05, 0.05, 0.02, 0.08, 0.01, 0.09), to keep the same proportions when S decreases. The dynamic devices have access to all the channels, and use learning algorithms. We simulate the network during 10 6 discrete time slots, during which each device transmits on average 1000 packets (i.e., the learning time is about 1000 steps, for each algorithm).

Figure 2 presents the evolution of the successful transmission rate, as a function of time. The two MAB algorithms, UCB 1 Thompson Sampling (TS), are compared against the naive random policy from below, and the two oracle policies (optimal and greedy) from above. The results are displayed when 10, 30, 50 and 100% of the traffic is generated by dynamic devices. We can see in Figure 2 that the TS algorithm (in red) outperforms the UCB 1 algorithm (in blue), when the number of end-devices is below 50%. When the number of end-devices is higher, both algorithms have almost the same performance, and perform well after very few transmissions (quick convergence). Moreover, we can see in Figures 2a,2b, and 2c that both have better success rate than the random policy and the probability of successful transmission is between the oracle optimal and oracle suboptimal policies. For instance, for 10% of dynamic devices, after about 1000 transmissions, using UCB 1 over the naive uniform policy improved the successful transmission rate from 83% to 88%, and using Thompson Sampling improved it to 89%. Increasing the number of end-devices decreases the gap between the optimal and random policies: the more dynamic devices, the less useful are learning algorithms, and basically for networks with only dynamic devices, the random policy is as efficient as the optimal one, as seen in Figures 2d and3.

To better assess the evolution of the optimal policy compared to the random one, we have displayed on Figure 3 the evolution of the gain, in term of successful transmissions rate, provided by the optimal oracle and the two learning policies, after 10 6 time slots, i.e., about 1000 transmissions for each object. We can see that when the proportion of end-devices is low (e.g., 1% of devices are dynamic), the optimal policy provides an improvement of 16% compared to random channel selection. The TS algorithm always provides near-optimal performance, but the UCB 1 algorithm has a lowest rate of convergence and performs consequently worse after 1000 transmissions, for instance it only provides a gain of 12% for the same proportion of dynamic devices (1%). Figure 3 also shows that learning keeps near-optimal performance even when the proportion of devices becomes large. Note that when this proportion increases, the assumptions of a stochastic MAB model are clearly violated, and there is no justification for the efficiency of TS and UCB 1 algorithms. Hence, it is surprising to have near optimal performance with stochastic MAB algorithms applied to partly dynamic and fully dynamic scenarios.

Conclusion

In this article, we proposed an evaluation of the performance of MAB learning algorithms in IoT networks, with a focus on the convergence of algorithms, in terms of successful transmission rates, when the proportion of intelligent dynamic devices changes. Concretely, increasing this probability allows to insert more objects in the same network, while maintaining a good Quality of Service. We show that UCB 1 and TS have near-optimal performance, even when their underlying i.i.d. assumption is violated by the many "intelligent" end-devices. This is both a surprising and a very encouraging result, showing that application of bandit algorithms tailored for a stochastic model is still useful in broader settings. The fully decentralized application of classic stochastic MAB algorithms are almost as efficient as the best possible centralized policy in this setting, after a short learning period, even though the dynamic devices can not communicate with each others, and do not know the system parameters. We will investigate this behavior in order to understand it better theoretically. We will also experiment more with adversarial algorithms, to confirm that they work less efficiently than stochastic bandit algorithms in our non-stochastic setting.

Moreover, for sake of simplicity we supposed that all devices use the same standard. Our future work will consider more realistic interference scenarios and IoT networks, with, e.g., non-slotted time, more than one base station etc.

Fig. 1 .

 1 Fig. 1. The considered time-frequency slotted protocol. Each frame is composed by a fix-duration uplink slot in which the end-devices transmit their packets. If a packet is well received, the base station replies by transmitting an Ack, after the ack delay.

In the experiments below, p is about 10 -3 , because in a crowded network p should be smaller than Nc/(S + D) for all devices to communicate successfully (in average).

This optimal policy needs an oracle seeing the entire system, and affecting all the dynamic devices, once and for all, in order to avoid any signaling overhead.

It worth noting that f is neither concave nor quasi-concave on [0, ∞) Nc[START_REF] Luenberger | Quasi-convex programming[END_REF][START_REF] Yaari | A note on separability and quasiconcavity[END_REF].

We tried similar experiments with other values for Nc and this repartition vector, and results were similar for non-homogeneous repartitions. Clearly, the problem is less interesting for homogeneous repartition, as all channels appear the same for dynamic devices, and so even with D small in comparison to S, the system behaves like in Fig.2d, where the performance of the five approaches are very close.

Acknowledgements This work is supported by the French National Research Agency (ANR), under the projects SOGREEN (grant coded: N ANR-14-CE28-0025-02) and BADASS (N ANR-16-CE40-0002), by Région Bretagne, France, by the French Ministry of Higher Education and Research (MENESR) and ENS Paris-Saclay.

Note: the simulation code used for the experiments in Section 5 is for MATLAB or GNU Octave, and is open-sourced under the MIT License, at: https://Bitbucket.org/scee_ietr/rl_slotted_iot_networks.


C Aldebert
email: clement.aldebert@univ-amu.fr

D Nerini
email: david.nerini@univ-amu.fr

M Gauduchon
email: mathias.gauduchon@univ-amu.fr

J C Poggiale

Does structural sensitivity alter complexity-stability relationships?

Keywords: food webs, dynamical systems, structural sensitivity, complexity-stability relationships, functional response Highlights

We study food webs of varying trophic complexity (number of species, connectance)

• Food web dynamics is modeled using either Holling's or Ivlev's functional response

Community persistence is mainly driven by trophic complexity • Occurrence of fluctuating dynamics is mainly driven by the functional response • Conclusions are robust to some changes in parameter values and model assumptions 1

Introduction

Predictions made by mathematical models can be sensitive to model formulation [START_REF] Anderson | Influence of grazing formulations on the emergent properties of a complex ecosystem model in a global ocean general circulation model[END_REF]Fulton et al., 2003b;Fussmann & Blasius, 2005, among others). However, this sensitivity hass rarely been tested in theoretical and operational ecosystem models [START_REF] Arhonditsis | Evaluation of the current state of mechanistic aquatic biogeochemical modeling[END_REF]Fulton et al., 2003a).

In ecological models with multiple interacting populations, phenomena observed at the community scale are usually represented by simplifying smaller scale processes. For instance, collective and individual behaviours as well as physiological processes involved in predation are collapsed into one function, the functional response [START_REF] Gentleman | Functional responses for zooplancton feeding on multiple resources: a review of assumptions and biological dynamics[END_REF][START_REF] Jeschke | Predator functional response: discriminating between handling and digesting prey[END_REF]. Numerous mathematical formulations of a given biological phenomenon are relevant in the sense that: (i) their properties and assumptions about underlying processes are consistent with the knowledge of the system to model, (ii) they equivalently fit empirical data [START_REF] Cordoleani | Structural sensitivity of biological models revisited[END_REF][START_REF] Mullin | Ingestion by planktonic grazers as a function of food concentration[END_REF]. Moreover, some of these functions may have the same mathematical properties (pointwise properties, monotonicity, convexity, ...). However, the choice of a particular function among relevant ones can affect the dynamics predicted by the same model.

Differences occur in predicted steady-state values, equilibrium vs. oscillating dynamics and in the system response to external disturbances [START_REF] Aldebert | Structural sensitivity and resilience in a predator-prey model with density-dependent mortality[END_REF]. Uncertainty due to this choice of a function is coined as "structural sensitivity" (sensu [START_REF] Cordoleani | Structural sensitivity of biological models revisited[END_REF].

Structural sensitivity has been theoretically studied in simple models with a few state variables, mainly predator-prey and food chain models [START_REF] Adamson | When can we trust our model predictions? Unearthing structural sensitivity in biological systems[END_REF][START_REF] Adamson | Defining and detecting structural sensitivity in biological models: developing a new framework[END_REF][START_REF] Fussmann | Community response to enrichment is highly sensitive to model structure[END_REF][START_REF] Gross | Enrichment and foodchain stability: the impact of different forms of predator-prey interaction[END_REF][START_REF] Myerscough | Stability, persistence and structural stability in a classical predator-prey model[END_REF]. The aim of this study is to extend these results to more complex models such as food webs. Previous results on predator-prey models may suggest that food web models are sensitive to the choice of type-II functional response. We propose to compare food web dynamics under changes in both functional response formulation and trophic complexity (number of trophic species and trophic links).

Trophic complexity is known to affect food web dynamics and stability. Complexity-stability relationships have been conceptually studied by [START_REF] Macarthur | Fluctuations of animal populations and a measure of commmunity stability[END_REF] and then more formally by [START_REF] May | Will a large complex system be stable[END_REF][START_REF] May | Qualitative stability in model ecosystems[END_REF]).

May's work has led to a long-standing debate which is still open after decades of field and theoretical researches [START_REF] Loreau | From populations to ecosystems: theoretical foundations for a new ecological synthesis[END_REF][START_REF] May | Unanswered questions in ecology[END_REF][START_REF] Mccann | The diversity-stability debate[END_REF]. Food webs exhibit a huge number of different structures. A relevant analysis of their common properties requires to reproduce their diversity. Numerous food webs with empirically consistent structural properties and a desired trophic complexity can be built by simple random models [START_REF] Cattin | Phylogenetic constraints and adaptation explain food-web structure[END_REF]Williams & Martinez, 2000, among others). These models have been used to statistically investigate complexity-stability relationships in food web models based on different ecological phenomena [START_REF] Brose | Allometric scaling enhances stability in complex food webs[END_REF][START_REF] Heckmann | Interactive effects of body-size structure and adaptative foraging on food-web stability[END_REF][START_REF] Kartascheff | Positive complexity-stability relations in food web models without foraging adaptation[END_REF][START_REF] Kartascheff | Why allometric scaling enhances stability in food web models[END_REF][START_REF] Plitzko | Complexity-stability relations in generalized food-web models with realistic parameters[END_REF][START_REF] Stouffer | Understanding food-web persistence from local to global scales[END_REF], 2011;[START_REF] Uchida | Relation between complexity and stability in food webs with adaptive behavior[END_REF][START_REF] Williams | Effects of network and dynamical model structure on species persistence in large model food webs[END_REF][START_REF] Williams | Stabilization of chaotic and non-permanent food-web dynamics[END_REF].

The questions we address in this paper are: (i) are dynamics predicted by a food web model more impacted by structural sensitivity or by trophic complexity? (ii) Does structural sensitivity alter complexitystability relationships?

Next section presents the studied food web model. It is an extension of a predator-prey model in which structural sensitivity has already been explored [START_REF] Aldebert | Structural sensitivity and resilience in a predator-prey model with density-dependent mortality[END_REF]. Structural sensitivity in this model is compared to the impact of trophic complexity in section 3.1. Observed results are then explained from the knowledge of predator-prey models (section 3.2) and their robustness to changes in the method used to fit functional responses is tested (section 3.3). Then, the relative importance of trophic complexity, functional response formulation and parameter values is estimated (section 3.4).

Next, robustness to changes in model assumptions is assessed, and complexity-stability relationships are compared to empirical findings (section 3.5). Paper ends with a more general discussion about structural sensitivity and modelling of biological systems (section 3.6).

Models

Food web structure

Food webs are composed by S species (sensu trophic species) and one resource. Species are linked by L trophic interactions, so that food web connectance is C = L/S 2 (directed connectance, [START_REF] Martinez | Artifacts or attributes? Effects of resolution on the Little Rock Lake food web[END_REF].

The niche model [START_REF] Williams | Simple rules yield complex food webs[END_REF] is used to randomly build numerous food webs with the desired number of species and connectance. The niche model generates quickly numerous food webs with patterns that are consistent with empirical data [START_REF] Allesina | A general model for food web structure[END_REF][START_REF] Cattin | Phylogenetic constraints and adaptation explain food-web structure[END_REF][START_REF] Williams | Simple rules yield complex food webs[END_REF]. It is based on the principle of ecological niche [START_REF] Hutchinson | Concluding remarks[END_REF]. A species i is characterized by a niche value n i uniformly drawn in the interval [0, 1], the niche axis.

The niche model is described in section 1 of Supporting Online Material (SOM). Food webs are made of distinct species, that are either a primary producer or a predator. Attribution of trophic links allows for cannibalism and trophic loops. We added a rejection step after food webs construction to avoid unrealistic patterns. We only studied food webs with a realized connectance that deviated at most by 0.01 of the expected one, that are connected (no disconnected parts), and in which all predators feed (as a prey or through a food chain) upon at least one primary producer.

Food web dynamics

Food web dynamics is modeled using a dynamical system of S differential equations. It is a bio-energetic model extended for a multi-species system [START_REF] Plitzko | Complexity-stability relations in generalized food-web models with realistic parameters[END_REF]Yodzis & Innes, 1992, among others). This deterministic model is continuous in time with unstructured populations. Each species i is described by its biomass B i , with dynamics given by the ordinary differential equation:

dB i dt = λq φ i B i + λ j∈Ri G φ i,j B i - j∈Ci G φ j,i B j -α i B i -β i B 2 i i = 1, ..., S. (1)
Right terms of model ( 1) handle respectively a gain in biomass by primary production, sum of gains by predation, sum of losses by predation, linear mortality and respiration, density-dependent mortality (intra-specific competition, diseases). Species i possesses a set of prey (predator) species denoted as R i (C i ). By definition, primary producers have R i = ∅ and top-predators have C i = ∅. The parameter λ is the assimilation efficiency. For the sake of simplicity, λ is assumed to be the same for all species.

Parameter α i is the linear mortality rate and parameter β i is the per-capita intra-specific competition rate of species i. The letter φ indicates the specific formulation used for the Holling-type II functional response G φ i,j . For simplicity, all species are assumed to have the same formulation. This one is either Holling's disc equation (1959;[START_REF] Holling | The functional response of predators to prey density and its role in mimicry and population regulation[END_REF] denoted as G H i,j or Ivlev's functional response (1955) denoted as G I i,j

(later called Holling's FR and Ivlev's FR):

G H i,j = a H i f i,j B j 1 + h H i a H i T i , G I i,j = 1 h I i 1 -exp -h I i a I i T i f i,j B j T i with T i = j∈Ri f i,j B j .
Both functional responses are extended for a predator with multiple prey species by assuming that it does not switch between preys [START_REF] Gentleman | Functional responses for zooplancton feeding on multiple resources: a review of assumptions and biological dynamics[END_REF]. For Holling's FR, parameters a H i and h H i are respectively the attack rate and the handling time of the predator. For Ivlev's FR, parameter 1/h I i is the maximal digestion rate and a I i h I i is the satiation coefficient of the predator. The total amount of prey available for species i is the weighted sum of its prey species biomass T i . The weighting parameter f i,j is constant and it can be considered as the foraging effort or the feeding preference of predator i for its prey species j (obviously, f i,j = 0 and G φ i,j = 0 if j / ∈ R i ). This means that the total functional response of a predator G φ,tot i

(T i ) = j∈Ri G φ i,j is a function of T i (figure 1) and that G φ i,j = G φ,tot i f i,j B j /T i . Both
functional responses also fulfills properties:

G φ,tot i ∈ C 2 , G φ,tot i (0) = 0, G φ,tot i (T i ) ≥ 0, G φ,tot i (T i ) > 0, G φ,tot i (T i ) < 0, lim Ti→+∞ G φ,tot i (T i ) < +∞,
with C 2 being the class of twice continuously differentiable functions. Other properties means that G φ,tot is null in absence of prey, increases with prey biomass, is concave and saturates at high prey biomass.

Functional response's parameters have the same mathematical meaning in both formulations:

G φ,tot i (0) = a φ i , lim Ti→+∞ G φ,tot i (T i ) = 1 h φ i .
Thus, a φ pred gives the slope of the functional response at the origin, and 1/h φ pred gives the asymptotic value of the functional response when it saturates at high prey biomass.

The term of primary productivity q φ i has the same equation as the functional response G φ i,j with a constant pool of resources B res :

q H i =          a H i B res 1 + h H i a H i B res if i ∈ P I 0 otherwise , q I i =          1 h I i 1 -exp -h I i a I i B res if i ∈ P I 0 otherwise.
This term vanishes for predators, with P I being the set of primary producers. Setting a constant pool of resource implies two assumptions. First, the time scale of the resource dynamics is such that it is always at quasi-equilibrium in comparison to the time scale of population dynamics. Second, the equilibrium value is independent of the food web dynamics, it only depends on the environment, which is supposed to be constant.

Let M i = 10 xni be the body mass of species i, with a scale parameter x. We assume that some parameter values scale allometrically across species [START_REF] Brown | Toward a metabolic theory of ecology[END_REF]: total prey biomass available (T i )

Figure 1. Functional responses used in the model. a: total functional response of species i (total amount of preys eaten by predator and by time unit) calculated as a function of its total prey biomass available, using Holling's disc equation (solid) or Ivlev's functional response (dashed). The latter is parameterized in order to minimize the weighted Euclidean distance between the two formulations. For the sake of visibility, only a part of the fitting range is shown. b: weighting function used to calculate the distance between formulations. See section 2.2 for formulations and section 2.3 for parameterization.

a φ i = a φ M -1/4 i , h φ i = h φ M 1/4 i , α i = αM -1/4 i , β i = βM -1/4 i .
As species abundance is quantified in terms of biomass, the scaling exponent -1/4 corresponds to an increase of individual metabolic rate with body mass power 3/4, divided by individual body mass. The 3/4 exponent is controversial. Empirical data and theoretical studies indicate values between 2/3 and 1 [START_REF] Kooijman | Dynamic energy budget theory for metabolic organisation[END_REF]. Nevertheless, the exact value of this exponent has a limited impact on species extinctions in similar food web models [START_REF] Kartascheff | Why allometric scaling enhances stability in food web models[END_REF].

Note that if two species have the same trophic interactions and parameter values, model ( 1) is sensitive to the aggregation of these identical species as no direct inter-specific competition is considered. It has no impact on our numerical results as identical (or infinitely close) species are infinitely rare in food webs built by the niche model. This problem in model consistency can be solved by considering an inter-specific competition term, with a competition strength that decreases with a measure of distance between species.

S A,ω a I , h I = d 2 A,ω G H,tot i , G I,tot i = A ω(T i ) G H,tot i (T i ) -G I,tot i (T i ) 2 dT i , (2)
with a weighting function ω(T i ) on a fitting range A of T i values. Total prey biomass ranges from 0 to constant resource biomass for primary producers, so we set A = [0, B res ]. For a fixed amount of total prey biomass j∈Ri B j , the resulting T i depends on |R i |. To balance this effect a priori, the distance is weighted by (figure 1):

ω P (T i ) = 1 S max T max N n=1
n with N = min ( T max /T i , S max ) and

T max = max(A), (3)
which is the frequency distribution of T i values for j∈Ri B j ∈ A and for a predator's number of prey species which ranges between one and S max = 60. Minimizing the cost function S [0,Bres],ω P using the simplex method by [START_REF] Nelder | A simplex method for function minimization[END_REF] give values a I = 3.17 and h I = 0.36. The sensitivity of our results to choices in A and ω(T i ) is discussed in section 3.3.

Let us recall that, when an equilibrium is reached, the right-hand side of equation ( 1) vanishes for all species, with B * 1 , B * 2 , ..., B * S being the species biomass at an equilibrium (and T * i the total prey biomass available of species i at equilibrium). The jacobian matrix J φ of system (1) is presented in appendix S1.

The algorithm presented in section 2 of SOM (with a discussion about its limits) is used to automatically determine the kind of asymptotic dynamics reached by a food web. The proposed algorithm classifies asymptotic dynamics in three types: dynamics in which at least one species goes extinct, equilibrium in which all species coexist, fluctuating dynamics (e.g. limit cycle, torus, strange attractor) in which all species coexist. The proportion of food webs in which all species coexist (at equilibrium or with fluctuating dynamics) is a measure of food web persistence. The proportion of persistent food webs (i.e. food webs in which all species coexist) which exhibit a fluctuating dynamics is a measure of food web variability. Food webs construction and simulation were achieved by C++ programs using GNU Scientific Library for C/C++ [START_REF] Galassi | GNU Scientific Library Reference Manual[END_REF]. Post-simulation analysis and figures were performed using R language (R Core Team, 2013).

Results and Discussion

Structural sensitivity vs. trophic complexity

Both functional responses predict similar patterns of food web persistence as a function of trophic complexity (figure 2a,c). Within this range, Ivlev's FR predicts less extinctions than Holling's FR does (from -11 % to +1 % food webs with at least one species extinction). With both functional responses, the proportion of non-persistent food webs increases with trophic complexity. This increase is higher with the number of species (×1.4 with Holling's FR and ×1.5 with Ivlev's FR) than with connectance (×1.2 and ×1.3). Functional responses are parameterized to predict quantitatively close predation fluxes, which explains their similar predictions of food web persistence.

Conversely, food web variability is mainly driven by functional response formulation within the tested range of trophic complexity (figure 2b,d). Ivlev's FR predicts 2.3 times more equilibrium dynamics (from 98 % to 99 % of persistent food webs reach an equilibrium) than Holling's FR does (from 25 % to 80 %). The proportion of equilibrium dynamics increases more with connectance (on average ×3.2 with Holling's FR and +1 % with Ivlev's FR) than with the number of species (×1.1 with Holling's FR, no trend with Ivlev's FR). Even if the effect of trophic complexity can change with other measures of food web persistence and variability, the latter is still higher with Ivlev's FR (section 5 in SOM).

Understanding structural sensitivity in food webs

Model (1) applied to a one-predator/one-prey system owns an equibrium which is stable with Ivlev's FR and unstable with Holling's FR in 26 % to 49 % of the parameter space explored in [START_REF] Aldebert | Structural sensitivity and resilience in a predator-prey model with density-dependent mortality[END_REF].

With Holling's FR, dynamics converge on a stable limit cycle. Equilibrium stability is different because the slope of Holling's FR at equilibrium is lower. The slope of the functional response at equilibrium in food webs can be described by its elasticity:

γ i := g φ i (1) ∈ [0, 1] with t i := T i T * i and g φ i (t i ) := G φ,tot i (t i T * i ) G φ,tot i (T * i )
,

where the normalized functional response equals 1 at equilibrium (g φ i (1) = 1). A high elasticity γ i stabilizes equilibria in predator-prey [START_REF] Aldebert | Structural sensitivity and resilience in a predator-prey model with density-dependent mortality[END_REF][START_REF] Yeakel | Generalized modeling of ecological population dynamics[END_REF] and food web models (section 6 of SOM).

The elasticity depends on prey biomass at equilibrium, which has a similar distribution between functional responses (figure 3). However, Ivlev's FR has a higher elasticity than Holling's FR for low prey biomass values (figure 4b). As a consequence, the realized γ i distribution in food web simulations is 1.6 times higher with Ivlev's FR (γ = 0.60) than with Holling's FR (γ = 0.37, figure 4a). With the latter, one can expect that stable equilibria are closer to a bifurcation threshold (with respect to γ i ).

Thus, more food webs are likely to be on the unstable side of this bifurcation with Holling's FR. This can explain why this function leads more frequently to fluctuating dynamics. This reasoning extends the previous mechanism found in predator-prey models [START_REF] Fussmann | Community response to enrichment is highly sensitive to model structure[END_REF] to the scale of complex

Holling's disc equation (b) fraction of persistent food webs with a stable equilibrium q q q q q q q q q q q q 20 30 40 Ivlev's (dashed) functional response. Computed from 160 000 food webs of varying complexity (4 levels of connectance and 4 numbers of species, 10 000 food webs by pair of values). Density estimates are realized using non-parametric kernel methods [START_REF] Simonoff | Smoothing methods in statistics[END_REF].

Structural sensitivity and the method used to fit the functional responses

As the slope of the functional response drives food web variability, one can think about different ways to parameterize Ivlev's FR. The relevance of these different approaches from a biological point of view is discussed in section 3.6. Here, we focus on robustness to changes in parameter values estimated for Ivlev's FR. Different parameter sets are estimated by minimizing the weighted distance ( 2): (i) over different ranges of T i values, and (ii) using either the weighting function ω P (3) (figure 1) or a uniform weighting function ω U (T i ). One last parameter set is estimated using the empirical distribution of total prey biomass observed in persistent food webs at equilibrium with Holling's FR (figure 3b).

Parameter sets that give a better fit to Holling's FR at low prey biomass (closer a φ value) predict a closer value of food web variability between functional responses (table 2 and figure 5). However, they also predict more distant primary productivity (≈ 1/h φ value) and predicted value of food web persistence.

Furthermore, equilibrium dynamics are still more frequent with Ivlev's FR within the tested range of Computed for predator species from 160 000 food webs of varying complexity (4 levels of connectance and 4 numbers of species, 10 000 food webs by pair of values). Density estimates are realized using non-parametric kernel methods [START_REF] Simonoff | Smoothing methods in statistics[END_REF].

trophic complexity (section 3 in SOM). The parameter set that leads to the closer food web variability (ω P within [0, 3.5]) does not lead to the best fit of functional responses within the range of prey biomass observed (column S ω U ,[0,3.5] in table 2), as it is optimized for even smaller biomass values. In other words, this parameter set optimizes only the slope at the origin (a φ ), as we know that the functional response slope drives food web variability. Only one parameter set is usually considered in previous studies on structural sensitivity in predator-prey models [START_REF] Aldebert | Structural sensitivity and resilience in a predator-prey model with density-dependent mortality[END_REF][START_REF] Cordoleani | Structural sensitivity of biological models revisited[END_REF][START_REF] Fussmann | Community response to enrichment is highly sensitive to model structure[END_REF][START_REF] Myerscough | Stability, persistence and structural stability in a classical predator-prey model[END_REF]. We checked here with different parameter sets that functional response formulation influences food web variability, even if this influence is lower with some fitting methods.

Structural sensitivity vs. parameter sensitivity

Even with the simplifying assumption that parameter values scale allometrically, a full parameter sensitivity analysis requires a high computational effort. Indeed, to be reasonably exhaustive it is necessary to study thousands of food webs with different numbers of species and connectance levels for each pa-Table 2. Impact of the method used to fit the functional responses. The columns indicate respectively: the weighting function used, the range of T i values within which the weighted distance between functions is minimized, the pair of parameters values obtained (a φ and h φ ) for Ivlev's functional response, the corresponding distance (computed with the same weighting function for comparison) between functions at low ([0, 3.5]) and high ([3.5, 500]) total prey biomass available, the predicted dynamics (proportion of food webs with extinction(s), proportion of persistent food webs reaching an equilibrium) in 30-species food webs with a connectance of 0.15. The number of food webs studied for each pair of parameters is enough to obtain 1 000 food webs with a stable positive equilibrium. The bold line corresponds to the parameters values used in the rest of the study. In addition, one parameter set for Ivlev's functional response has been estimated using the empirical distribution of T i observed in persistent food webs at equilibrium with Holling's disc equation (figure 3b). Dynamics predicted by Holling's disc equation are indicated for comparison.

weighting function [3.5,500] (discretized by steps of 10 %) and re-estimate the parameters for Ivlev's FR by fitting the two functions.

T i ∈ A a φ h φ S ω U ,[0,3.5] S ω U ,
The relative variation considered is comparable to uncertainties in parameter estimation from empirical data [START_REF] Cordoleani | Structural sensitivity of biological models revisited[END_REF]. Food webs with the same numbers of species and connectance levels as in figure 2 are studied for each parameter set and each functional response for a total amount of ≈ 6.10 7 food webs studied.

Within the tested range of trophic complexity, food web persistence is strongly correlated (0.76) to the number of species, whereas food web variability is strongly correlated (0.82) to functional response formulation (3). Both persistence and variability are also correlated (0.37 and 0.36) to connectance.

In 2). b and c: food web dynamics (30 species, connectance of 0.15) predicted for these parameter sets. Gray labels and arrows summarize the main trend of results: a better fit at low prey biomass (a) leads to a closer predicted food web variability (c), but also to more distant functions at other prey biomass values (a) and more distant predicted food web persistence (b). correspond to a quantitatively similar variation measured between the functional responses [START_REF] Adamson | When can we trust our model predictions? Unearthing structural sensitivity in biological systems[END_REF][START_REF] Adamson | Defining and detecting structural sensitivity in biological models: developing a new framework[END_REF][START_REF] Cordoleani | Structural sensitivity of biological models revisited[END_REF][START_REF] Wood | Super-sensitivity to structure in biological models[END_REF]. Indeed, for a quantitatively similar variation, a change of formulation may affect the function shape and lead to a higher change in functional response slope near equilibrium, i.e. in equilibrium stability.

Model assumptions, complexity-stability relationships and biological systems

Different combinations of model assumptions are made by recent studies on food webs, often for simplifying reasons [START_REF] Brose | Allometric scaling enhances stability in complex food webs[END_REF][START_REF] Gross | Generalized models reveal stabilizing factors in food webs[END_REF][START_REF] Heckmann | Interactive effects of body-size structure and adaptative foraging on food-web stability[END_REF][START_REF] Kartascheff | Why allometric scaling enhances stability in food web models[END_REF][START_REF] Plitzko | Complexity-stability relations in generalized food-web models with realistic parameters[END_REF][START_REF] Stouffer | Understanding food-web persistence from local to global scales[END_REF], 2011;[START_REF] Williams | Effects of network and dynamical model structure on species persistence in large model food webs[END_REF]. The impact of these assumptions is detailled in section 4 of SOM and summarized in table 4. In all cases, equilibrium dynamics are more frequent with Ivlev's FR than with Holling's FR. In addition, both functional responses show the same qualitative effect of each assumption.

Table 3. Correlations between food web dynamics, trophic complexity, the functional response used and its parameterization. We simulated food webs of varying complexity (4 numbers of species times 4 connectance levels) with different functional response formulation (Holling's or Ivlev's FRs) and parameterization (5 values of a φ times 5 values of h φ ), leading to n = 800 combinations. For each one, food web dynamics is summarized by the proportion of food webs with extinction(s) and the proportion of persistent food webs reaching a stable equilibrium. A positive correlation with the functional response corresponds to an increase with Ivlev's FR. NS means that the corresponding p-value is higher than 0.01 (Student's test on the regression coefficient). For each combination, the number of food webs studied is enough to obtain 1 000 food webs with a stable positive equilibrium. The total number of food webs simulated is ≈ 6.10 7 . number of species connectance a φ h φ functional response extinction(s) 0.76 0.37 0.12 0.07 -0.17 equilibrium NS 0.36 -0.22 0.07 0.82 Table 4. Model assumptions and complexity-stability relationships. Three changes in model assumptions are tested: studying food webs with a fixed number of five primary producer species (PP), setting the body mass of primary producers to 1 and deleting cannibalistic links. For each changes, results (proportion of food webs with extinction(s) and proportion of persistent food webs that reach an equilibrium) are summarized for Holling's and Ivlev's functional responses. Trophic complexity (number of species S and connectance C) is varied within the same range as in figure 2. The obtained range of results and the mean impact of trophic complexity are shown. The "-" indicates that the impact of the number of species is not discussed because the number of persistent food webs studied has been decreased for computational reasons. These changes in model assumptions show that food web persistence decreases with connectance due to a decrease of the proportion of primary producer species in food webs made by the niche model (SOM, figure 2). However, persistence is increased by the connectance per se, i.e. the number of pathways where energy can flow. Furthermore, connectance increases the occurrence of equilibrium dynamics and improves species survival by keeping their dynamics far from the extinction threshold. More equilibrium dynamics occur at high connectance because it increases the number of weak trophic links. Weak links stabilize equilibria by dampening predator-prey oscillations [START_REF] Mccann | Weak trophic interaction and the balance of nature[END_REF].

The impact of the number of species on system dynamics can be related to general questions on the effect of biodiversity on ecosystem functioning [START_REF] May | Unanswered questions in ecology[END_REF][START_REF] Mccann | The diversity-stability debate[END_REF]. [START_REF] Jiang | Different effects of species diversity on temporal stability in single-trophic and multitrophic communities[END_REF] have made a meta-analysis of empirical studies on the effect of biodiversity. Both observational and experimental studies on multi-trophic communities indicate that the temporal variability of state variables (biomasses, processes, etc.) decreases at the population level (15 studies), and even more at the community level (25 studies). These trends are predicted by the model, which predicts more equilibrium dynamics (figure 2 and table 4) and a lower temporal variability of population biomass (section 5 in SOM) in larger persistent food webs. In addition, larger food webs have a lower variability of their predicted total biomass (data not shown). Note that in the model we used, there is no feedback of population dynamics on the trophic structure, which is set a priori. Thus, the link between diversity and ecosystem persistence can be better explored using evolutionary models that build persistent networks of interaction based on population dynamics [START_REF] Drossel | The influence of predator-prey population dynamics on the long-term evolution of food web structure[END_REF][START_REF] Loeuille | Evolutionary emergence of size-structure food webs[END_REF][START_REF] Rossberg | Food webs: experts consuming families of experts[END_REF].

Structural sensitivity and modelling biological systems

We have shown that quantitatively close mathematical functions can lead to qualitatively different food web dynamics. However, the amount of differences in dynamics depends on the way the functions were parameterized. In general, parameterization is realized by either fitting the functions to data on the process itself [START_REF] Anderson | Influence of grazing formulations on the emergent properties of a complex ecosystem model in a global ocean general circulation model[END_REF][START_REF] Cordoleani | Structural sensitivity of biological models revisited[END_REF][START_REF] Fussmann | Crossing the Hopf bifurcation in a live predator-prey system[END_REF][START_REF] Poggiale | How far details are important in ecosystem modelling: the case of multi-limiting nutrients in phytoplankton-zooplankton interactions[END_REF],

or fitting the full model predictions (other parameters can be optimized at the same time) to data on the temporal evolution of the system [START_REF] Canale | Experimental and mathematical modeling studies of protozoan predation on bacteria[END_REF][START_REF] Kooi | The transient behaviour of food chains in chemostat[END_REF].

Here, we used the first method in different ways. The most classical ones, namely a fit with a uniform weight over the range of biomass observed or a fit using the empirical biomass distribution, lead to significant differences in model dynamics between functional responses. These differences become lower only if the fit is done so that functions are close in term of property that drives model dynamics (here the functional response slope at the origin). So, model dynamics have to be known a priori to be able to use this approach, which is in fact more similar to the second method of fit.

The second method focuses on the dynamics predicted by the model. Different functions can be parameterized in order to predict close dynamics. However, these dynamics are specific to the set of environmental conditions (i.e. other parameter values) used. If these conditions are changed (e.g. resource availability, mortality due to external factors), the different functions are likely to predict different dynamics as they are outside the range of optimization. The underlying idea is similar to the method of generalized modelling [START_REF] Gross | Generalized models as a universal approach to the analysis of nonlinear dynamical systems[END_REF] and the approach proposed by [START_REF] Adamson | When can we trust our model predictions? Unearthing structural sensitivity in biological systems[END_REF].

These approaches do not specify model formulation and describe only the local properties of the model (like the slope of a function) near an equilibrium. So, results on equilibrium stability are independent of a specific formulation. However, these approaches provide no information on the system dynamics far from this equilibrium, like the existence of alternative stable states. In conclusion, if different functions are parameterized to predict similar dynamics in a given situation, this is likely to work only in the vicinity of this given situation (e.g. a given equilibrium). Choosing to do so or to use the first method depends on the purpose of the study.

Stable equilibria and limit cycles can coexist in the food web model applied to a one-predator/oneprey system, due to density-dependent mortality [START_REF] Aldebert | Structural sensitivity and resilience in a predator-prey model with density-dependent mortality[END_REF]. Density-dependent mortality allows the coexistence of different asymptotic states that depend on functional response formulation.

As a consequence, functional response formulation drives model predictions in situations where external disturbances and recovery policies are applied (resilience, hysteresis phenomena). Such situations have not been investigated in food webs in this study, as we did one dynamical simulation per food web. Multiple simulations with different initial conditions are required for each food web to try to detect different asymptotic dynamics. Indeed, the existence of such dynamics cannot be fully determined using bifurcation analysis in systems with tens of state variables like food webs. However, quantifying uncertainties in food web model predictions in situations with external disturbances would be an important step toward more accurate model predictions.

Conclusion

We investigated the sensitivity of a food web model to the choice of functional response formulation (Holling's or Ivlev's FRs). We found a little effect of functional response formulation on food web persistence, whereas food web variability is significantly lower with Ivlev's functional response. Functional response slope at equilibrium explains this lower variability. In addition, food web variability is more driven by functional response formulation than by its parameter values and the tested range of trophic complexity (number of species, connectance). However, complexity-stability relationships are not qualitatively affected by functional response formulation. These conclusions are robust both to different combinations of model assumptions (cannibalism, primary production) and to different fitting methods to parameterize functional responses.

Because of intrinsic data variability and because model formulation is always a simplified representation of complex biological processes, the choice of functional response may remain uncertain for many species. The results demonstrate that this uncertainty in the formulation of a food web model can lead to uncertainties in the type of asymptotic dynamics it predicts. In addition, uncertainties in predicted system resilience in case of external disturbances are known to arise in some predator-prey models, and so are likely to occur in food web models. The quantification of these potential uncertainties in food web resilience may be a challenging way of research toward more accurate model predictions.

Figure 2 .Figure 3 .

 23 Figure 2. Food web dynamics as a function of trophic complexity predicted by Holling's (a, b) and Ivlev's functional responses (c, d). a, c: proportion of food webs with at least one species extinction. b, d: proportion of persistent food webs with a stable positive equilibrium. The number of food webs studied for each pair of parameters is enough to obtain 10 000 food webs with a stable positive equilibrium.

Figure 4 .

 4 Figure4. Elasticity of the functional response: estimated density in food webs reaching a positive equilibrium (a) and as a function of total prey biomass available (b) with Holling's (plain) or Ivlev's (dashed) functional response. Computed for predator species from 160 000 food webs of varying complexity (4 levels of connectance and 4 numbers of species, 10 000 food webs by pair of values). Density estimates are realized using non-parametric kernel methods[START_REF] Simonoff | Smoothing methods in statistics[END_REF].

Figure 5 .

 5 Figure5. Impact of the method used to fit the functional responses. a: Holling's disc equation (plain), Ivlev's functional response with the parameter set used in the study (bold dashed) and with parameter sets obtained with different fitting methods (thin dashed, see table2). b and c: food web dynamics (30 species, connectance of 0.15) predicted for these parameter sets. Gray labels and arrows summarize the main trend of results: a better fit at low prey biomass (a) leads to a closer predicted food web variability (c), but also to more distant functions at other prey biomass values (a) and more distant predicted food web persistence (b).

Table 1 .

 1 Parameter values used in the food web model. Parameter values from Heckmann et al.(2012) were estimated from empirical data sets (up to > 700 organisms from unicellular eukaryotes to plants and mammals, which span 20 orders of magnitude in body mass,[START_REF] Brose | Allometric scaling enhances stability in complex food webs[END_REF][START_REF] Brown | Toward a metabolic theory of ecology[END_REF] using allometric scaling or set to values similar to other studies for comparison (like[START_REF] Kartascheff | Positive complexity-stability relations in food web models without foraging adaptation[END_REF][START_REF] Kartascheff | Why allometric scaling enhances stability in food web models[END_REF].

	biological meaning	parameter	value	source	unit
	mortality rate	α	0.3	(Heckmann et al., 2012)	time -1
	per-capita competition rate	β	0.5	(Heckmann et al., 2012) biomass -1 time -1
	assimilation efficiency	λ	0.65	(Heckmann et al., 2012)	-
	resource biomass	B res	500	(Heckmann et al., 2012)	biomass
	magnitude of body mass range	x	8	(Plitzko et al., 2012)	-
	Holling's disc equation:
	attack rate	a H	6	(Heckmann et al., 2012) biomass -1 time -1
	handling time	h H	0.35	(Heckmann et al., 2012)	time
	Ivlev's functional response:
	maximal consumption rate	1/h I	1/0.36	see text	time -1
	satiation coefficient	a I h I	3.17 × 0.36	see text	biomass -1
	2.3 Parameterization and numerical study

Parameter values are displayed in table 1. The range of trophic complexity (20 to 60 species, connectance of 0.10 to 0.30) is set to include values used in previous studies

[START_REF] Gross | Generalized models reveal stabilizing factors in food webs[END_REF][START_REF] Kartascheff | Positive complexity-stability relations in food web models without foraging adaptation[END_REF][START_REF] Kartascheff | Why allometric scaling enhances stability in food web models[END_REF] Plitzko et al., 2012, among others)

for comparison. The number of simulations is limited by assuming that a species i has the same feeding preference f i,j = 1/|R i | for all its prey species j ∈ R i (with |R i | being the number of prey species of species i). This choice implies that T i is the average biomass of prey species j∈Ri B j / |R i | (weak generalist model, sensu

[START_REF] Williams | Effects of network and dynamical model structure on species persistence in large model food webs[END_REF]

.

Parameters in Ivlev's FR are chosen to minimize the weighted Euclidean distance between both functional responses:

Acknowledgments

The authors thank four anonymous reviewers for their very constructive comments on earlier versions of the manuscript. We also thank Bastien Bellemin-Noel for language proof-reading. The authors acknowledge the staff of the "Cluster de calcul intensif HPC" Platform (OSU Institut Pythéas, Aix-Marseille Université, INSU-CNRS) for providing computing facilities. The PhD scholarship of C.A. was funded by the French Ministry for Education and Research.

Appendix S1: the jacobian matrix

∂B j i, j = 1, ..., S, j = i