Angelo Luiz

Maxime Steffenel

Denis Martinasso

Trystram

Luiz Angelo

Steffenel Crestic-Syscom

Maxime Martinasso
email: maxime.martinasso@imag.fr

Denis Trystram
email: denis.trystram@imag.fr

Assessing Contention Effects of All-to-All Communications on Clusters and Grids

Keywords: total exchange, network contention, performance modeling, computational

One of the most important collective communication patterns used in scientific applications is the complete exchange, also called All-to-All. Although efficient algorithms have been studied for specific networks, general solutions like those available in wellknown MPI distributions (e.g. the MPI_Alltoall operation) are strongly influenced by the congestion of network resources. In this paper we present an integrated approach to model the performance of the Allto-All collective operation, which consists in identifying a contention signature that characterizes a given network environment, using it to augment a contentionfree communication model. This approach, assessed by experimental results, allows an accurate prediction of the performance of the All-to-All operation over different network architectures with a small overhead. We also discuss the problem of network contention in a grid environment, studying some strategies to minimize the impact of contention on the performance of an Allto-All operation.

Index Terms-MPI, all-to-all,

I. Introduction

One of the most important collective communication patterns for scientific applications is the total exchange [START_REF] Christara | An efficient transposition algorithm for distributed memory computers[END_REF] (also called All-to-All), in which each process holds n different data items that should be distributed among the n processes, including itself. An important example of this communication pattern is the All-to-All operation, where all messages have the same size m.

Although efficient All-to-All algorithms have been studied for specific networks structures like meshes, hypercubes, tori and circuit-switched butterflies, general solutions like those found in well-known MPI distributions rely on direct point-to-point communications among the processes. Because all communications are started simultaneously, architecture independent algorithms are strongly influenced by the saturation of network resources and subsequent loss of packets -the network contention.

In the first part of this paper we present a new approach to model the performance of the All-to-All collective operation. Our strategy consists in identifying a contention signature that characterizes a given network environment. Using such contention signature, we are able to accurately predict the performance of the All-to-All operation, with an arbitrary number of processes and message sizes. To demonstrate our approach, we present experimental results obtained with different network architectures (Fast Ethernet, Gigabit Ethernet and Myrinet). We believe that this model can be extremely helpful on the development of application performance prediction frameworks such as PEMPIs [START_REF] Edson | PEMPIs: A new metodology for modeling and prediction of MPI programs performance[END_REF] and task schedulers.

Later, we focus on the performance modeling of Allto-All in a grid environment. We observe that contention at the wide-area level not only exists but represents a major component on the communication performance. We propose therefore an algorithm adapted for intercluster message exchanges. This algorithm minimizes the impact of network contention and, by extension, allows a better performance modeling. Indeed, we show that this algorithm not only outperforms the traditional Allto-All algorithms on a grid environment but also allows an accurate performance modeling using only elements derived from our local-area models.

This paper is therefore organized as follows: Section II presents a survey of performance modeling under communication contention. Section III presents the network models used in this paper and we formalize the total exchange problem, as well as some performance lower bounds. In Section IV we propose a strategy to characterize the contention signature of a given network and for instance, to predict the performance of the All-to-All operation. Section V validates our model against experimental data obtained on different network architectures (Fast Ether-net, Gigabit Ethernet and Myrinet). In Section VI we provide a study case for predicting the performance of the All-to-All algorithm on a grid platform. Finally, Section VII presents some conclusions and the future directions of our work.

II. Related Works

In the All-to-All operation, every process holds m × n data items that should be equally distributed among the n processes, including itself. Because general implementations of the All-to-All collective communication rely on direct point-to-point communications among the processes the network can easily become saturated, and by consequence, degrades the communication performance. Indeed, Chun and Wang [START_REF] Tam | Performance Studies of High-Speed Communication on Commodity Cluster[END_REF][6] demonstrated that the overall execution time of intensive exchange collective communications are strongly dominated by the network contention and congestive packet loss, two aspects that are not easy to quantify. As a result, a major challenge on modeling the communication performance of the All-to-All operation is to represent the impact of network contention.

Unfortunately, most communication models like those presented by Christara et al. [START_REF] Christara | An efficient transposition algorithm for distributed memory computers[END_REF] and Pjesivac-Grbovic et al. [START_REF] Pjesivac-Grbovic | Performance analysis of MPI collective operations[END_REF] do not take into account the potential impacts of network contention. Indeed, these works usually represent the All-to-All operation as parallel executions of the personalized one-to-many pattern [START_REF] Johnssonn | Optimum broadcasting and personalized communication in hypercubes[END_REF], as presented by the linear point-to-point model below, where α is the start-up time (the latency between the processes), 1 β is the bandwidth of the link, m represents the message size in bytes and n corresponds to the number of processes involved in the operation:

T = (n -1) × (α + βm) (1)
The development of contention-aware communication models is relatively recent, as shown by Grove [START_REF] Grove | Performance Modelling of Message-Passing Parallel Programs[END_REF]. For instance, Adve [START_REF] Vikram | Analysing the Behavior and Performance of Parallel Programs[END_REF] presented one of the first models to take into account the effects of resource contention. This model considers that the total execution time of parallel programs is the sum of four components, namely:

T = t computation + t communication +t resource-contention + t synchronization (2)
While conceptually simple, this model was non-trivial in practice because of the non-deterministic nature of resource contention, and because of the difficulty to estimate average synchronization delays.

While the non-deterministic behavior of the network contention is a major obstacle to modeling communication performance, some authors suggested a few techniques to adapt the existing models. As consequence, Bruck [START_REF] Bruck | Efficient algorithms for all-to-all communications in multiport message-passing systems[END_REF] suggested the use of a slowdown factor to correct the performance predictions. Similarly, Clement et al. [START_REF] Clement | Network performance modelling for PM clusters[END_REF] introduced a technique that suggested a way to account contention in shared networks such as non-switched Ethernet, consisting in a contention factor γ that extends the linear communication model T:

T = α + β × m × γ ( 3
)
where γ is equal to the number of processes. A restriction on this model is that it assumes that all processes communicate simultaneously, which is only true for a few collective communication patterns. Anyway, in the cases where this assumption holds, they found that this simple contention model enhanced the accuracy of their predictions for essentially zero extra effort.

The use of a contention factor was supported by the work of Labarta et al. [START_REF] Labarta | DiP: A parallel program development environment[END_REF], that intent to approximate the behavior of the network contention by considering that if there are m messages ready to be transmitted, and only b available buses, then the messages are serialized in m b communication waves. Also, König et al. [START_REF] König | Analysis of gossiping algorithms with restricted buffers[END_REF] have shown indeed that some All-to-All algorithms that are optimal with unlimited buffers become less efficient when communications depend on restricted buffers sizes.

A similar approach was followed by Jeannot et al. [START_REF] Jeannot | Two fast and efficient message scheduling algorithms for data redistribution through a backbone[END_REF], who designed scheduling algorithms for data redistribution through a backbone. In their work, they suppose that at most k communications can be performed at the same time without causing network contention (the value of k depending on the characteristics of the platform). Using the knowledge of the application transfer pattern, they proposed two algorithms to schedule the messages transfers, performing an application-level congestion control that in most cases outperforms the TCP contention control mechanism.

Most recently, some works aimed to design contentionaware performance models. For instance, LoGPC [START_REF] Csaba | LoGPC: Modeling network contention in message-passing programs[END_REF] presents an extension of the LogP model that tries to determine the impact of network contention through the analysis of k-ary n-cubes. Unfortunately, the complexity of this analysis makes too hard the application of such model in practical situations.

Another approach to include contention-specific parameters in the performance models was introduced by Chun [START_REF] Tam | Performance Studies of High-Speed Communication on Commodity Cluster[END_REF]. In his work, the contention is considered as a component of the communication latency, and by consequence, the model uses different latency values depending on the message size. Although easier to use than LoGPC, this model does not take into account the number of messages passing in the network nor the link capacity, which are clearly related to the occurrence of network contention.

III. Definitions

A. Network Models

In this work we assume that the network is fully connected, which corresponds to most current parallel machines with distributed memory.

Communication Model: The links between pairs of processes are bidirectional, and each process can transmit data on at most one link and receive data on at most one link at any given time.

Transmission Model: We use Hockney's notation [START_REF] Hockney | The communication challenge for MPP: Intel paragon and meiko cs-2[END_REF] to describe our transmission model. Therefore, the time to send a message of size w i,j from a process p i to another process p j , is α + w i,j β, where α is the start-up time (the communication latency between the processes) and 1 β is the bandwidth of the link. As in this paper we assume that all links have the same latency and bandwidth, and because we only investigate the regular version of the MPI_Alltoall operation where all messages have the same size m, ∀i, ∀j, w i,j = m, and therefore the time to send a message from a process p i to a process p j is α + mβ.

Synchronization Model: We assume an asynchronous communication model, where transmissions from different processes do not have to start at the same time. However, all processes start the algorithm simultaneously. This synchronization model corresponds to the execution of the MPI_Alltoall operation, used as reference in this work.

B. Problem Definitions

In the total exchange problem, n different processes hold each one n data items that should be evenly distributed among the n processes, including itself. Because each data item has potentially different contents and sizes according to their destinations, all processes engage a total exchange communication pattern. Therefore, a total exchange operation will be complete only after all processes have sent their messages to their counterparts, and received their respective messages.

Formally, the total exchange problem can be described using a weighted digraph dG(V, E) of order n with V = {p 0 , ..., p n-1 }. This digraph is called a message exchange digraph or MED for short. In a MED, the vertices represent the process nodes, and the arcs represent the messages to be transmitted. An integer w(e) is associated with each arc e = (p i , p j ), representing the size of the message to be sent from process p i to process p j . Note that there is not necessarily any relationship between a MED and the topology of the interconnection network.

The port capacity of a process for transmission is the number of other processes to which it can transmit simultaneously. Similarly, the port capacity for reception is the number of other processes from which it can receive simultaneously. We will concentrate on the performance modeling problem with all port capacities restricted to one for both transmitting and receiving. This restriction is well-known in the literature as 1-port full-duplex.

C. Notation and lower bounds

In this section, we present theoretical bounds on the minimum number of communications and on the bandwidth for the general message exchange problem. The number of communications determines the number of start-ups, and the bandwidth depends on the message weights.

Given a MED dG(V ; E), we denote the in-degree of each vertex p i ∈ V by ∆ r (p i ), and the out-degree by ∆ s (p i ). Let ∆ r = max pi∈V {∆ r (p i )} and ∆ s = max pi∈V {∆ s (p i )}. Therefore, we obtain the following straightforward bound on the number of start-ups.

Claim 1. The number of start-ups needed to solve a message exchange problem on a digraph dG(V ; E) without message forwarding is at least max(∆ s , ∆ r ).

Given a MED dG(V, E), the bandwidth bounds are determined by two obvious bottlenecks for each vertexthe time for it to send its messages and the time for it to receive its messages. Each vertex p i has to send messages with sizes {w i,j | j = 0 . . . n -1}. The time for all vertices to send their messages is at least t s = max i n-1 j=0 w i,j β. Similarly, the time for all vertices to receive their messages is at least t r = max j n-1 i=0 w i,j β. Claim 2. The time to complete a personalized exchange is at least max{t s , t r }.

We can combine the claims about the number of startups and the bandwidth when message forwarding is not allowed.

Claim 3. If message forwarding is not allowed, and either the model is synchronous or both maxima are due to the same process, the time to complete a personalized exchange is at least max(∆ s , ∆ r ) × α + max{t s , t r }.

Because in this paper we do not assume messages forwarding, the fan-in and fan-out of a process must be (n -1). Further, as we consider messages to be the same size and the network to be homogeneous, we can simplify Claim 3 so that the following bound holds.

Proposition 1. If message forwarding is not allowed, and all messages have size m, and both bandwidth and latency are identical to any connection between two different processes p i and p j , the time to complete a total exchange is at least

(n -1) × α + (n -1) × βm.
Proof. The proof is trivial, as the time to complete a total exchange is at least the time a single process needs to send one message to each other process.

IV. Contention Signature Approach

To cope with this problem and to model the contention impact on the performance of the All-to-All operation, we adopt an approach similar to Clement et al. [START_REF] Clement | Network performance modelling for PM clusters[END_REF], which considers the contention sufficiently linear to be modeled. Our approach, however, tries to identify the behavior of the Allto-All operation with regard to the theoretical lower bound (Proposition 1) on the 1-port communication model. In our hypothesis, the network contention depends mostly on the physical characteristics of the network (network cards, links, switches), and consequently, the ratio between the theoretical lower bound and the real performance represents a "contention signature" of the network. Once identified the signature of a network, it can be used in further experiments to predict the communication performance, provided that the network infrastructure does not change.

Initially, we consider communication in a contention-free environment. In this case, a process that sends messages of size m to n-1 processes needs at least (n-1)×α+(n-1) × mβ time units. Further, by the properties of the 1port communication model, the total communication time of the All-to-All operation must be at least (n -1) × α + (n-1)×mβ time units if all processes start communicating simultaneously, as stated by Proposition 1.

In the case of the All-to-All operation, however, the intensive communication pattern tends to saturate the network, causing message delays and packet loss that strongly impact on the communication performance of this collective communication. In this network congestion situation, traditional models such as those presented by Christara [START_REF] Christara | An efficient transposition algorithm for distributed memory computers[END_REF] do not hold anymore, even if the communication pattern has not changed. Therefore, our approach to model the performance of the MPI_Alltoall operation despite network contention consists on determining a contention ratio γ that express the relationship between the theoretical performance (lower bound) and the real completion time. For simplicity, we consider that this contention ratio γ is constant and depends exclusively on the network characteristics. Therefore, the simplest way to integrate this contention ratio γ in our performance model would be as follows:

T = (n -1) × (α + mβ × γ) (4)

A. Non-linear aspects

Although the performance model augmented by use of the contention ratio γ improves the accuracy of the predictions, we observe nonetheless that some network architectures are still subject to performance variations according to the message size. To illustrate this problem, we present in Fig. 1, a detailed mapping of the communication time of the MPI_Alltoall operation in a Gigabit Ethernet network. We observe that the communication time does not increase linearly with the message size, but instead, present a non-linear behavior that prevents our model to accurately predict the performance when dealing with small messages.

To cope with this non-linearity, we propose an extension of the contention ratio model to better represent this phenomenon when messages are sufficiently large. Hence, we augment the model with a new parameter δ, which depends on the number of processes but also on a given message size M . As a consequence, the association of different equations helps to define a more realistic performance model for the MPI_Alltoall operation, as follows:

T = (n -1) × (α + mβ × γ) if m < M (n -1) × (α + mβ × γ + δ) if m ≥ M (5)
V. Validation To validate the approach proposed in this paper, this section presents our experiments to model the performance of MPI_Alltoall operation using three network architectures, Fast Ethernet, Gigabit Ethernet and Myrinet. As previously explained, our approach compares the expected and real performance of the MPI_Alltoall operation using a sample experiment with n ′ nodes; the relationship between these two measures allows us to define the γ and δ parameters that characterize the "network contention signature".

To obtain these parameters, we compare the sample data obtained from both theoretical lower bound and experimental measure, while varying the message size. Indeed, the lower bound comes from Proposition 1, with parameters α and β obtained from a simple point-to-point measure. The parameters γ and δ are obtained through a linear regression with the Generalized Least Squares method, comparing at least four measurement points in order to better fit the performance curve.

The different experiments presented in this paper represent the average of 100 measures for each set of parameters (message size, number of processes), and were conducted over two clusters of the Grid'5000 network1 :

The The GdX (GriD'eXplorer) cluster, operated by INRIA-Futurs. This cluster includes 216 nodes with dual AMD Opteron processors at 2 GHz running Debian Linux kernel 2.6.8 and a Broadcom Gigabit Ethernet network.

A. Fast Ethernet

Taking as basis the measured performance for a 24 machines network, we were able to approximate the performance of the Fast Ethernet network with a contention ratio γ = 1.0195. Indeed, this relatively small difference must be considered in the light of the retransmission policy: although the communication latency (and therefore the timeouts) is relatively small (around 60 µs), the reduced bandwidth of the links minimizes the impact of the retransmission of a lost packet. More important, we observe that the experimental measures behave like an affine equation, showing a start-up cost usually not considered by the traditional performance model which corresponds to the δ parameter proposed in our model. Therefore, we determined δ = 8.23 ms for messages larger than M = 2 kB, which means that each simultaneous communication induces an overload of 8.23 ms to the completion time of the All-to-All operation. Applying both γ and δ parameters to the performance model we were able to approximate our predictions from the performance of the MPI_Alltoall operation with an arbitrary number of processes, as demonstrate in Fig. 2a. We observe indeed that our error rate is usually smaller than 10% when there are enough processes to saturate the network, as presented in Fig. 2b.

B. Gigabit Ethernet

To compute the contention ratio γ and a start-up cost δ, we use sample data for an arbitrary number of processes. Indeed, we chose in this example the results for an execution of the All-to-All operation with 40 processes (one by machine). Using linear regression on these data we obtain γ = 4.3628 and δ = 4.93 ms (to be used only for messages larger than M = 8 kB). As a result, the performance predictions from our model correspond to the curve presented on Fig. 3a. As in the case of the Fast Ethernet network, the error rate is quite small when the network becomes saturate, even when we consider different message sizes (Fig. 3b).

C. Myrinet

Although the two previous experiments give important proofs on the validity of our modeling method, they share many similarities on both network architecture and transport protocol (TCP/IP). To ensure that our method is not bounded to a specific infrastructure, we chose to validate our performance model also in a Myrinet network, using the gm transport protocol. Because of the Myrinet+gm stack differs considerably from the Ethernet+TCP/IP stack, any systematic behavior introduced into our sampling data by these architectures should be exposed. Indeed, the Myrinet network differs from Ethernetbased architectures due to an start-up cost almost inexistent (one of the main characteristics of the Myrinet+gm stack). Indeed, we were able to fit the performance of a 24-processes All-to-All operation using only the contention ratio γ = 2, 49754 (as the linear regression pointed a startup cost δ smaller than 1 microsecond).

Nevertheless, when applying this factor to an arbitrary number of machines, as presented in Fig. 4a, we observe that our predictions do not follow the experimental data as observed before with Fast Ethernet and Gigabit Ethernet. Actually, a close look at the error rate (Fig. 4b) indicates that network saturation occurs only when there are more then 40 communicating processes (evidenced by the constant error rate from that point). These results demonstrate the limitations of our approach: while a contention ratio may provide precise performance predictions, it depends on the data used to define the network signature. By using reference data from a partially saturated network we are subjected to inaccurate approximations (even if they are better than the contention unaware predictions).

VI. Problem of Total Exchange between Two Clusters

Multicluster and computational grids are popular heterogeneous environments used by the HPC community, and part of our efforts concentrate on the development of efficient algorithms and performance models for these environments. More specifically, we must deal with communications costs that differ significantly if they are local or remote. Without loss of generality, let us assume two clusters C 1 and C 2 with respectively n 1 nodes and n 2 nodes, as represented in Figure 5. A network, called a backbone, interconnects the two clusters. We assume that a cluster use the same network card to communicate to one of its node or to a node of another cluster. Based on that topology inter cluster communications are never faster than communication within a cluster. Let us suppose that an application is running and using both clusters (for example, a code coupling application). One part of the computation is performed on cluster C 1 and the other part on cluster C 2 . During the application, data must be exchanged from C 1 to C 2 using the alltoall pattern. Altogether, this means that we will have to transfer (n 1 + n 2 ) 2 messages over different network environments. The data of all these messages are different Measured and predicted performance for the standard MPI_Alltoall in a grid but the size of the messages are the same and is given and called m (in bytes). Several MPI libraries (OpenMPI, MPICH2, etc.) implement the allltoall routine assuming that all the nodes are on the same clusters, which means that all communications have the same weight. However, in our case, some messages are transferred within a cluster (from a node of C 1 to a node of C 1 or from C 2 to C 2 ) or between the two clusters.

From the homogeneous cluster prediction model presented in the first part of this paper, which allows a quite accurate representation on the performance of localarea networks, we initially proceeded our multi-cluster performance modelling by composing local (contentionaware) and remote communications. Estimate completion time corresponds to the sum of the contention-corrected local-area cost (T Cn ) and the wide-area communication cost (a function of wide-area latency α w and bandwidth β w ) :

T = max(T C1 , T C2 ) + max(n 1 , n 2 ) × (α w + β w × m) (6)
Unfortunately, this simple strategy fails to represent the operation of the MPI_Alltoall in a grid. Figure 6 compares estimate and measured completion time of the MPI_Alltoall implementation from OpenMPI in a grid with two clusters of 30 machines each. From this experience we observe that not only contention affects also wide-area links but also that local area communications play a small role in the overall performance. Of course, one could try to define additional parameters for the wide-area communications, but the final model would be too complex to be useful in real situation. Instead, we decided to address this problem by redefining the All-to-All problem against the challenges that characterize a grid environment.

A. Minimizing the impact of contention on the backbone

When dealing with wide-area networks, the most important factor to be considered is the time a message takes to be delivered. Indeed, in addition to the geographical distance, message are subjected to network protocols heterogeneity, message routing and transient interferences on the backbone.

Actually, popular algorithms for collective communications on grids (such as the ones implemented in PACX MPI [START_REF] Gabriel | Distributed Computing in a Heterogenous Computing Environment[END_REF] , GridMPI [10] and MagPIe [START_REF] Kielmann | Network performance-aware collective communication for clustered wide area systems[END_REF]) try to minimize communications over the wide-area network by defining a single coordinator in every cluster, which participates in the inter-cluster data transfers across the wide-area backbone. By minimizing the number of WAN communication steps, we reduce the probability of inducing contention and accumulating transmission delays on the messages.

However, a single communication between each cluster is an approach inappropriate for the MPI_Alltoall operation. First, it induces additional communication steps to/from the cluster coordinator, which becomes a bottleneck. Second, this approach is not optimal concerning the usage of the wide-area bandwidth, as wide-area backbones are designed to support simultaneous transfers [START_REF] Casanova | Network modeling issues for grid application scheduling[END_REF]. Hence, in order to improve the performance in a WAN, we need to change the MPI_Alltoall algorithm strategy.

B. The Local Group(LG) algorithm

To cope with this problem, we try to minimize wide-area communication steps in a different way. Actually, most of the complexity of the All-to-All problem resides on the need to exchange different messages through different networks (local and distant). The traditional implementation of the MPI_Alltoall operation cannot differentiate these networks, leading to poor performances. However, if we assume that communications between clusters are slower than intra-clusters ones, it might be useful to collect data in the local level before sending it in parallel through the backbone, in a single communication step.

As a consequence, we propose in [START_REF] Jeannot | Fast and efficient total exchange on two clusters[END_REF] a grid-aware solution which performs on two phases. In the first phase only local communications are performed. During this phase the total exchange is performed on local nodes on both cluster and extra buffers are prepared for the second (inter-cluster) phase. During the second phase data are exchanged between the clusters. Buffers that have been prepared during the first phase are sent directly to the corresponding nodes in order to complete the total exchange.

More precisely, our algorithm called Local Group or simply LG works as follow. Without loss of generality, let us assume that cluster C 1 has less nodes than C 2 (n 1 ≤ n 2 ).

Nodes are numbered from 0 to n 1 + n 2 -1, with nodes from 0 to n 1 -1 being on C 1 and nodes from n 1 to n 1 +n 2 -1 being on cluster C 2 . We call M i,j the message (data) that has to be sent from node i to node j. The phases are sumup in Algorithm 1.

1) First phase: During the first phase, we perform the local exchange: Process i sends M i,j to process j, if i and j are on the same cluster. Then it prepares the buffers for the remote communications. On C 1 data that have to be send to node j on C 2 is first stored to node j mod n 1 . Data to be sent from node i on C 2 to node j on C 1 is stored on node ⌊i/n 1 ⌋ × n 1 + j.

= i + n 1 × s on C 2 (if j < n 1 + n 2 ). More precisely i sends M k,j to j where k ∈ [0, n 1 ] and j sends M k,i to i where k ∈ [n 1 × s, n 1 × s + n 1 -1].
3) Example: Suppose that n 1 = 3 and n 2 = 7. What happens to the message M 7,2 (i.e the messages that goes from node 7 on cluster C 2 to node 2 on C 1 )? This is illustrated in Figure 7. During the first phase it is stored on node ⌊7/3⌋ × 3 + 2 = 8 on C 2 . Then during the second phase it is sent to node 2 during the step s = 2: node 8 sends M 6,2 , M 7,2 and M 8,2 to node 2. During this step nodes 1 and 7 and node 0 and 6 exchange data as well, while in the previous step node 0 and 3, 4 and 1 and 5 and 2 exchange data and in the last step only 0 and 9 exchange data.

C. Comparison with the standard Total Exchange algorithm

As our algorithm tries to minimize the number of intercluster communications between the clusters, we need only 2 × max(n 1 , n 2 ) messages of size m × min(n 1 , n 2 ) in both directions against 2 × n 1 × n 2 messages of size m in the traditional algorithm. For instance, the exchange of data between two clusters with the same number of processes will proceed in one single communication step of the second phase. At the other hand, if n 2 ≫ n 1 , the total number of communication steps will be similar to the traditional algorithm.

An important difficulty that arises in the implementation of the algorithm is how to determine if a node belong to the first cluster or to the second cluster. Indeed, this information is the only extra information required to run the algorithm compared to the standard one. We propose several approaches:

• a simple solution is to extend the MPI API by adding a multi-cluster alltoall call where the users passes the value of n 1 and n 2 assuming that the first n 1 nodes of the machine file belong to cluster 1 and the n 2 nodes to cluster 2. The problem of such an approach is to port old MPI code to the multi-cluster version. • At the beginning of the execution MPI can run some network benchmark in order to try to guess the topology of the infrastructure. This is what already does OpenMPI in the case of multiple network interfaces. • GridMPI has an interesting approach that consists in launching a different MPI "client" on each cluster, each client being interconnected to the other through a GridMPI server. Based on that, for each MPI collective communication call GridMPI is able to determine the cluster a node belongs to. • another solution consists in modifying the MPI machine file by adding keywords that tells to which cluster a node belongs to, such as in MagPIe [START_REF] Kielmann | Network performance-aware collective communication for clustered wide area systems[END_REF]. As our algorithm minimizes the number of inter-cluster communications, it is also wide-area optimal since it ensures that a data segment is transferred only once between two clusters separated by a wide-area link. Additionally, wide-area transmissions pack several messages together, reducing the impact of transient interferences on the backbone. Hence, Figure 8 presents a comparison between the traditional algorithm used by OpenMPI and the LG algorithm. We observe that LG improves the performance of the MPI_Alltoall operation, reaching over than 50% of performance improvement comparing to the traditional strategy.

D. Modeling approach

As shown above, the algorithm we propose to optimize All-to-All communications in a grid environment rely on the relative performances of both local and remote networks. Indeed, we extend the total exchange among nodes in the same cluster in order to reduce transmissions through the backbone. This approach has two consequences for performance prediction: first, it prevents contention in the wide-area links, which are hard to model. Second, the transmission of messages packed together is less subjected to network interferences. For instance, we can design a performance model by composing local-area predictions obtained with our contention ratio model and wide-area predictions that can be easily obtained from traditional methods. Hence, an approximate model would consider the following parts, where T Cn corresponds to Equation 5:

T = max(T C1 , T C2 ) + ⌈n 2 /n 1 ⌉ × (α w + β w × m × n 1 ) (7)

E. Experimental validation

To validate the algorithm we propose in this paper, this section presents our experiments to evaluate the performance of the MPI_Alltoall operation with two clusters connected through a backbone.

These experiments were conducted over two clusters of the Grid'5000 platform2 , one located in Nancy and one located in Rennes, approximately 1000 Km from each other. Both clusters are composed of identical nodes (dual Opteron 246, 2 GHz) locally connected by a Gigabit Ethernet network and interconnected by a private backbone of 10 Gbps. All nodes run Linux, with kernel 2.6.13 and OpenMPI 1.1.4. The measures were obtained with the broadcast-barrier approach [START_REF] Bronis | Accurately measuring MPI broadcasts in a computational grid[END_REF].

To model the communication performance of both intercluster and intra-cluster communications we use the parameterised LogP model (pLogP) [START_REF] Kielmann | Network performance-aware collective communication for clustered wide area systems[END_REF]. The pLogP parameters for both local and distant communications were obtained with the method described in [START_REF] Kielmann | Fast measurement of LogP parameters for message passing platforms[END_REF]. To model the contention at the local level we used γ = 2.6887 and δ = 0.005039 for M >= 1KB parameters obtained from the method of the least squares as described in [START_REF] Luiz | Modeling network contention effects on alltoall operations[END_REF].

Therefore, in Figure 9 we compare the performance predictions obtained with Equation 7 against the effective completion time of the LG algorithm. We observe that prediction fit with a good accuracy to the real MPI_Alltoall performance prediction on a 30+30 grid

LG Algorithm Contention-aware estimation Fig. 9. Performance predictions for the LG algorithm execution times, which is not possible with the traditional MPI_Alltoall algorithm. Indeed, the new algorithm minimizes the impact of distant communications, concentrating the contention problems at the local level. Because we are able to predict the performance of local communications even under contention, we can therefore establish an accurate performance model adapted to grid environments.

VII. Conclusions and Future Works

In this paper we address the problem of modeling the performance of Total Exchange communication operations, usually subject to important variations caused by network contention. Because traditional performance models are unable to predict the real completion time of an Allto-All operation, we try to cope with this problem by identifying the contention signature of a given network. In our approach, two parameters γ and δ are used to augment a linear performance model in order to fit the performance of the MPI_Alltoall operation. Because these parameters characterize the network contention and are independent of the number of communicating processes, they can be used to accurately predict the communication performance when communications tend to saturate the network. Indeed, we demonstrate our approach through experiments conducted on popular network architectures, Fast Ethernet, Gigabit Ethernet and Myrinet.

We also address the problem of modeling the MPI_Alltoall operation in a grid environment represented by two clusters. We demonstrate that network contention also affects wide-area communications but, instead of developing a new complex performance model, we observed the behaviour of MPI_Alltoall and proposed a gridaware algorithm that reduces drastically the impact of contention, improving drastically the communication performance. Additionally, we are able to predict the overall performance of this new algorithm by simply extending our local-area contention-aware model, which represents an important advantage when comparing with the traditional MPI_Alltoall implementation.

Our efforts now are concentrating on the development of grid-optimized algorithms for the generalized total exchange problem (represented by the MPI_Alltoallv op-eration), where processes exchange messages of different sizes. By developing algorithms that minimize the network contention we may be able to improve the performance of the operation at the same time as we allow a simple and accurate performance modeling.

Fig. 1 .

 1 Fig. 1. Non-linearity of communication cost with small messages

  icluster2 cluster, located at INRIA-Rhone-Alpes, composed of 104 dual Itanium2 nodes at 900 MHz, used for the experiments with the Fast Ethernet network (5 Fast Ethernet switches -20 nodes per switch -interconnected by 1 Gigabit Ethernet switch) and the Myrinet 2000 network (one 128 ports M3-E128 Myrinet switch). All machines run Red Hat Enterprise Linux AS release 3, with kernel version 2.4.21.

Fig. 2 .

 2 Fig. 2. Performance prediction on a Fast Ethernet network

Fig. 3 .

 3 Fig. 3. Performance prediction on a Gigabit Ethernet network

Fig. 4 .Fig. 5 .

 45 Fig. 4. Performance prediction on a Myrinet network

  Fig.[START_REF] Tam | Realistic communication model for parallel computing on cluster[END_REF].Measured and predicted performance for the standard MPI_Alltoall in a grid

Fig. 7 .

 7 Fig. 7. Example of the 2 phases of the algorithm

Fig. 8 .

 8 Fig. 8. Performance comparison between OpenMPI and LG algorithms

Algorithm 1

 1 The LG (Local Group) algorithm when n 1 ≤ n 2 // Local Phase for i = {0, ..., (n 1 + n 2 ) -1} do in parallel for j = {0, ..., (n 1 + n 2 ) -1} do if i < n 1 // the sender is on C 1 send M i,j to j mod n 1 else // the sender is on C 2 if j ≥ n 1 // the receiver is on C 2 send M i,j to j else // the receiver is on C 1 send M i,j to ⌊i/n 1 ⌋ × n 1 + j // Inter-cluster Phase for s = {1, ..., ⌈n 2 /n 1 ⌉} for i = {1, ..., n 1 -1} do in parallel if (i + s × n 1 < n 1 + n 2 ) exchange messages between i and j = i + s × n 1 2) Second phase: During the second phase only n 2 intercluster communications occurs. This phase is decomposed in ⌈n 2 /n 1 ⌉ steps with at most n 1 communications each. Steps are numbered from 1 to ⌈n 2 /n 1 ⌉ During step s node i of C 1 exchange data stored in its local buffer with node j

http://www.grid5000.fr/

http://www.grid5000.fr/

Acknowledgments

We are grateful to the anonymous referees for many comments and helpful suggestions which helped us improve the focus of the paper.

Experiments presented in this paper were carried out using the Grid'5000 experimental testbed, an initiative from the French Ministry of Research through the ACI GRID incentive action, INRIA, CNRS and RENATER and other contributing partners (see https://www.grid5000.fr).


Daniel Montoya
email: daniel.montoya@sete.cnrs.fr

Bart Haegeman

Sabrina Gaba

Claire De Mazancourt

Michel Loreau

Habitat fragmentation and food security in crop pollination systems

Keywords: 75, Figures: 5, Tables: 1 Habitat fragmentation, pollination, biodiversity, food security, stability, ecosystem services, global change

come   L'archive ouverte pluridisciplinaire

Introduction

Ensuring stable food supplies for a growing population is one of the UN Sustainable Development Goals (UN2017). Biodiversity-based approaches to food security suggest that safeguarding species diversity is fundamental to increase yields and guarantee stable yields. This is supported by theory and data reporting positive effects of biodiversity on the mean values of various ecosystem functions and services, such as biomass production, nutrient cycling and crop pollination [START_REF] Lefcheck | Biodiversity enhances ecosystem multifunctionality across trophic levels and habitats[END_REF][START_REF] Woodcock | Meta-analysis reveals that pollinator functional diversity and abundance enhance crop pollination and yield[END_REF].

Additionally, biodiversity can stabilize ecosystem service supply by providing an insurance against environmental fluctuations [START_REF] Loreau | Biodiversity as spatial insurance in heterogeneous landscapes[END_REF], which are predicted to intensify under global change [START_REF] Giorgi | Regional climate information-evaluation and projections[END_REF][START_REF] Saltz | The impact of increased environmental stochasticity due to climate change on the dynamics of Asiatic wild ass[END_REF][START_REF] Fischer | Robust spatially aggregated projections of climate extremes[END_REF]. Biodiversity insurance effects have been observed in agriculture, where a greater diversity of crops is associated with increased year-to-year temporal stability of total production (Winfree & Kremen 2009;[START_REF] Renard | National food production stabilized by crop diversity[END_REF]. However, most biodiversity in intensively-managed agricultural landscapes is found in the remaining fragments of (semi)natural habitat not converted into crop land, and the effects of such non-crop biodiversity on the provision and stability of crop production are not clearly understood. This has led to a growing concern over the large-scale conversion of natural habitats into crop land and their effects on biodiversity and food production [START_REF] Godfray | Food security and sustainable intensification[END_REF]; UN2017; [START_REF] Réchauchère | Sustainable Agriculture Reviews 30: Environmental Impact of Land Use Change in Agricultural Systems[END_REF].

Agricultural systems are often highly fragmented with areas of intensive cultivation interspersed among remnant patches of seminatural habitat. This loss and fragmentation negatively affect biodiversity and many ecosystem functions and services [START_REF] Haddad | Habitat fragmentation and its lasting impact on Earth's ecosystems[END_REF][START_REF] Rybicki | Habitat fragmentation and species diversity in competitive communities[END_REF]. Despite this, fragments of natural habitat continue to supply important services. The spatial coexistence of crops and natural land creates an opportunity for spillover effects [START_REF] Rand | Spillover edge effects: the dispersal of agriculturally subsidized insect natural enemies into adjacent natural habitats[END_REF], a situation where ecological interactions extend across habitats boundaries and propagate ecological functions. In some cases, fragmentation can increase ecosystem service supply [START_REF] Martin | The interplay of landscape composition and configuration: new pathways to manage functional biodiversity and agroecosystem services across Europe[END_REF][START_REF] Haan | Predicting landscape configuration effects in agricultural pest suppression[END_REF]), e.g. if fragmentation of natural habitat for pollinators optimizes interspersion with crop land to maximize crop pollination [START_REF] Brosi | Optimal design of agricultural landscapes for pollination services[END_REF]). But fragmentation can also reduce ecosystem service supply if biodiversity decreases significantly in the remnant fragments of natural habitat [START_REF] Haddad | Habitat fragmentation and its lasting impact on Earth's ecosystems[END_REF]. For most services, however, we do not know how fragmentation affects their provision in fragmented landscapes, and this limits our capacity to manage biodiversity-based ecosystem service provision, e.g. crop pollination, and food security in human-dominated landscapes.

Recent research has revealed strong, non-linear effects of land-use change on crop production at multiple spatial scales. For example, theoretical studies agree on the humpshaped relationship between animal pollination-dependent crop production and the fraction of remnant natural land within intensive farming systems [START_REF] Braat | The Cost of Policy Inaction: The case of not meeting the 2010 biodiversity target[END_REF][START_REF] Mitchell | Strong and non-linear effects of fragmentation on ecosystem service provision at multiple scales[END_REF][START_REF] Montoya | Trade-offs in provisioning and stability of multiple ecosystem services in agroecosystems[END_REF]. Empirical research on the stability of animal pollination dependent crop production shows that stability decreases with agricultural intensification and the degree to which crops depend on animal pollination (Garibaldi et al. 2011a(Garibaldi et al. , b, 2014;;[START_REF] Deguines | Largescale trade-off between agricultural intensification and crop pollination services[END_REF]. Changes in crop production stability of animal-pollinated crops also depend on the spatial composition and structure of agricultural landscapes, such as the amount of natural land cover [START_REF] Montoya | Trade-offs in provisioning and stability of multiple ecosystem services in agroecosystems[END_REF] and the isolation of crops from natural land (Garibaldi et al. 2011b). Overall, both theoretical and empirical studies suggest that improved management of agricultural landscapes should increase the amount and stability of animal pollination-dependent crop production, and that an understanding of how the spatial pattern of land-use changefragmentation -impacts ecosystem services is key to achieve this goal. However, none of these studies have simultaneously combined crop dynamics at different scales and spatially-explicit landscapes to investigate the effects of land-use change on biodiversity and crop production services. Further, while multiple factors can influence crop production in fragmented agroecosystems, our current understanding of how they interact and determine crop production stability, and thus food security, remains incomplete.

There is general consensus that increased food production is necessary, but not sufficient, to achieve food security [START_REF] Godfray | Food security and sustainable intensification[END_REF], and that agriculture, especially in the global change context, should also aim at stabilizing crop production over time [START_REF] Schmidhuber | Global food security under climate change[END_REF][START_REF] Montoya | Reconciling biodiversity conservation, food production and farmers' demand in agricultural landscapes[END_REF]. Bearing this in mind, we here extend a model of crop dynamics into a spatially-explicit landscape to investigate how habitat loss and fragmentation, i.e. the amount and spatial configuration of seminatural habitat, influence the mean provision and stability of crop production in agricultural landscapes. We focus on animal-dependent crop production because (i) animal crop pollination is a key agricultural service that depends on biodiversity (pollinator animals), and (ii) worldwide agriculture is shifting towards more animal pollination-dependent food production systems [START_REF] Aizen | How much does agriculture depend on pollinators? Lessons from long-term trends in crop production[END_REF][START_REF] Breeze | Agricultural policies exacerbate honeybee pollination service supply-demand mismatches across Europe[END_REF]. Because the way food is produced worldwide threatens the existence of much of the world's biodiversity that contributes to crop pollination and food security, we explore how changes in biodiversity following land-use change affect animal-dependent crop production in fragmented agroecosystems. Thus, our model accounts for a variety of potential drivers of crop production mean and stability, such as loss and fragmentation of natural habitat and biodiversity, that are difficult to address collectively in empirical studies. Specifically, we address two questions: (i) How does the spatial pattern of land-use change influence the provision and stability of animal -dependent crop production in agroecosystems? (ii) How does biodiversity in fragmented landscapes influence animal-dependent crop production and food security?

Materials and Methods

To study the effects of habitat fragmentation on biodiversity and crop production, we extended a non-spatial model of crop dynamics [START_REF] Montoya | Trade-offs in provisioning and stability of multiple ecosystem services in agroecosystems[END_REF]) to a spatiallyexplicit landscape. Our model has two types of patches: crop land and seminatural habitat.

Crop land is used to grow annual crops with varying degrees of dependence on animal pollination, whereas seminatural habitat shelters biodiversity, including 'wild' plants and pollinators. The model does not consider managed honeybees as they do not depend on seminatural habitat for nesting, and generally pollinate less efficiently than non-managed pollinators [START_REF] Garibaldi | Wild pollinators enhance fruit set of crops regardless of honey bee abundance[END_REF]. Naturalized Apis species are implicitly considered in the model as they nest in seminatural habitat and forage across the landscape. This model represents intensively-managed agricultural systems, where crop land does not host significant levels of biodiversity, and spatial heterogeneity is broadly defined by two patch types.

The model investigates the expected biodiversity (i.e. species richness) and crop production, with a special focus on animal-dependent crop production, at the farm-(i.e. crop yield per area) and landscape-levels (i.e. the magnitude and stability of crop production). We distinguish between two additive ecosystem services associated with crop production: the production that results from wild animal pollination (hereafter animal-dependent production), and the production that is independent from animal pollination (animal-independent production), but can be wind-or self-pollinated.

Pollinators are assumed to be generalist central-place foragers that nest in seminatural habitat [START_REF] Gill | Protecting an ecosystem service: Approaches to understanding and mitigating threats to wild insect pollinators[END_REF]), yet move across the landscape to forage on either flowering crops or 'wild' plants, or both. Crop land and seminatural habitat are therefore linked by the pollinators' foraging movement. For simplicity, the model assumes similar foraging movements across pollinators. The main difference between the non-spatial and spatial models lies on the way pollinators link these two patch types: whereas in the non-spatial model pollinators move globally and reach all crop land, the spatial model imposes restrictions to pollinator's movement based on the spatial structure of the agricultural landscape (i.e. distribution of crop land and seminatural habitat) and the foraging range of pollinators (i.e. spillover). A conceptual representation of the model, including nonspatial and spatial components, as well as key model assumptions, is provided in Fig. 1.

In what follows, we first describe the spatial agroecosystem model and then present a simpler, mean-field approximation of the model that we use to produce the results.

Spatial agroecosystem model

The addition of space allows: (i) to explore a continuous gradient of land-use change patterns, from completely random to highly aggregated, that encompasses a wide range of fragmentation scenarios, (ii) to investigate the effects of fragmentation on pollinator diversity and crop production services, and (iii) to study how variations in the strength of spillover from seminatural habitat to crop land -i.e. the distance-decay of pollination flow -affects animal-dependent crop production. To do this, we simulated a continuous gradient of land-use change patterns, and therefore fragmentation, based on the aggregation degree of seminatural habitat. A complete description of the model parameters can be found in Table 1. Information on model parameterization is in Appendix S1.

Landscape pattern generation. The landscape consisted of a two-dimensional lattice (25×25 cells) where individual cells can have either of two states: crop land or seminatural habitat. We generated land-use change patterns by iteratively creating crop land cells in a landscape that consisted initially only of seminatural land. In a single step of the algorithm only one semi-natural habitat cell is selected and converted. The pattern generation is controlled by a parameter w that determines the clustering degree of the land-use change pattern (see Appendix S1). Hence, for w = 0 all seminatural land cells had the same relative probability to be chosen, leading to a fully random, unclustered pattern. Larger values of w resulted in more aggregated patterns. Therefore, variation in the value of w allowed us to produce a continuous gradient of land-use change patterns, and therefore fragmentation, based on the aggregation degree. For each land-use change pattern, we characterised fragmentation of the remaining seminatural habitat by quantifying mean fragment size, number of fragments, mean fragment perimeter, and perimeter: area ratio.

Fragmentation directly affects pollinator diversity (through the estimation of the metapopulation capacity) and the Landscape Pollination Potential (LPP), which in turn influence crop production services (see below).

General model dynamics. Pollinators, 'wild' plants and crop production are represented by Equations 1-3. In these equations, the first term defines growth and is modelled using the discrete Ricker equation, whereas the second and third terms add environmental and demographic stochasticity to the model: ! " ($ + 1) = ! " ($)exp (-.," ($) 01 -. 2 (3) " 4 5 6 7) + 8 . 9 : . 9 ($)! " ($) +

; 4 < = 4 < (3) >. 2 (3) ! " ($) (1)
? " ($ + 1) = ? " ($)exp (-@," ($) 01 -@ 2 (3) " A 5 6 7) + 8 @ 9 : @ 9 ($)? " ($) + is the total crop production at the end of the growing season t. rP,k (t), rW,k (t), rC,k (t) are the pollinators', 'wild' plants' and crop's per capita growth rates. ZC is the crop production independent of animal pollination, which allows varying degrees of animal-dependence production. Ck is thus the sum of animal-dependent and -independent parts of crop production. kP and kW are the carrying capacities of pollinators and 'wild' plants, respectively, per unit area. A1 is the area of a single cell; A is total landscape area; A[1wsn] is the total crop land area, and Awsn is total seminatural area. In the non-spatial model, k and A1 were not defined. Stochasticity is included in the form of σ u(t), where σ 2 is the variance ( e environmental or d demographic) and u(t) are independent Gaussian random variables with zero mean and unit variance, of either pollinators (σPuP), 'wild' plants (σWuW) or crops (σCuC) (see model stochasticity section below).

; A < = A < (3) >@ 2 (3) ? " ($) (2) B " ($) = C D EF G + -G," ($)H(1 + 8 G 9 : G 9 ($)) (3
Growth rate of pollinators and plants. Since pollinators are assumed to be generalist central-place foragers, pollinator's growth rate depends on the availability of resources ('wild' plants and crops) in the neighborhood. Plant and pollinator uptake of resources, as well as the animal-dependent crop production, follow a saturating, type II functional response, widely supported and consistent with real biological examples (Thebault & Fontaine 2010;[START_REF] Holland | Consumer-resource dynamics of indirect interactions in a mutualism-parasitism food web module[END_REF][START_REF] Holland | Population ecology of mutualism[END_REF]. Growth rates are thus defined by the following Monod/Michaelis-Menten equations:

-.," ($) = J 4 [ ∑ M 2N N∉P A N (Q)R∑ M 2N S N (Q) N∈P U ] W 4 X ∑ M 2N A N (Q) N∉P R∑ M 2N S N (Q) N∈P U (4) -@," ($) = J A ( ∑ M 2N 4 N (Q) N∉P U ) W A X ∑ M 2N 4 N (Q) N∉P U (5) -G," ($) = J S Y ∑ M 2N 4 N (Q) N∉P U Z W S X ∑ M 2N 4 N (Q) N∉P U (6)
, where εkl is the distance-decay function representing the decrease of ecosystem service flow from seminatural habitat to crop land (i.e. spillover effects, see below); L is the set of crop land cells; aP and aW are the maximum growth rates of pollinators and 'wild' plants, respectively; bP and bW are half-saturation constants of pollinators and 'wild' plants, respectively; aC is the maximum crop production derived from pollination, and bC is the half-saturation constant of crops.

Because crops vary in their dependency on animal pollination, different crop types can respond differently to landscape fragmentation. In our model, the degree to which crops depend on animal pollination is measured by αC/(ZC+αC), and reflects the part of total crop production dependent on animal pollination -i.e. animal-dependent production.

Landscape crop production is estimated by summing up the individual contribution of each cell k. Thus, for total crop production we obtained B($) = ∑ B " ($) "∈[ . We assume that crops are harvested yearly; hence, mean crop production represents the temporal mean of the yearly averaged crop production across the agricultural landscape. Crop yield per unit of agricultural area is calculated by dividing total crop production by crop land area.

Model stochasticity. To investigate the effects of fragmentation on crop production stability, the model includes environmental and demographic stochasticity (second and third terms, Eqs. 1-3). Stochasticity is included in the form of σ u(t), where σ 2 is the variance and u(t) are independent Gaussian random variables with zero mean and unit variance, of either pollinators (σPuP), 'wild' plants (σWuW) or crops (σCuC). Environmental stochasticity (σ e u e (t)) reflects variation in weather variables, e.g. temperature and rainfall, whereas demographic stochasticity (σ d u d (t)) emerges from stochastic variation in individuals' births and deaths, commonly observed in nature and important in small populations. Crops are sown at high densities, and thus we assume demographic stochasticity is prevented in crops. These two sources of stochasticity add biological realism to the model, as they allow variation in the biomass of pollinators, 'wild' plants and crops, useful to assess stability. We used 1/CV 2 as a metric of stability (invariability), where CV is the temporal coefficient of variation of crop production. This measure of stability is commonly used in ecological studies [START_REF] Loreau | Biodiversity and ecosystem stability: a synthesis of underlying mechanisms[END_REF].

Distance-decay of pollinators (spillover effects).

The spatially-explicit nature of the model allows to study how variations in the strength of spillover effects influences crop pollination. Ecosystem service distance-decay (εkl) affects the flow of pollination to crop land: the further crops are from seminatural habitat, the more difficult it is for pollinators to reach the crops, and thus the smaller the effect of pollinators on crop production. To model spillover effects, we used a logistic distance-decay function adapted from [START_REF] Mitchell | Strong and non-linear effects of fragmentation on ecosystem service provision at multiple scales[END_REF] (Appendix S1). This function is consistent with both theoretical predictions [START_REF] Ries | Ecological responses to habitat edges: mechanisms, models, and variability explained[END_REF][START_REF] Rand | Spillover edge effects: the dispersal of agriculturally subsidized insect natural enemies into adjacent natural habitats[END_REF][START_REF] Ricketts | Landscape effects on crop pollination services: are there general patterns[END_REF]) and empirical observations [START_REF] Farwig | Isolation from forest reduces pollination, seed predation and insect scavenging in Swiss farmland[END_REF][START_REF] Mitchell | Forest fragments modulate the provision of multiple ecosystem services[END_REF][START_REF] Martins | Pollination services are mediated by bee functional diversity and landscape context[END_REF] of the effects of habitat edges and distance-to-habitat fragment gradients on ecosystem service provision. The logistic distance-decay function has two parameters: (i) the distance from seminatural cell edge at which the pollination flow equals one half of its initial value -50% decay distance -(dm), and (ii) the distance over which the pollination flow decreases from 90% to 10% of its initial value -90% to 10% decay distance (∆d).

Fragmentation, biodiversity and crop pollination. To consider the effects of fragmentation on biodiversity and crop production, we included two relationships in our model: (i) the dependence of pollinator diversity on the amount and distribution of seminatural habitat, and (ii) the dependence of crop production on pollinator diversity. In the first case, we used the Species-Fragmented Area Relationship (SFAR; [START_REF] Hanski | Species-fragmented area relationship[END_REF]), which extends one of the oldest known and most documented patterns in ecology -the species-area relationship (SAR) -to fragmented landscapes. Despite recent debate has ensued on the relative importance of habitat loss versus fragmentation on species diversity [START_REF] Rybicki | Habitat fragmentation and species diversity in competitive communities[END_REF][START_REF] Fletcher | Is habitat fragmentation good for biodiversity?[END_REF][START_REF] Fahrig | Is habitat fragmentation bad for biodiversity?[END_REF], empirical evidence shows that larger and more connected fragments of natural habitat in general host more biodiversity than smaller, more isolated fragments [START_REF] Haddad | Habitat fragmentation and its lasting impact on Earth's ecosystems[END_REF]. In agricultural landscapes, this means that different fragmentation patterns will affect biodiversity and crop production in different ways. In an agricultural landscape, the SFAR can be modeled as a power law: \ = ] ^(C_ ^`) a b , where S is pollinator diversity, Awsn is the total area of seminatural habitat, and zs is the power-law exponent. The prefactor ] c is not a constant as in the standard SAR; rather, it decreases in more fragmented landscapes (see Appendix S1 for more details).

For the second relationship -the dependency between crop production and pollinator diversity -, we created a dependence of pollinator's carrying capacity (kP) on biodiversity following a power law: d . = ] . \ a 4 , where S is the number of pollinator species estimated by the SFAR, and ] . , e . are parameters of the power law [START_REF] Liang | Positive biodiversity-productivity relationship predominant in global forests[END_REF]O'Connor et al. 2017). Although pollinators can differ in their pollination efficiency [START_REF] Matsuki | Pollination efficiencies of flowervisiting insects as determined by direct genetic analysis of pollen origin[END_REF][START_REF] Kleijn | Delivery of crop pollination services is an insufficient argument for wild pollinator conservation[END_REF][START_REF] Willcox | Deconstructing pollinator community effectiveness[END_REF]), such differences do not dilute the positive relationship between pollinator diversity and crop production and yield, as reported in meta-analytical studies [START_REF] Woodcock | Meta-analysis reveals that pollinator functional diversity and abundance enhance crop pollination and yield[END_REF]. Additionally, we considered the ability of pollinator diversity to provide an insurance against environmental fluctuations.

Following the biodiversity insurance hypothesis [START_REF] Loreau | Biodiversity as spatial insurance in heterogeneous landscapes[END_REF], pollinator diversity insures crop production because many pollinators better maintain crop production if pollinators differ in their responses to environmental variation. We followed the ecological literature and made pollinators' environmental stochasticity inversely related to their diversity (8 . 9 = f . \ g h ; Tilman 1999). Demographic stochasticity acts at the individual level, and in the same manner for conspecifics and heterospecifics; thus, there is no insurance effect for demographic stochasticity.

Non-spatial model results: summary

The relationship between mean crop production and the proportion of seminatural habitat is hump-shaped, and the height and position of the hump depends on crop pollination dependence (αC/(ZC+αC)) and the crop relative requirement for pollinator densities (βC/kP) (Fig. S1; [START_REF] Montoya | Trade-offs in provisioning and stability of multiple ecosystem services in agroecosystems[END_REF]). Higher and lower values of these factors, respectively, shift maximum production to higher fractions of seminatural habitat. On the other hand, the stability of animal-dependent crop production is generally an increasing function of proportion of seminatural habitat and is also determined by the same factors: higher crop pollination dependence and crop relative requirement for pollinators decreases stability.

Mean-field approximation of the spatial model

Given the computational demands of the model, we developed a mean-field approximation which replaces the detailed spatial flows between seminatural habitat and crop land cells by their values spatially-averaged over the landscape. More precisely, we substituted the pollination flow decay coefficients (Eqs. 4-6) by the average i (i.e. the average value of i "j when taking a random cell k ∈ L and a random cell l ∉ L):

i = D (Dkl bm )`n D l bm `n ∑ i "j "∈[,j∉[ (7)
This quantity has two complementary interpretations (see Appendix S2 for the derivation). First, it is the fraction of seminatural habitat from which a crop land cell can be pollinated, averaged over all crop land cells. Second, it is the fraction of crop land that is reachable by pollinators from a seminatural habitat cell, averaged over all seminatural cells. Taken together, these two interpretations can be summarized by the term Landscape Pollination Potential or LPP (LPP replaces i hereafter; 0≤LPP≤1).

The mean-field approximation is a very accurate description of the ecosystem service dynamics in agricultural landscapes (Fig. 2; Appendix S2), and shows that the spatial effects of fragmentation on animal-dependent crop production are determined by LPP.

Thus, LPP is a metric that captures fragmentation and pollinator spillover effects within agricultural landscapes. To consider the spatial structure of land-use change, the term βC/kP of the non-spatial model must be replaced by

D [.. W S " 4 (8)
where βC/kP is the crop relative requirement for pollinators (described above).

Results

The results here presented correspond to the mean-field approximation, as it provides a very accurate description of the model dynamics (Fig. 2). We present three groups of results: (i-ii) the effects of LPP and biodiversity on crop production, and (iii) the net effects of fragmentation on animal-dependent production.

Effects of LPP on crop production

When LPP = 1, fragmentation effects are negligible and crop dynamics are identical to those of the non-spatial model (Fig. 3A-C, dark-blue lines). In this case, we retrieve the same mechanisms of the non-spatial model [START_REF] Montoya | Trade-offs in provisioning and stability of multiple ecosystem services in agroecosystems[END_REF]. The effects of fragmentation kick off when LPP < 1. Lower LPP reduces the carrying capacity of pollinators effectively; that is, the effects of LPP on the mean and stability of crop production can be fully taken into account by changing the value of the pollinators' carrying capacity from kP to LPP*kP (Eq. 8). LPP-driven reduction in pollinators' carrying capacity decreases the mean provision of animal-dependent production (Fig. 3A)

and crop yield per area (crop yield per area increases with seminatural habitat, but larger LPP saturates this relationship earlier; Fig. 3C). The same is true for the stability of animal-dependent production, except at small fractions of seminatural habitat and/or small values of LPP (Fig. 3B). LPP has no effect on animal-independent crop production as it does not depend on animal pollination and, therefore, on seminatural habitat; its mean value decreases linearly with seminatural habitat, whereas its stability is solely determined by environmental stochasticity (Fig. 3D-E). Because of this, animalindependent crop production is not further considered in subsequent results. Whereas the full complexity of the spatial fragmentation effects (i.e. those not mediated by biodiversity) on animal-dependent production are captured by LPP (Fig. S2), we did not find any clear, consistent effect of specific fragmentation metrics (Fig. S3), except for 'Number of crop land cells > distance threshold'. This is because the latter variable uses distance-decay parameter dm as the reference threshold, and it is thus correlated with LPP.

However, the dependence of both metrics on dm does not really explain their correlation:

LPP measures the fraction of crop land within a distance dm of seminatural habitat and quantifies the interspersion of seminatural habitat and crop land at the scale of the pollination spillover distance dm, thus capturing an important dimension of fragmentation (i.e. interspersion) that is not or poorly considered by the standard fragmentation metrics.

Finally, the effect of LPP on animal-dependent production increases with the degree to which crops depend on animal pollination: higher pollination dependence shifts maximum production to higher fractions of seminatural habitat at landscape and local scales, and the stability of animal-dependent production increases faster (Fig. S4), consistent with [START_REF] Montoya | Trade-offs in provisioning and stability of multiple ecosystem services in agroecosystems[END_REF].

LPP depends on two factors: fragmentation and spillover, or the distance-decay of ecosystem service flow. High aggregation (low fragmentation) and fast distance-decay result in lower LPP (LPP<1), which in turn reduces animal-dependent crop production.

These two factors interact: only when the flow of pollinators to crop land is limited (fast distance-decay) aggregation patterns influence animal-dependent production. When no restrictions exist in the flow of pollinators to crop land, LPP is maximum and fragmentation does not affect animal-dependent production (LPP ≈ 1).

Effects of biodiversity on crop production

In parallel to LPP, pollinator diversity also impacts animal-dependent crop production.

As expected, biodiversity enhances animal-dependent production at local (i.e. yield per area) and landscape scales (Fig. 4A,C), and stabilizes animal-dependent production (Fig. 4B) by increasing the pollinators' carrying capacity, and by reducing the response of animal-dependent production to environmental fluctuations. A higher biodiversity effect (larger e . ) increases both mean animal-dependent production and its stability, as well as yield per area.

Net effects of fragmentation on crop production

Fragmentation has a dual effect on animal-dependent crop production. On one hand, aggregation of seminatural fragments decreases pollination by lowering LPP, which in turn reduces the pollinators' carrying capacity. On the other hand, aggregation benefits biodiversity (Fig. S5), which in turn increases pollinators' biomass and animal pollination. The net effect of fragmentation on animal-dependent production depends on the distance-decay of ecosystem service flow (dm) and the proportion of seminatural habitat remaining. When the decay distance is low (Fig. 5, first row), fragmentation effects tend to be positive for mean animal-dependent production and yield per area because the fraction of crop land within reach from non-crop land areas is higher. Yet, the stability of animal-dependent production decreases due to the lower biodiversity in fragmented landscapes, except at high fractions of seminatural habitat where the impact of fragmentation is minimum. Conversely, when the decay distance is high, seminatural fragments are perceived as more connected and animal-dependent production is not limited by space. In this case, fragmentation becomes irrelevant, or even negative, due to the lower biodiversity in fragmented landscapes (Fig. 5, last row).

Discussion

Our The mean-field approximation shows that the effects of space on animal-dependent production can be interpreted in the same terms as varying the pollinator's carrying capacity in the non-spatial model, i.e. a lower LPP reduces the pollinator's carrying capacity and the provision of animal-dependent production.

Our results suggest that understanding the factors that affect Landscape Pollinator Potential is a fundamental step towards food security in animal-dependent agriculture. If no restrictions exist in the flow of pollinators to crop land, LPP is maximum and the spatial structure of land-use change does not affect crop dynamics. In this situation, seminatural fragments are perceived as more connected and the provision and stability of animal-dependent production is not conditioned by the spatial configuration, i.e. spatial and non-spatial models converge. However, agricultural landscapes are fragmented to some extent and the foraging ranges of most organisms are local [START_REF] Darvill | Use of genetic markers to quantify bumblebee foraging range and nest density[END_REF][START_REF] Zurbuchen | Maximum foraging ranges in solitary bees: only few individuals have the capability to cover long foraging distances[END_REF][START_REF] Geib | Bumble bee nest abundance, foraging distance, and host-plant reproduction: implications for management and conservation[END_REF], which produces higher aggregation and weaker spillover effects, thus reducing LPP. Such reductions in LPP affect crop production by (i) decreasing mean animal-dependent production and total yield per area, and (ii) decreasing the stability of animal-dependent production along the seminatural habitat gradient.

Because LPP can have important implications for crop production, it would be interesting to operationalize this metric in empirical studies and landscape management. The estimation of LPP in real farming systems would require data on the aggregation of seminatural habitat, and on the spillover of pollinators to adjacent crops. The former can be obtained with GIS processing of aerial pictures or satellite images. For the latter, information on foraging distances of existing pollinator species (this comes from ecological censuses) combined with experimental studies could be used to reveal species' foraging patterns and how the flow of pollinators to adjacent crop land decays with distance (examples of this type of studies are [START_REF] Ries | Ecological responses to habitat edges: mechanisms, models, and variability explained[END_REF][START_REF] Rand | Spillover edge effects: the dispersal of agriculturally subsidized insect natural enemies into adjacent natural habitats[END_REF][START_REF] Ricketts | Landscape effects on crop pollination services: are there general patterns[END_REF]. This information will be useful to design agricultural landscapes for high LPP. Although measuring LPP in empirical studies demands more data than simple, traditional fragmentation metrics, this investment is worthy as LPP can be a better predictor of crop production.

Producing food requires land, and increasing the land devoted to farming reduces the land devoted to biodiversity conservation. This and other recent empirical studies show that pollinator diversity can increase food production [START_REF] Catarino | Bee pollination outperforms pesticides for oilseed crop production and profitability[END_REF][START_REF] Dainese | A global synthesis reveals biodiversity-mediated benefits for crop production[END_REF][START_REF] Woodcock | Meta-analysis reveals that pollinator functional diversity and abundance enhance crop pollination and yield[END_REF], and lead to lower variability in agricultural productivity. This is relevant in the global change context, as biodiversity can stabilize animaldependent crop production by providing an insurance against environmental fluctuations, which are predicted to intensify under global change [START_REF] Giorgi | Regional climate information-evaluation and projections[END_REF][START_REF] Saltz | The impact of increased environmental stochasticity due to climate change on the dynamics of Asiatic wild ass[END_REF][START_REF] Fischer | Robust spatially aggregated projections of climate extremes[END_REF]. The response of biodiversity to land-use change depends on the amount and the spatial structure of seminatural habitat. For example, although the effects of fragmentation on biodiversity are stronger at low-intermediate fractions of seminatural habitat -typical of intensive farming systems -, aggregation favours biodiversity within seminatural habitat fragments. The stabilizing effect of biodiversity and its role in food security is increasingly supported, even at crop levels [START_REF] Renard | National food production stabilized by crop diversity[END_REF]. Our results add to this view and point to biodiversity conservation as one key policy to achieve food security.

Our findings are consistent with previous studies that found non-linear effects of fragmentation on ecosystem services [START_REF] Keitt | Habitat conversion, extinction thresholds and pollination services in agroecosystems[END_REF][START_REF] Mitchell | Strong and non-linear effects of fragmentation on ecosystem service provision at multiple scales[END_REF], and provide a theoretical basis of the effects of fragmentation on the stability of animal-dependent crop production. Fragmentation has a dual effect on animal-dependent production. On one hand, aggregation decreases animal-dependent production by reducing the Landscape Pollination Potential. On the other hand, aggregation increases animal-dependent production by maintaining higher biodiversity, especially when the proportion of seminatural habitat within the agricultural landscape is medium or low. The net effects of aggregation on animal-dependent crop production depend on the strength of spillover effects. These results have management implications (e.g. land sharing-sparing debate; [START_REF] Fischer | Land sparing versus land sharing: moving forward[END_REF][START_REF] Grass | Land-sharing/-sparing connectivity landscapes for ecosystem services and biodiversity conservation[END_REF], as the goals of different landscape managers can be conditioned by the way that natural land is converted into crops. For example, maintaining a large number of seminatural fragments may be a better strategy at multiple spatial scales than maintaining a few large fragments when pollinator flow to crop land is low. Yet, this strategy may increase the temporal variability of animal-dependent production at low-intermediate proportions of seminatural habitat, reflecting a trade-off between ecosystem service mean and stability. Conversely, larger fragments of seminatural habitat have higher pollinator diversity when the fraction of seminatural habitat is low or intermediate, and higher biodiversity can stabilize animal-dependent production. These results agree with recent claims that the land sharing-sparing dichotomy lends itself to overly simplistic policy prescriptions [START_REF] Kremen | Reframing the land-sparing/land-sharing debate for biodiversity conservation[END_REF], and suggest that management decisions for food security should consider factors such as pollinators' spillover, the amount and spatial aggregation of seminatural habitat and the animal pollination-dependence of crops. Our results also provide recommendations for landscapes with both animal-independent and -dependent crops. In such cases, for a given proportion of seminatural habitat, management should maximize the Landscape Pollination Potential while preserving pollination diversity, favouring mean crop production and its stability. This could be achieved by increasing fragmentation of seminatural habitat in the vicinity of animal-dependent crops, but without making semi-natural habitat patches too small so that biodiversity is not negatively affected. Animalindependent production only depends on how much seminatural habitat there is, and is not affected by fragmentation because it does not depend on animal pollinators and seminatural habitat. as well as the effect of pollinator diversity on animal-dependent crop production [START_REF] Liang | Positive biodiversity-productivity relationship predominant in global forests[END_REF]O'Connor et al. 2017).

Our model has several limitations. First, our model focuses on intensive farming systems, where crop land does not host important biodiversity levels; other types of agriculture -e.g. organic farming -allow moderate levels biodiversity to thrive within crop land, and can modify the results reported here [START_REF] Clough | Combining high biodiversity with high yields in tropical agroforests[END_REF]). The model does not consider fine-scale seminatural patches (e.g. wildflower strips, hedgerows) that might support high pollinator diversity by providing complementary resources within a short range from crops, nor it considers habitat quality differences within seminatural fragments; although the latter can affect the abundance and composition of pollinators, it would not qualitatively affect the positive effect of pollinator diversity on animal-dependent production. Also, the observation that biodiversity loss has either none (stability) or positive (mean) effects on animal-independent crop production may change if organisms responsible for other services, e.g. pest control, are included. Besides, although seminatural habitat had no effect on the stability of animal-independent production, this may change if environmental stochasticity of crops increases with decreasing amounts of seminatural habitat, as suggested by studies linking seminatural habitat to climate regulation, natural hazard regulation and water flow regulation services [START_REF] Harrison | Identifying and prioritizing services in European terrestrial and freshwater ecosystems[END_REF]). In addition, negative perturbations larger than the ones considered in our model may be possible; although these so-called black-swan events are rare (e.g. they affect ~3% of insect populations; [START_REF] Anderson | Black-swan events in animal populations[END_REF], they could heavily decrease pollinator diversity, with negative impacts on animal-dependent crop production. Differences in pollination efficiency may affect LPP and could be captured using functional group analysis (Mayfield et al. 2019); yet, the positive diversity-function relationship still holds even considering such differences [START_REF] Woodcock | Meta-analysis reveals that pollinator functional diversity and abundance enhance crop pollination and yield[END_REF]. Also, differences in flight ability and range across pollinators influence pollinator diversity within agricultural landscapes (e.g. InVEST Crop pollination model, [START_REF] Lonsdorf | Modelling pollination services across agricultural landscapes[END_REF], and this may in turn affect LPP, through its effects on pollinator spillover, and the provision of animaldependent production. Ensuring stable food supplies is a challenge that may require multiple solutions. Policies to increase production, changing diets, irrigation, crop diversity, and tolerance of crops to drought, are proposed as stability-enhancing solutions [START_REF] Rosenzweig | Potential impact of climate change on world food supply[END_REF][START_REF] Lobell | Prioritizing climate change adaptation needs for food security in 2030[END_REF][START_REF] Godfray | Food security and sustainable intensification[END_REF][START_REF] Bailey | Extreme Weather and Resilience of the Global Food System. Final Project Report from the UK-US Task force on Extreme Weather and Global Food System Resilience[END_REF][START_REF] Renard | National food production stabilized by crop diversity[END_REF].

By addressing multiple drivers of crop production in spatially-explicit agroecosystems, our study provides a theoretical basis of the effects of fragmentation on the mean and stability of animal-dependent crop production, with strong consequences for food production and food security. These results are highly relevant in the global change context, and given the worldwide trends in agriculture, which shifts towards more pollinator-dependent crops. Table 1 |

Figures and Tables

  analysis reveals a variety of effects of land-use change on biodiversity and crop production. Using a mean-field approximation, our model suggests that (i) the full complexity of the fragmentation-induced spatial effects on animal-dependent crop production is captured by one factor -the Landscape Pollination Potential of the remaining seminatural land (LPP) -which determines the mean and stability of animaldependent production; (ii) fragmentation reduces biodiversity and increases LPP, thus impacting animal-dependent production in opposite directions; and (iii) the net effects of fragmentation on animal-dependent production depend on the strength of pollinators' spillover to crop land, the proportion of seminatural habitat remaining, and the degree to which crops depend on animal pollination. The loss of seminatural land has contrasting effects on agricultural landscapes: biodiversity decreases, animal-independent crop production increases, while animaldependent crop production is maximized at intermediate fractions of seminatural habitat. But fragmentation can modify these relationships in two ways. First, land-use change can produce multiple patterns of aggregation of seminatural habitat fragments. These patterns, combined with the strength of the spillover of pollinators to crop land, determine the Landscape Pollination Potential of seminatural land, which is one main driver of food production in animal-dependent agriculture. The second type of effects are mediated by biodiversity, as the level of aggregation of seminatural habitat affects the pollinator richness. Such purely spatial and biodiversity-mediated effects modify the carrying capacities of pollinators, which ultimately determine animal-dependent crop production.

  Parameters and variables of the model. An asterisk in the unit column indicates that the value of that parameter depends on fragmentation.

Figure 1 |

 1 Figure 1 | Conceptual diagram of our modelling framework. Green boxes represent

Figure 2 |

 2 Figure 2 | Mean-field approximation vs Exact solution. (a) Mean animal-dependent

Figure 3 |

 3 Figure 3 | Effects of landscape composition and LPP on crop production. (A) Mean

  )

	Pk represents pollinators, Wk 'wild' plants and Ck crop production in cell k, with ! " ($) =
	? " ($) = 0 if k is a crop land cell, and B " ($) = 0 if k is seminatural habitat. One unit of
	time t corresponds to one growing season, P(t) and W(t) can be interpreted as total
	biomass of pollinators and 'wild' plants over growth season t, respectively, whereas C(t)

  Aside from Landscape Pollination Potential (see above), what type of empirical data could be used to calibrate key parameters of the model? Carrying capacities of pollinators and wild plants can be approximated from species biomass and ecological censuses of species numbers (Appendix S1), yet this assumes that communities are saturated. On the

	other hand, animal pollination-dependencies of the large majority of crops is available
	from pollinator exclusion experiments (e.g. Klein et al. 2007). Pollinator spillover
	depends on the distance-decay function, and its parameters values require data on
	pollinator movement patterns (distance, range), which can be obtained via experiments
	and observational studies (e.g. tracking of 'marked' pollinator individuals with remote
	devices; e.g. Ratnayake et al. 2021). Finally, experimental studies provide information to
	calibrate biodiversity estimates in fragmented landscapes (SFAR function; Appendix S1),

  Therefore, differences in movement/flight patterns of pollinators should be considered in future extensions of the model, e.g. large bees such as Bombus species are powerful long-distance pollinators, whereas smaller bees forage locally (e.g. Lasioglossum, Halictus). Finally, our model focuses on wild central-place pollinators (i.e. all types of wild bees, including bumble bees and solitary bees), whose presence and abundance directly depend on seminatural habitat. Non-bee pollinators are excluded as they have diverse nesting habits (e.g. many flies lack central nest locations, and others depend on floral resources only during adult life). However, non-bee pollinators respond less negatively to land-use changes[START_REF] Rader | Non-bee insects are important contributors to global crop pollination[END_REF], and their role in stabilizing animal-dependent production is important in the event of bee declines. Honey bee colonies are used to substitute wild pollinators, yet, with a few exceptions (e.g. honeybees complement and sometimes efficiently pollinate crops of the family Cucurbitaceae o

Solanacea;

[START_REF] Macias-Macias | Contribution of native bees and Africanized honey bees (Hymenoptera:Apoidea) to Solanaceae crop pollination in tropical México[END_REF]

, the pollination services of wild pollinators cannot be compensated by managed bees because (1) pollinator-dependent crop land grows more rapidly than the stock of honey bee colonies

[START_REF] Lindström | Experimental evidence that honeybees depress wild insect densities in a flowering crop[END_REF]

), (2) wild insects generally pollinate crops more efficiently than honeybees

(Garibaldi et al. 2013), and (3)

honeybees may depress wild pollinator densities

[START_REF] Winfree | Native bees provide insurance against ongoing honey bee losses[END_REF]

). Wild central-place foragers thus remain a fundamental group of crop pollinators in agriculturally-dominated landscapes

[START_REF] Potts | Safeguarding pollinators and their values to human well-being[END_REF]

.

Acknowledgments DM was funded by the European Union and Institut National de Recherche pour l'Agriculture, l'Alimentation et l'Environnement (INRAE) in the framework of the Marie-Curie FP7 COFUND People Program, through the award of an AgreenSkills/AgreenSkills+ fellowship, and by and the FRAGCLIM Consolidator Grant, funded by the European Research Council under the European Union's Horizon 2020 research and innovation programme (grant agreement number 726176). This work was supported by the TULIP Laboratory of Excellence (ANR-10-LABX-41) and by the BIOSTASES Advanced Grant funded by the European Research Council under the European Union's Horizon 2020 research and innovation program (grant agreement no.

666971). We also thank three anonymous referees for their insightful comments.

Data availability statement

Data available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.8sf7m0cn8 [START_REF] Montoya | Habitat fragmentation and food security in crop pollination systems[END_REF].

). On the other hand, biodiversity impacts the stability of animal-dependent crop production both indirectly -increasing the carrying capacity of pollinators -and directly -reducing the response of animal-dependent production to environmental fluctuations (8 . 9 = f . \ g ⁄ ).

For each ecosystem service, the plots compare two scenarios: (i) a scenario where biodiversity has no effect on animal-dependent production (e . = 0, q = 0), represented by the grey shadows, versus (ii) a scenario where biodiversity has an effect on animaldependent production (e . = 0.26, q = ½) [START_REF] Tilman | The ecological consequences of changes in biodiversity: a search for general principles[END_REF][START_REF] Liang | Positive biodiversity-productivity relationship predominant in global forests[END_REF]