Laetitia Dubois
email: laetitia.dubois.2008@polytechnique.org

Jérôme Mathieu
email: jerome.mathieu@upmc.fr

Nicolas Loeuille
email: nicolas.loeuille@upmc.fr

The manager dilemma: optimal management of an ecosystem service in heterogeneous exploited landscapes

Ecosystem services are defined as benefits produced by ecological communities, supporting human welfare. Because sustainable agriculture relies on such ecosystem services, finding the optimal management -which optimizes both the surface dedicated to human activities and the delivery of ecosystem services -is particularly critical. Ecosystem services heavily depend on the presence and activity of organisms, especially ecosystem engineers. In order to find the proportion and the spatial aggregation of exploited areas that optimize an ecosystem service, we developed three complementary metapopulation models of a keystone species in an exploited landscape. We considered both anthropic and ecological constraints, by modelling the simultaneous management of two variables: the yield of human activities and the ecosystem service provided by the metapopulation. We also investigate how this optimal management can drive the metapopulation close to extinction, and how two key ecological traits of species -population growth and dispersal rates -can mitigate such extinction risks. The two spatially implicit metapopulation models show that the optimal management is a trade-off, benefits often being optimized for intermediate surfaces of exploitation. This optimal surface depends on the ecological traits and on the degree of disturbance incurred by human activities. Spatially explicit simulations suggest that optimal management is further improved when the spatial distribution of human activities is fragmented.

Keywords

Ecosystem services; Landscape; Heterogeneity; Metapopulation; Ecosystem engineer; Source-sink dynamics.

Ecosystem services refer to the benefits human populations derive, directly or indirectly, from ecosystem functions [START_REF] Costanza | The value of the world's ecosystem services and natural capital[END_REF]. They are most usually provided by ecological communities embedded in a given landscape and have positive effects on human welfare.

A landscape is a mosaic of habitat patches [START_REF] Dunning | Ecological processes that affect populations in complex landscapes[END_REF]. The landscape can be either marine or terrestrial. The patch is the basic building block. A metapopulation (see below) can disperse within the landscape.

A patch is here assumed to be either exploited by human beings, either non-exploited. A patch can be occupied or not by the population of interest, which provides the ecosystem service.

Heterogeneity refers to a spatial characteristic of the landscape [START_REF] Turner | Landscape Ecology: The Effect of Pattern on Process[END_REF]. Here, the landscape is an environmental mosaic, in which heterogeneity has two components: the amplitude of the environmental alteration due to human activitiescompositional heterogeneity [START_REF] Fahrig | Functional landscape heterogeneity and animal biodiversity in agricultural landscapes[END_REF], and the spatial aggregation of the anthropized patches -configurational heterogeneity [START_REF] Fahrig | Functional landscape heterogeneity and animal biodiversity in agricultural landscapes[END_REF].

Fragmentation (as the reverse of spatial aggregation) is one component of spatial heterogeneity. The more fragmented the landscape, the less aggregated the exploited -or non-exploited -patches.

A metapopulation is a population of populations, connected by dispersers [START_REF] Levins | Some demographic and genetic consequences of environmental heterogeneity for biological control[END_REF]. Here, we assume that the metapopulation is made of local populations that live either in exploited patches or in non-exploited patches.

An ecosystem engineer species directly or indirectly modulates the availability of resources to other species by causing physical or chemical changes in biotic or abiotic materials [START_REF] Jones | Organisms as ecosystem engineers[END_REF].

A species niche is defined from its relationships with the biotic and abiotic elements of its environment. In other words, each species has a niche which is the intersection of all of the ranges of tolerance under which it can live [START_REF] Hutchinson | Concluding remarks[END_REF].

A population with a positive growth rate and which loses more emigrants than it receives immigrants is called a sourcepopulation. Conversely, a population with a negative growth rate, with a greater immigration than emigration is called a true sink population [START_REF] Pulliam | Sources, sinks and population regulation[END_REF]. It differs from a pseudo-sink population that has positive growth rates, but still receives more immigrants than produces emigrants due to spatial heterogeneities in environmental conditions (Watkinson and Sutherland, 1995). Here, non-exploited patches are sources and exploited patches can be either true sinks or pseudo-sinks.

Introduction

The concept of ecosystem services arises from the need to define the relationships between human welfare and ecosystems. Ecosystem services can be defined as the goods and services that human populations derive from ecosystems [START_REF] Costanza | The value of the world's ecosystem services and natural capital[END_REF]. Several types of ecosystem services can be distinguished: provisioning services, such as water or wood; regulating services -such as water purification; supporting services -such as biological cycles; and cultural services -like ecotourism and aesthetic value. Most of them heavily depend on organisms' presence or activity. Some species delivering ecosystem services may be unaffected or even be favored by human activities within an agricultural landscape [START_REF] Eriksson | Species pools in cultural landscapes -niche construction, ecological opportunity and niche shifts[END_REF]. In such instances, economical and ecological benefits are positively correlated, and the optimal solution is to devote the whole landscape to the considered activity. However, in most situations, human exploitation of ecosystems is harmful to species delivering ecosystem services, either reducing their abundance or their biodiversity. The 2005 Millennium Ecosystem Assessment pointed out the irreversible changes humans have caused over the last fifty years on ecosystems, resulting in a threat of extinction on 10 to 30 % of mammal, bird, and amphibian species. For instance, in intensively managed agricultural landscapes, many pollinator populations have decreased, causing large decreases in the provisioning of the pollination ecosystem service [START_REF] Biesmeijer | Parallel Declines in Pollinators and Insect-Pollinated Plants in Britain and the Netherlands[END_REF][START_REF] Potts | Global Pollinator Declines: Trends, Impacts and Drivers[END_REF]. When human activities harm species providing ecosystem services, a dilemma emerges: which surface should be exploited, considering both human welfare and species conservation? How should we organize human activity to keep it sustainable?

While landscapes are increasingly disturbed by human activity, conservation of engineers appears all the more urgent.

Several constraints -biotic and abiotic -affect the outcome of this dilemma. On the abiotic side, environmental spatial heterogeneity is a key factor. Two components play an important role. First, the amplitude of the environmental alteration due to human activities, expressed as the degradation of the species ecological niche (i.e., the reduction in the species growth rate). Second, the spatial aggregation of the human perturbation on the environment, affecting the ecological dynamics at landscape scales.

Based on this idea, some theoretical bioeconomics models take into account as many abiotic constraints as possible in order to optimize the human exploitation of a given resource [START_REF] Sanchirico | Bioeconomics of Spatial Exploitation in a Patchy Environment[END_REF].

Many empirical studies illustrate the critical role of spatial heterogeneities in the management of ecosystem services. Fragmentation has for instance been shown to increase extinction risks in disturbed landscapes [START_REF] Fahrig | Conservation of Fragmented Populations[END_REF]. Landscape structure also strongly affects the spatial distribution of wild populations. For instance, the abundance and structure of vole populations widely change depending on landscape type -from a village to an open field [START_REF] Delattre | Effect of landscape structure on common vole (microtusarvalis) distribution and abundance at several space scales[END_REF]. Landscape heterogeneity management is therefore at the heart of agricultural practice issues, as exemplified by the land sharing vs. land sparing debate [START_REF] Green | Farming and the Fate of Wild Nature[END_REF]. This debate tackles the future food security.

It discusses how management option trades-off exploitation and species conservation in space, by opposing farming the entire surface with wildlife friendly techniques -i.e. land sharing -vs. farming intensively some land whilst other land remains as a nature reserve -i.e. land sparing. Such options indirectly question the issue of landscape heterogeneity: is it optimal to reduce such heterogeneities (by exploiting a large group of homogeneous nature-friendly patches) or to strongly exploit fewer patches? Biological aspects inherent to the ecosystem engineer species are equally important. Variation in dispersal capacities or reproduction rates modulates the consequences of human modifications for the persistence and functioning of such populations. Metapopulation [START_REF] Hanski | Single-species metapopulation dynamics: concepts, models and observations[END_REF] and metacommunity [START_REF] Leibold | The metacommunity concept: a framework for multi-scale community ecology[END_REF]) models incorporate such aspects, and offer good opportunities to determine suitable strategies for a sustainable management of ecosystem services accounting for both conservation and economic issues. The role of non-exploited patches is highlighted when it comes to agricultural landscapes [START_REF] Burel | Habitat quality and connectivity in agricultural landscapes: The role of land use systems at various scales in time[END_REF] because the spatial variations in habitat quality may increase extinction risks for many species [START_REF] Hanski | Single-species metapopulation dynamics: concepts, models and observations[END_REF] but also provide adjacent ecosystems with services due to population spillovers (Tylianakis et al., 2007;[START_REF] Loeuille | Eco-Evolutionary Dynamics of Agricultural Networks: Implications for Sustainable Management[END_REF]. Dispersal, combined with environmental autocorrelation, define the environmental "grain" under which the species demography, evolution and management should be considered [START_REF] Macarthur | The limiting similarity, convergence, and divergence of coexisting species[END_REF]. Adapting spatial models also requires considering how the ecosystem service is delivered by populations. Ecosystem services can be due to the presence of the species, independently of its density (presence-dependent ecosystem service: e.g., [START_REF] Byers | Using ecosystem engineers to restore ecological systems[END_REF] or proportional to the species density (density-dependent ecosystem service: e.g., [START_REF] Hodgson | Comparing organic farming and land sparing: optimizing yield and butterfly populations at a landscape scale[END_REF].

Several models tackle the issue of the optimal exploited proportion of a landscape. This is for instance a classical topic concerning the design of marine reserves, to improve the fishing efficiency and preserve wild populations [START_REF] Baskett | Designing Marine Reserves for Interacting Species: Insights from Theory[END_REF][START_REF] Gaines | Marine Reserves Special Feature: Designing Marine Reserve Networks for Both Conservation and Fisheries Management[END_REF], White et al., 2008). Pollination services can be similarly modeled, integrating ecological and economic constraints to define the optimal design of a landscape [START_REF] Brosi | Optimal design of agricultural landscapes for pollination services[END_REF]. Such previous studies have often focused on specific situations with a detailed description of a particular ecosystem service with its inner constraints. Here, we intend to develop models that are more general in order to focus on the two constraints we view as key in determining ecosystem service management: species dispersal and spatial environmental heterogeneity. Particularly, we aim at understanding how these two constraints interact.

We consider a species in a heterogeneous landscape made of two different types of patches: exploited and non-exploited. Individuals disperse throughout the landscape, and provide an ecosystem service which positively affects human activities. The exploitation activity has a negative impact on the species survival. We take the point of view of a landscape manager who wishes to optimize his welfare or landscape utility. This welfare depends on two factors: the economic yield of the landscape exploitation, and the ecosystem service provided by the studied species. Using three complementary metapopulation models, we tackle the following questions:

1. What range of exploitation intensity allows the survival of the metapopulation? 2. What proportion of the landscape should the manager exploit to optimize both his yield and the ecosystem service? 3. Can this optimal management be dangerous for the metapopulation in terms of extinction probability? 4. What is the optimal spatial aggregation of non-exploited patches within the landscape? Model (1) is spatially implicit and considers a "presence-dependent" ecosystem service. Model (2) is a two-patch model that investigates the consequences of density-dependence for ecosystem service management ("density-dependent" ecosystem service). Model (3) uses a spatially explicit approach that allows us to investigate how the aggregation of human activities constrains the optimal management.

Materials and Methods

Let us suppose that the economic performance is proportional to the exploited surface, modulated by the ecosystem service. For a manager running such a system, the optimization of the utility is defined by the following equation 1:

max 𝑥 {𝐸 * (𝑥) = 𝑥𝑠 * (𝑥)}
Equation 1

E is the manager's utility, x the proportion of the surface dedicated to the activity and s the ecosystem service produced by the population (all parameters and variables are detailed in Table 1). The symbol * represents the equilibrium of the metapopulation, given the constraint x. The manager's utility follows the cardinal utility theory, such that the utility provided to the manager by several inputs can be quantified (Walras, 1874). This hypothesis of measurability enables the comparison of different utility values. The exploitation activity is assumed detrimental to the species delivering the ecosystem service, so that the manager faces a trade-off between direct production of goods and conservation of the ecosystem service. Increasing the exploited surface harms the metapopulation and therefore reduces the ecosystem service. s* is thus a decreasing function of x. Conversely, if the manager under-exploits the landscape, he reduces the economic gains of the exploitation activity, assumed to be proportional to x. The problem of the manager is to determine the optimal exploited proportion of the landscape x opt that optimizes equilibrium utility E*. The function E should follow the Cobb-Douglas model, often used in economics, which models the relationship between inputs and output welfare [START_REF] Douglas | A theory of production[END_REF]. The two-factor Cobb-Douglas function expression is: y = k 1 x k2 1 x k3 2 , where {x 1 , x 2 } are the factors, k 1 , k 2 and k 3 are constants, here equal to 1. To optimize E, the simplest choice has been made: the optimization of the two-variable product. We detail in the online appendix the effects of changing the Cobb-Douglas function.

We successively deal with two issues: (1) determining the optimal proportion to exploit x, (2) determining the optimal degree of aggregation of exploited patches. The first issue is tackled with spatially implicit models that allow for an entirely analytic solution of the problem, whereas the second requires a spatially explicit model, studied through numerical simulations. The three following models are illustrated on figure 1.

3.1.Spatially Implicit Landscape

We consider a landscape that is spatially implicit. It includes two types of interconnected patches that form a metapopulation -a population of populations [START_REF] Levins | Some demographic and genetic consequences of environmental heterogeneity for biological control[END_REF]. Because the ecosystem service can depend on the presence of the species (presence-dependent ecosystem service) or on the local activity of each individual of the species (density-dependent service), two different metapopulation models are considered: a presence-absence model first [START_REF] Levins | Some demographic and genetic consequences of environmental heterogeneity for biological control[END_REF], then a density-dependent source-sink model [START_REF] Pulliam | Sources, sinks and population regulation[END_REF]. A typical example of a present-dependent ecosystem service is the restoration of an ecosystem induced by the presence of a key engineer species on the landscape [START_REF] Byers | Using ecosystem engineers to restore ecological systems[END_REF]. Alternatively, the provision of fish biomass by designing marine reserves may illustrate density-dependent ecosystem services [START_REF] Costello | Optimal harvesting of stochastic spatial resources[END_REF].

Model (1): presence-dependent ecosystem service

Suppose that the landscape is made of two types of patches, non-exploited (a) and exploited (b). Time variations of patch occupancy can be written (see also figure 1):

� 𝑝 𝑎 ′ (𝑡) = 𝑐(𝑝 𝑎 + 𝑝 𝑏 )(ℎ 𝑎 -𝑝 𝑎 ) -𝑒 𝑎 𝑝 𝑎 𝑝 𝑏 ′ (𝑡) = 𝑐(𝑝 𝑎 + 𝑝 𝑏 )(ℎ 𝑏 -𝑝 𝑏 ) -𝑒 𝑏 𝑝 𝑏 𝑠 ′ (𝑡) = 𝑠 0 �𝛼 𝑝 𝑎 + 𝑝 𝑏 � -𝑚 𝑠

System of equations 2

Assumptions for the first two equations are those of the classical Levins metapopulation model:

• We consider an infinitely spread network of patches.

• A patch can be in two different states: occupied or unoccupied by the population.

• There are two different types of occupied patches: non-exploited (a) and exploited (b). Exploited patches cover a fraction h b of the landscape and the remaining h a is left unexploited (h a + h b = 1).

The fraction h b therefore represents the exploitation strategy of the manager he wishes to optimize.

p a is then the proportion of non-exploited patches occupied by the focal species while p b is the proportion of occupied exploited ones.

• Within each patch, only colonization and extinction processes -at rates c and e respectively -rule the population dynamics.

• We assume that exploitation deteriorates the conditions for the focal species, extinction rates being higher in exploited patches (e a < e b ).

The third equation of the system of equations 2 defines the dynamics of the ecosystem service, given the following assumptions:

• In the case of species extinction (p a and p b = 0), the ecosystem service decreases exponentially toward zero at rate m. The survival of the metapopulation is therefore needed for an optimal management of the landscape; otherwise E tends asymptotically toward zero -see equation 1. Such an exponential decrease allows the modelling of contrasted scenarios. For an ecosystem service that reach zero immediately when the species disappears (e.g., pollination), m can be chosen very high.

Conversely, for an ecosystem service having a certain inertia (e.g. underground temperature regulation by fossorial organism), m can be taken smaller.

• The ecosystem service is provided by the metapopulation at the patch-level throughout the whole landscape. It thus increases with the total proportion of occupied patches p a and p b , at growth rate s 0 .

• To allow for a comparison between the presence/absence model ( 1) and the density-dependent one

(2), the ecosystem service is produced the same way in each case. In model ( 1), the ecosystem service provision thus depends on the presence of the population at the patch-level on the whole landscape, including the non-exploited patches. The relative influence of non-exploited patches compared to exploited patches is expressed by α. This parameter however affects none of our results

qualitatively and we set it to 1 in our analytical computations.

We use this model to solve the optimization problem -equation 1 -for x = h b .

Note that we do not model the reverse interaction between the species and its own niche because model ( 1) is ignoring the feedbacks of the ecosystem service on the engineers. In an alternative model, we consider niche construction processes (appendix 2), where the ecosystem service has an impact not only on the human utility, but also on the engineer species niche.

Model (2), density-dependent ecosystem service

In this second model, the delivery of the ecosystem service is dependent on species density. This model highlights the consequences of density dependence for an optimal management. Time variations of densities can then be written (see also figure 1):

� 𝑁 𝑎 ′ (𝑡) = 𝑟 𝑎 𝑁 𝑎 (1 -𝑁 𝑎 𝐾 𝑎 ⁄ ) + 𝑑 (𝑁 𝑏 -𝑁 𝑎 ) 𝑁 𝑏 ′ (𝑡) = 𝑟 𝑏 𝑁 𝑏 (1 -𝜀𝑁 𝑏 𝐾 𝑏 ⁄ ) + 𝑑 (𝑁 𝑎 -𝑁 𝑏 ) 𝑠 ′ (𝑡) = 𝑠 0 �𝛼 𝑁 𝑎 + 𝑁 𝑏 � -𝑚 𝑠
System of equations 3

Such a system requires the following assumptions:

• The landscape is made of two patches, contrasted in size and / or productivity as modeled through carrying capacities (K a and K b ).

• There are two types of patches: non-exploited (a) and exploited (b). We track the number of individuals in each type of occupied patches: N a in the non-exploited patch and N b in the exploited one.

• Individuals disperse between the two patches at rate d.

• Landscape compositional heterogeneity is modeled by the proportion of land-uses. The effects of the different land-uses on the dynamics are modeled in two ways. The exploited patch is less suitable for the species, affecting intrinsic growth rates (r b < r a with r a > 0) and carrying capacities (K b < K a ). The non-exploited patch is therefore a source in the metapopulation. Conversely, the exploited patch has more immigrants than emigrants being either a true sink (r b < 0) or a pseudo-sink (r b > 0,

K b < K a )
. The landscape has a total carrying capacity K, the non-exploited patch has a carrying capacity K a and the exploited patch a carrying capacity

K b = K -K a .
The third equation determines the dynamics of the ecosystem service according to the following assumptions:

• In case of species extinction, the ecosystem service decreases exponentially at rate m. The survival of the metapopulation is therefore needed for an optimal management of the ecosystem service.

• Contrary to the previous model ( 1), in the present case of a density-dependent ecosystem service, provision depends on the local activity of individuals so that they are produced locally. As in model

(1) however, the ecosystem service may be produced both in the exploited and non-exploited parts of the landscape. The ecosystem service thus increases with the number of individuals in the exploited patch N b and the abundance in the non-exploited patch N a , at growth rate s 0 . The relative influence of the non-exploited patch compared to the exploited patch is expressed by α. This parameter however affects none of our results qualitatively (see appendix 5, which shows the effect of a variation of the parameter α's value). As a result, we set it to 0 in our numerical computations, modelling a situation in which the service produced in the unexploited patch does not diffuse spatially to influence the patch used by the manager.

We use here K b as a proxy of the surface devoted to the exploitation activity and use the model (3) to solve the optimization problem for x = K b / K.

Model (3): spatial aggregation and ecosystem service provisioning

Finally, we study a spatially explicit landscape to tackle how the spatial aggregation of exploited patches affects the optimization problem.

To do so, we apply the system of equations 3 over a n² (n = 30) cell grid to represent a heterogeneous landscape, and study its dynamics in discrete time -see equations 4 and 5, with s 0 and m defined as in the system of equations 3. The grid is set on a torus to avoid edge effects. Cells can be of two different types: non-exploited or exploited. Variables are as defined previously, with i and j coordinates of a cell on the grid, t the time, t* the time at which equilibrium is reached. On odd time steps, reproduction and competition take place, while on even time steps, dispersal takes place (see also figure 1). All technical data and code are given in appendix 3.

Variations in densities can then be modeled:

⎩ ⎨ ⎧ 𝑁 𝑡+1 𝑖,𝑗 = 𝑟 𝑎 𝑁 𝑡 𝑖,𝑗 -𝜅 𝑎 𝑁 𝑡 𝑖,𝑗 2 𝑜𝑜 𝑎 𝑝𝑎𝑡𝑐ℎ 𝑜𝑜 𝑡𝑡𝑝𝑒 (𝑎) 𝑁 𝑡+1 𝑖,𝑗 = 𝑟 𝑏 𝑁 𝑡 𝑖,𝑗 -𝜅 𝑏 𝑁 𝑡 𝑖,𝑗 2 𝑜𝑜 𝑎 𝑝𝑎𝑡𝑐ℎ 𝑜𝑜 𝑡𝑡𝑝𝑒 (𝑏) 𝑁 𝑡+2 𝑖,𝑗 = 𝑑 ′ �𝑁 𝑡+1 𝑖+1,𝑗 + 𝑁 𝑡+1 𝑖-1,𝑗 + 𝑁 𝑡+1 𝑖,𝑗+1 + 𝑁 𝑡+1 𝑖,𝑗-1 � -𝑑 𝑁 𝑡+1 𝑖,𝑗
System of equations 4

𝑠 𝑡+1 = 𝑠 0 �∑ 𝑁 𝑡 𝑖,𝑗 𝑝𝑎𝑡𝑝ℎ𝑒𝑒 𝑜𝑜 𝑡𝑡𝑝𝑒 𝑏 � + 𝑠 𝑡 (1 -𝑚)
Equation 5

• First, populations grow following a logistic equation in each cell. For a given cell a, the growth rate r a is reduced by κ a N a ², a density-dependent competition term.

• Individuals disperse. Dispersal happens at random, in the nearest four neighbor cells.

• Quality of a cell i, q i , has as a positive effect on the growth rate r i of the species and reduces competition coefficient κ i , (equations 6). Such an approach allows for modelling simultaneously growth rates and the carrying capacities based on a single parameter q i .

� 𝑟 𝑖 (𝑞 𝑖 ) = 𝛽 1 ln�𝛽 2 𝑞 𝑖 � 𝜅 𝑖 (𝑞 𝑖 ) = 𝜅 𝑚𝑖𝑚 + (𝜅 𝑚𝑎𝑥 -𝜅 𝑚𝑖𝑚 ) exp (-𝛽 3 𝑞 𝑖 )
System of equations 6

• A favorable environment therefore allows a better reproduction and less intraspecific competition.

• Note that, given system of equations 6: (i) when the quality is close to zero, the growth rate tends toward increasingly negative values (ii) when the quality equals 1, the growth rate is zero, 1 being therefore the limit between source and sink patches.

Quality q b of exploited cells is chosen to mimick different situations: real sink (r b < 0), non-sink

(r b = 0) or pseudo-sink (r b > 0, K b < K a )
. Unexploited cells quality q a always allows positive growth rates, such patches being therefore sources in the metapopulation (figure 6a). Dispersal rate d varies from 0 to 1.

Two spatial components are studied: the proportion h b of the landscape dedicated to exploitation and cell aggregation on the other hand. For a given exploited proportion h b , the (1 -h b ) n² cells are as-signed to unexploited type (a) and h b n² remaining to "exploited" type (b).Landscapes of different levels of aggregation are built (figure 2) following [START_REF] Schlather | An introduction to positive definite functions and to unconditional simulation of random fields[END_REF][START_REF] Schlather | Simulation and Analysis of Random Fields[END_REF]. Cell aggregation is determined by parameters C and π that respectively represent the sill and range of the variogram function γ:

𝛾 : 𝑙 → 𝐶 [1 -exp (-3 𝑙 𝜋)] � Equation 7
The range of the variogram π modulates cell aggregation. The greater the range, the more aggregated the cells of the same type. Conversely, the lower the range, the more fragmented the landscape. Aggregation depends on the exploited proportion h b and aggregation π [START_REF] Griffith | A Casebook for Spatial Statistical Data Analysis[END_REF]. For instance, for an exploited proportion close to 50 %, the aggregation value provides a direct measure of the landscape autocorrelation. The further h b from 0.5 the higher the aggregation, regardless of π.

For each range value π and exploited proportion h b , simulations are run using R until the equilibrium is reached. Equilibrium utility E* is averaged over all the replicates, and we show how E* varies depending on these two parameters (figure 5). We determine which pair (h b ,π) optimizes the utility and then deduce the optimal landscape structure. Therefore here x = {h b ,π} given the optimization problem (equation 1).

Simulation plan:

We adopted an ambitious simulation plan (1 200 000 simulations) that systematically varies the heterogeneity of the landscape and source-sink relationships among exploited and non-exploited patches.

• Quality q b of exploited cells varies in the set {0.01, 1.0, 1.5}.

• Unexploited cell quality q a is 2.

• Dispersal rate d varies between 10 -3 and 1, each intermediate power of 10 being considered

• Parameter π varies between 0 and 1, with a step of 0.02 (50 configurations) and parameter h b varies between 0 and 1, with a step of 0.02 (50 configurations), for a total of 2,500 scenarios (50*50).

• For each scenario of simulations, we did 40 replicates.

This extensive simulation approach allows us to feel confident regarded the robustness of the qualitative results we discuss later.

Results

4.1.Spatially Implicit Landscape

Model (1), presence-dependent ecosystem service

We compute analytically the equilibrium of system 2. Two equilibria exist: a trivial equilibrium {p a * = 0, p b * = 0, s* = 0}, at which the species is extinct, and a non-zero equilibrium, {p a *, p b *, s*}:

⎩ ⎪ ⎪ ⎨ ⎪ ⎪ ⎧ 𝑝 𝑎 * = 𝑝 (𝑒 𝑎 -2 𝑒 𝑏 ℎ 𝑎 -2 𝑒 𝑎 ℎ 𝑏 )+𝑒 𝑎 �𝑒 𝑏 -𝑒 𝑎 +�(𝑝-𝑒 𝑎 +𝑒 𝑏 ) 2 +4 𝑝 (𝑒 𝑎 -𝑒 𝑏 ) ℎ 𝑏 � 2 𝑝 (𝑒 𝑎 -𝑒 𝑏 ) 𝑝 𝑏 * = 𝑝 �𝑒 𝑏 -2 𝑒 𝑏 ℎ 𝑏 +2 𝑒 𝑎 ℎ 𝑏 �-𝑒 𝑏 �𝑒 𝑎 -𝑒 𝑏 +�(𝑝-𝑒 𝑎 +𝑒 𝑏 ) 2 +4 𝑝 (𝑒 𝑎 -𝑒 𝑏 ) ℎ 𝑏 � 2 𝑝 (𝑒 𝑎 -𝑒 𝑏 ) 𝑠 * = 𝑒 0 {𝑝-𝑒 𝑎 -𝑒 𝑏 +�(𝑝-𝑒 𝑎 +𝑒 𝑏 ) 2 +4 𝑝 (𝑒 𝑎 -𝑒 𝑏 ) ℎ 𝑏 } 2 𝑝 𝑚

System of equations 8

By definition, the proportion of exploited occupied patches p b * (resp. p a *) must be positive and not greater than the exploited proportion h b (resp. h a ). The negative root is not considered there, as it makes neither ecological nor economic sense. The ecosystem service s* must also be positive. We also checked analytically the stability of the equilibrium, using the Routh-Hurwitz criterion (details in appendix 4). Such an analysis shows that the two equilibria cannot be simultaneously stable.

As reported on figure 3, two scenarios are to be distinguished. If colonization is greater than extinction (e a < e b < c), then the non-zero equilibrium is stable and the metapopulation maintained (figure 3b). If colonization is greater than extinction only in non-exploited patches (e a < c < e b ) (figure 3a), the non-zero equilibrium is stable if and only if the exploited proportion of the landscape does not exceed a critical value h bc (equation 9). In this second case, exploited patches have a negative effect on the metapopulation survival (c < e b ) and limit the exploited proportion to prevent extinction.

ℎ 𝑏𝑝 = 𝑒 𝑏 (𝑝-𝑒 𝑎 ) 𝑝(𝑒 𝑏 -𝑒 𝑎 )
Equation 9

This critical threshold h bc decreases with extinction rates e a and e b and increases with c. Lower extinction rates and/or higher colonization rates thus allow a higher level of exploitation of the landscape, by making the system more suitable for the species. The critical proportion h bc also increases with the ratio e a / e b . e a / e b is a measure of the similarity between the two types of patches. The closer e a / e b is to 1, the more similar the two types of patches, and the greater the exploited surface can be. In other words, slighter modifications induced by exploitation enable the manager to exploit a larger surface.

We now determine the values of h b opt that are solutions to the manager dilemma (equation 1), i.e. that optimize E* = s*(h b ) h b . Only one solution exists, h b opt :

ℎ 𝑏 𝑜𝑝𝑡 = �(𝑒 𝑎 +𝑒 𝑏 -𝑝) 2 [(𝑝-𝑒 𝑎 ) 2 +(𝑝-𝑒 𝑎 ) 𝑒 𝑏 +𝑒 𝑏 2 ]-[(𝑝-𝑒 𝑎 ) 2 +4 (𝑝-𝑒 𝑎 ) 𝑒 𝑏 +𝑒 𝑏 2 ] 9 𝑝 (𝑒 𝑏 -𝑒 𝑎 )
Equation 10

h b opt decreases with e b and increases with c. Lower impacts of exploitation, and/or higher colonization rate thus favor the maintenance of the metapopulation, and more surface can then be exploited. Increasing e a has contrasted effects on the optimal proportion h b opt .

• Increasing e a increases the similarity between the two types of patches, which increases h b opt .

• Increasing e a reduces the survival of the population, which decreases h b opt .

The first of the two effects is dominant when e a < e b < c (h b opt then increases with e a ). The second effect is dominant when e a < c < e b (h b opt then decreases with e a ).

From comparison of equations 9 and 10, note that the optimal proportion and critical proportion are always such that h b opt < h bc . Setting an optimal management in itself never leads to extinction. However, the optimal management can be dangerously close to extinction threshold h bc . To illustrate this, we compute Δh, the difference between h bc and h b opt :

𝛥ℎ = �(𝑒 𝑎 +𝑒 𝑏 -𝑝) 2 [(𝑝-𝑒 𝑎 ) 2 +(𝑝-𝑒 𝑎 ) 𝑒 𝑏 +𝑒 𝑏 2 ]+𝑝(𝑝-𝑒 𝑎 )[(𝑝-𝑒 𝑏 )²+5 (𝑒 𝑎 -𝑝) 𝑒 𝑏 +𝑒 𝑏 2 ] 9 𝑝 (𝑒 𝑏 -𝑒 𝑎 )
Equation 11

When 𝛥ℎ is small, any mistake in the assessment or implementation of the optimal management strategy, and any environmental stochastic event can lead to extinction. Panel (c) on figure 3 illustrates two managements of the ecosystem service, one riskier than the other. Under conditions 2, the optimal management comes closer to metapopulation extinction than under conditions 1: Δh2 < Δh1.

A sensitivity analysis shows that Δh is minimal for highest |e b -e a | and smallest |c-e a |. Under these conditions, the similarity index e a / e b is low. In other words, when the exploited patches and the nonexploited ones are significantly different in terms of quality, the optimal management may involve large risks for the survival of the metapopulation.

Model (2), density-dependent ecosystem service

The dynamics of the density-dependent model is qualitatively similar to the dynamics of model ( 1).

An analytical solution exists, but it is cumbersome and we prefer here to illustrate the results using a graphical approach (figure 4). As in model ( 1), we always obtain a single strategy K b opt that optimizes utility E. E* is maximal for a given value of K b / K, named K b opt / K. This optimal strategy depends on the source-sink relationships existing among patches.

When the exploited patch is a true sink (figure 4b), the non-zero equilibrium is lower than when the patch is a pseudo-sink (figure 4a): values of N* a , N* b and s* in true sinks are lower than in pseudo-sinks. For true sinks, again in agreement with model (1), a limit of viability exists. Figure 4b represents a case where the extinction limit is not reached. Stable existence requires that K bc > K b , K bc being the critical carrying capacity of the exploited patch above which the metapopulation goes extinct. K bc increases with respect to d and r a . Thus, under low dispersal or low intrinsic growth rate in unexploited patches (a), the metapopulation goes extinct.

Note that pseudo-sink scenarios cannot occur for all exploitation strategies K b (figure 4a). Pseudosinks require that r b > 0, and K b < K a . As K b = K -K a these conditions entail that K b / K < 0.5. When the exploited patch is a pseudo-sink, K b opt = min (0.5 K ; max E*). In the case of a non-sink scenarios (r b = 0), replacing r b = 0 in equations 3, we find s*(

K b ) = (K -K b ) K b s 0 / m, so K b opt = 0.5 K.
The results of a variation of parameter α are provided in appendix 5.

Comparison between the results of model (2) and model (1)

As pointed out by Levins when introducing his classical metapopulation model [START_REF] Levins | Some demographic and genetic consequences of environmental heterogeneity for biological control[END_REF], on which model ( 1) is based, the equations are structurally identical to the classical logistic equation, on which model ( 2) is based. Therefore, posing r i = c -e i and K i =(c -e i )/ c, for i = {a, b} establishes strong links between the two models. Under these conditions, studying the optimal occupancy using model (1) or the optimal carrying capacity division using model ( 2) are equivalent procedures.

A difference however exists due to our definition of ecosystem services in the two models. In model

(1), we considered an ecosystem service produced at landscape scale, due to the presence of the metapopulation ("presence-dependent") while in model ( 2) we assumed that the ecosystem services are produced by local individual activity ("density-dependent"). Equivalence between the two optimization processes therefore requires that we turn off this difference (i.e., setting α in model ( 1) equal to α in model ( 2)).

4.2.Spatially Explicit Landscapes: study optimal aggregation

We compute the utility at equilibrium, thanks to the ecosystem service at equilibrium, which is given by the following expression, according to equation 5:

𝑠 * = 𝑠 0 /𝑚 � � 𝑁 𝑡 * 𝑖,𝑗 𝑝𝑎𝑡𝑝ℎ𝑒𝑒 𝑜𝑜 𝑡𝑡𝑝𝑒 𝑏 �

Influence of the exploited cells quality

The different columns of figure 5 show the effect of varying exploited cell quality on the optimization issue. For quality values lower than 1 (true sinks, left column), the maximum utility (white dot) is obtained for a low aggregation of 0.04 (π opt ) and for an exploited proportion between 0.64 and 0.68 (h b opt ), depending on the dispersal rate and other fixed parameters. Optimally exploited landscapes imply a large exploitation of the landscape and that this activity is spatially fragmented. A linear regression between h b opt and π shows a positive and significant linear relationship between the optimal proportion and the aggregation (figure 5, left column, white lines).

For exploited patches quality values superior or equal to 1 (non-sink and pseudo-sink, middle and right columns of figure 5), two remarkable results must be stressed. First, variations in utility are mainly driven by the area devoted to exploitation (i.e., variations in h b ). Aggregation has comparatively little effect. Second, the optimal space devoted to exploitation h b opt is largely increased in such scenarios (compare with the left column, true sink scenarios). In the results presented on figure 5, this value is even very close to one, i.e. a complete exploitation of the landscape. Note however that such extreme values are reached because only exploited cells produce an ecosystem service. If we allow the non-exploited cells to also contribute to the ecosystem service (appendix 6), h b opt values become less extreme, but still increase with increased quality of exploited cells (increased q b ).

Influence of the dispersal rate

The no dispersal case can be treated analytically. In the case true sink scenario, the population then disappears in exploited cells so the utility is zero whatever the aggregation and the exploited proportion. In the case of pseudo-sinks, the ecosystem service is no longer coming from non-exploited sources cells, so the utility is maximized when exploiting the whole landscape:

h b opt = 1.
Comparing the rows on figure 5, one might think that dispersal has no influence on the results. Indeed, dispersal little impacts the optimal management of the ecosystem service. However, it has a strong effect on the expected level of the manager's welfare E (figure 6b). For true sink scenarios, a strong positive relationship exists between the maximal utility and dispersal (figure 6b, dotted curve). A strong difference indeed exists between exploited cells and non-exploited ones -between sinks and sources, so that local population density and related ecosystem services are closely linked to the contribution of immigrating individuals. Conversely, the higher the exploited cells quality, the lower the effect of dispersal on the maximum utility value (no sink and pseudo-sink scenarios, dashed and solid curves on figure 6b). The higher the quality of exploited cells, the smaller the difference among exploited vs. unexploited patches, the less important the migration effects for the provisioning of the ecosystem service.

Convergence of models (2) and (3) at low aggregation

For large number of patches in the simulation model ( 3), we observe that the associated optimization problem becomes equivalent to the optimization problem associated with model ( 2). This requires that dispersal movements most often link patches of type (a) with patches of type (b), i.e., that the aggregation of the landscape associated with model ( 3) is minimum.

Indeed, consider the total exploited carrying capacity in model ( 3), from which the ecosystem service emerges, given the total capacity of the landscape. This can be written:

𝑥 = 𝑚 𝑏 𝐾 𝑏 𝑚 𝑎 𝐾 𝑎 +𝑚 𝑏 𝐾 𝑏 ,
where n i is the total number of patches of type i and K i their carrying capacity.

Low aggregation only happens when (1) π is small; (2) patches are in approximately equal numbers (i.e., 𝑜 𝑎 ≈ 𝑜 𝑏 ). Under such conditions, x reduces to:

𝑥 = 𝐾 𝑏 𝐾 𝑎 +𝐾 𝑏 ,
which is exactly the optimization problem posed in model (2) (see part 4.1 in Results section).

Discussion

The three models we propose give different, but complementary perspectives on the management of ecosystem services, showing that spatial environmental variations and dispersal levels interact to determine the optimal strategy to be adopted, in terms of the surface to be devoted to the activity, and of the level of aggregation to be used for such activities. The spatially explicit model (3) shows that this compromise is particularly critical when exploited patches are true sinks. Spatially implicit models ( 1) and ( 2) demonstrate the existence of a non-trivial optimal proportion (neither 0 nor 100 %), and show the critical role of variations in quality and in the dispersal rate linking patches (c in the model ( 1), d in the model ( 2)). There is thus a qualitative agreement between the results of models ( 1), ( 2) and ( 3). If we consider more precisely the characteristic pattern of model ( 3) and the E-h b bell curves of the models ( 1) and ( 2), strong links exist between the different models. The projection of utility values at equilibrium E* onto the h b axis (model ( 1)) or K b / K (model ( 2)) is comparable to the pattern drawn on the grid figure 5, for a given aggregation.

Beyond these common points, the three models are complementary in the way they account for spatial components of the dilemma. Considering the two density-dependent models (2) and (3), dispersal processes are not strictly equivalent. Within the spatially implicit model (2), dispersal occurs between two different patches, so all the metapopulation movements are source-sink and not source-source or sinksink. In the spatially explicit model (3), individuals disperse locally. Given a level of aggregation, many movements will link similar patches. The optimal position of sources compared to sinks is then critical to determine the ecosystem service provision. We show that utility is then optimized when the activity is fragmented. Exploited patches are then separated from each other and not combined into a single large patch. This yields large contact limits, which allows a better diffusion of the ecosystem service throughout the landscape, associated with a better maintenance of the metapopulation.

Bioeconomic models of renewable resource exploitation have however pointed out that the population distribution on a patchy landscape influences the effort required to harvest the resource. Considering three patches, one source and two sinks, [START_REF] Sanchirico | Bioeconomics of Spatial Exploitation in a Patchy Environment[END_REF] conclude that the total effort is lower than in the case of three non-interconnected patches. Indeed, as the population disperses towards low profitability patches (sinks), biomass and effort levels are enhanced in the two sink-patches, whereas they decrease in the source-patch in response. Such results can be linked to those of our spatially explicit model (3). Dispersal has proved to increase the gains in the overall landscape in the case of a true sink scenario, for an unchanged effort.

One could question the strong assumption made on dispersal in the model (3), for individuals only disperse in the nearest four neighbor cells. Alternative scenarios are possible, for instance individuals also moving diagonally on the chessboard or moving more than one cell away at each time step. Both scenarios lead to higher dispersal abilities, resulting in a decrease of the environmental grain (MacArthur and [START_REF] Macarthur | The limiting similarity, convergence, and divergence of coexisting species[END_REF], which decreases the effect of aggregation on the optimal landscape management (see figure 5) and modifies the expected ecosystem service (see figure 6). The random movement assumption could also be challenged, as individuals may move preferentially towards non-exploited cells (a). This asymmetry is equivalent to an increase of the relative weight of non-exploited cells over the landscape, which results in a higher optimal proportion of exploited cells.

Our three models may be linked to various empirical studies dealing with the optimal management of ecosystem services. The presence-absence model may be relevant in the case of species that have a broad impact even at low density [START_REF] Byers | Using ecosystem engineers to restore ecological systems[END_REF]. We can then make the approximation that the presence of the species is enough to provide the ecosystem service. To improve productivity (water use) of a desert landscape, porcupines digging pits or ant colonies that form mounds contribute to the creation of hydrological wells are particularly adapted, and efficient even at low densities [START_REF] Byers | Using ecosystem engineers to restore ecological systems[END_REF]. In such situations, our spatially implicit model (1) can provide a useful rule of thumb to determine the proportion of the landscape that should be preserved for these species to act. However, in many instances, the provision of ecosystem services is density-dependent. For instance pollination services will crucially depend on the abundance of the various guilds of pollinators. Regarding pollination, empirical data highlight the possible use of landscape strategies to optimize both the performance of an exploited field and the butterfly abundance [START_REF] Hodgson | Comparing organic farming and land sparing: optimizing yield and butterfly populations at a landscape scale[END_REF]. Because the service measured is tightly linked to butterfly densities, this example illustrates a density-dependent ecosystem service close to models (2) or (3). The study suggests that the intensity of exploitation critically determines the optimal exploitation strategy. When land-use intensification has negative effects on the ecosystem service provided the manager needs to insert natural reserves into his exploitation, highlighting the compromise between biodiversity and production. Alternatively, if the choice of extensive cropping is made, inserting natural reserves at the margins of the field appears efficient enough. We find a similar result in our spatially explicit model: if the quality of an exploited landscape is high enough (pseudo-sink), exploiting the whole landscape leads to an optimal management.

Next to the pollination service, food provisioning from marine reserves rely on fish biomass and is therefore typically linked to densities. Establishing marine reserves allows fish populations to reproduce. Through a spatially explicit model of harvesting a renewable resource (e.g. fishing) results that are consistent with the results of the model (3) can be obtained [START_REF] Costello | Optimal harvesting of stochastic spatial resources[END_REF]):

• To optimize the net value of harvest, one should preserve some areas using biological reserves.

• If biological reserves are connected with anthropic areas, the landscape is heterogeneous in terms of resource (fish) availability, with a possible positive effect on profit. This agrees with results of model (3): more randomly distributed unexploited patches increase overall utility. [START_REF] Costello | Optimal harvesting of stochastic spatial resources[END_REF] consider a different modelling approach (stochasticity affecting both growth rate and the dispersal) but results are still consistent with our model (3).

In terms of agricultural management, our spatially explicit model may be used to discuss the pros and cons of land sharing vs. land sparing. Land sharing proposes that future yields should involve a wildlife-friendly agriculture, with enlarged agricultural areas, while land sparing suggests that the manager dilemma should be solved by maintaining restricted exploited areas, with increased local yields.

Therefore, the two options differ in the proportion of the landscape devoted to agriculture, and in the magnitude of the spatial contrast in source/sink dynamics, two aspects we consider here. [START_REF] Green | Farming and the Fate of Wild Nature[END_REF] suggest that there is not only one simple solution to solve this trade-off. For instance, the optimal type of farming depends on the impact of agricultural yield on the species' density.

If the relationship between population density and crop yield is convex, their model proves that land sparing allows for greater population conservation with a higher yield. This is in agreement with the results of model (3): introducing unexploited islands into an exploited landscape is an optimal management conciliating ecosystem services and metapopulation maintenance. Our results however suggest that a large contact between natural and exploited patches improves this management and that decreasing differences of quality (e.g., pseudo-sink situations) may even allow a full exploitation of the landscape. Our results thus support the view that land sharing may be a better solution than land sparing.

Vandermeer and Perfecto (2005) question some assumptions and conclusions of [START_REF] Green | Farming and the Fate of Wild Nature[END_REF].

They point out the absence of interpatch migration and more generally that interactions among land-uses are ignored. In this regard, the present study helps to bridge such gaps. In the model (1), the colonization process for instance determines patch occupancy, hence the total ecosystem service that is provided. In models ( 2) and ( 3), individuals disperse between the different patches at rate d. Figure 5, stressing the role of dispersal rate on the optimal management, confirms the need to integrate metapopulation aspects into agricultural management, as proposed by Vandermeer and Perfecto.

More generally, in terms of concrete implications, the presented models may help managers to understand the key parameters affecting the optimal exploitation strategy. We assume the landscape manager controls the area he devotes to exploitation and part of its impacts on natural populations in exploited patches (e.g. intensity of the exploitation, temporal pace of crop, use of chemical products, mechanical or manual way of harvesting). The first two models allow a "rule of thumb" determination of the amount of landscape that should be left unexploited depending on the type (presence or density-dependent) of ecosystem service considered. The manager can influence this optimal area by modulating its impacts on patch quality. The third model adds the information that aggregation of natural patches has less effect than the proportion dedicated to exploitation. However, a fragmented distribution of natural patches clearly improves the optimal management. Because we kept the three models simple and general, various limitations affect our results.

From an economic point of view, we based our analysis on a multiplicative utility function E that is simplified in terms of costs. Indeed, we do include a direct cost of exploitation, part of the landscape being attributed to natural habitats and therefore not exploited, but we did not include any extra cost linked to the spatial structure of the landscape. Cost structure however depends on specific features of the ecosystem service or of the landscape considered [START_REF] Bockstael | Ecological economic modeling and valuation of ecosystems[END_REF]. For instance, our results highlight that an integration of natural areas within the exploited landscape is the best solution as it facilitates the diffusion of ecosystem services. However, in an agricultural context, the fragmentation of a field reduces the effectiveness of the exploitation, incurring additional costs we did not consider here. For marine ecosystems, the surveillance of fishing areas becomes more difficult when they are intertwined with marine reserves. If a harvest cost is introduced in the profit function, then if fish density is high, such harvest costs become relatively less important. Remarkably, this constraint also leads to an optimal structure with marine reserves (White et al., 2008).

The three models focus on a single metapopulation that provides a single ecosystem service, or positively correlated ecosystem services. In general, it is expected that ecosystem services would emerge from the entire dynamics of ecological networks, and therefore depends on the community structure and species diversity [START_REF] Mitchell | Linking Landscape Connectivity and Ecosystem Service Provision: Current Knowledge and Research Gaps[END_REF]. Such aspects could be studied using metacommunity models. Either the metacommunity provides a single ecosystem service and species interact with each other, either each species in the metacommunity provides a distinct ecosystem service. In this latter situation, the various ecosystem services are coupled through the multispecies interactions that exist within the network. As a result, trade-offs may emerge between the different ecosystem services, that have to be managed simultaneously [START_REF] Loeuille | Eco-Evolutionary Dynamics of Agricultural Networks: Implications for Sustainable Management[END_REF]. Such optimal management of multiple ecosystem services can significantly improve the utility provided to human beings, especially if the species involved strongly interact (White et al., 2012). They also significantly complexify the optimization approach compared to the relatively simple scenarios we show here.

One of the key assumptions which is made regarding the manager's utility is that the economic performance of the landscape increases with respect to the local ecosystem service produced. Without this critical assumption, the problem raised by this study is no longer relevant. In fact, if the ecosystem service produced has a negative impact on the manager's utility, the problem is solved, the optimal management strategy being the exploitation of the whole landscape. A typical example is the case where the metapopulation present on the landscape fights against the manager's crop. [START_REF] Byers | Using ecosystem engineers to restore ecological systems[END_REF] remind that accidentally introducing an undesirable species on an environment one wants to restore might counter the restoration. The interaction between the manager's welfare and the species function is a keystone of the models.

Conclusions

We model the provision of an ecosystem service by a metapopulation on a more or less anthropized landscape. We solve a double optimization problem: the quantity of ecosystem service provided by the metapopulation and the yield of the human manager's exploitation i.e. the exploited proportion of the landscape. When it comes to optimizing both the manager and the metapopulation welfare, it necessarily leads to a trade-off. First we find that overall benefits are optimized for intermediate surfaces of exploitation, meaning that there is no trivial solution. Then we show that this trade-off depends on ecological traits and on the degree of human disturbance on the landscape. At last, we prove that the optimal management is improved when human activities are spatially fragmented. The utility scale is provided on the bar on the right of each grid; the darker the colour the higher the utility at equilibrium. The white dotted line is the regression of the maximum utility values for each value of aggregation. The white lozenge is the maximum value on the grid, that is to say the optimal mode of management (hb opt ,πopt). Solid line, qb = 1.5, dashed, qb = 1, dotted qb = 10 -2 . Table

Table

Fig. 1 .

 1 Fig. 1. Conceptual diagram of the three different models. The first two models are illustrated on subplots a and b. The third one is presented on subplot c.

Fig. 2 .

 2 Fig. 2. Representation of a range of nine landscapes according to hb, the exploited proportion of the landscape and π, the level of aggregation of patches of the same type. Toric landscapes are shown in two dimensions. Each landscape is represented with a square of 30 patches side. In light grey, the non-exploited patches, in dark grey the exploited patches.

Fig. 3 .

 3 Fig. 3. Representation of the utility at equilibrium E* with respect to the exploited proportion hb. On (a) and (b), illustration of two solutions of the optimization problem. On (c), illustration of the compromise between the manager's welfare and the engineer conservation. The utility at equilibrium E* is plotted with respect to hb for m = 0.15 and s0 = 0.70. hbc is represented with a dotted line.(a) ea < c < eb. c = 0.20, ea = 0.05 and eb = 0.50. hb opt = 0.45. (b) ea < eb < c. ea = 0.30, c = 0.50 and eb = 0.45. hb opt = 0.68.(c) The arrows represent Δh. In grey, the optimal management is riskier than in black. Black: c = 0.20, ea = 0.05 and eb = 0.50. Distance Δh1 is equal to 0.38. Grey: c = 0.15, ea = 0.10 and eb = 0.70. Δh2 is then equal to 0.19.

Fig. 4 .

 4 Fig. 4. Illustration of numerical results of model (2). In each case, isoclines {Na, Nb} are represented on the left panel -the intersection illustrates the positive stable equilibrium of the system of equations. In dotted, Nb, in solid, Na. The same scale is used on X-axis and Y-axis. Unity is individuals. The utility function at equilibrium with respect to Kb / K is represented on the right panel. ra = 0.4, d = 0.1 (a) Pseudo-sink, rb = 0.1. In grey is represented the zone where Kb / K cannot take any value. (b) True sink, rb = -0.1. In every subplot, α = 0.

Fig. 5 .

 5 Fig. 5. Grids of results representing the utility at equilibrium E* with respect to the quality of the exploited patches qb and the dispersal rate d. qa is set to 2, s0 = 0.80 and m = 0.30. Data represented on each grid are utilities at equilibrium with π the aggregation in abscissa and hb the exploited proportion of the landscape in ordinate. The utility scale is provided on the bar on the

Fig. 6 .

 6 Fig. 6. (a) Representation of the competition coefficient and the intrinsic growth rate of the population with respect to the quality of a patch. Increasing curve: the intrinsic growth rate (left scale). Decreasing curve: competition coefficient (right scale). The parameters chosen are those used in the simulations: β1 = 1, β2 = 1, β3 = 0.5 log (log (2) -0.2) / 0.8, κmin = 2.10 -4 and κmax = 1.10 -3 . (b) Illustration of the influence of dispersal on the maximum value of utility at equilibrium. Representation of the maximum utility with respect to the dispersal rate, for three values of the quality of the exploited patches. Logarithmic scales.

Fig. 1

 1 Fig. 1

Fig. 2

 2 Fig. 2

Fig. 3

 3 Fig. 3

  Fig. 4

  Fig. 6

  Fig. S2

  Fig. S3

  Fig. S4

Fig. S5

  Fig. S5

Table 1 . Table of parameters and variables of the optimization problem and the different models.

 1 Tylianakis, J.M., Tscharntke, T.,Lewis, O.T., 2007. Habitat Modification Alters the Structure of Tropical Host-parasitoid Food Webs. Nature. 445 (7124), 202-205. doi:10.1038/nature05429. Vandermeer, J., Perfecto, I., 2005. The Future of Farming and Conservation. Science. 308, 1257. Verhulst, P.F., 1838. Notice sur la loi que la population poursuit dans son accroissement in: Correspondance mathématique et physique, vol. 10, pp. 113-121. Walras, L., 1874. Éléments d'économie politique pure, ou théorie de la richesse sociale. Watkinson, A.R., Sutherland, W.J., 1995. Sources, sinks and pseudo-sinks. Journal of Animal Ecology. 64, 126-130. White, C., Kendall, B., Gaines, S., Siegel, D., Costello, C., 2008. Marine reserve effects on fishery profit. Ecology Letters.11, 370-379. doi:10.1111/j.1461-0248.2007.01151.x. White, C., Costello, C., Kendall, B.E., Brown, C.J., 2012. The value of coordinated management of interacting ecosystem services.

	Ecology Letters.15, 509-519. doi: 10.1111/j.1461-0248.2012.01773.x.

1

  Variables of the optimization problem Equal to 0.5 log (log(2) -0.2) / 0.8 in all the simulations𝜅 𝑚𝑖𝑚Equal to 2. 10 -4 in all the simulations 𝜅 𝑚𝑎𝑥 Equal to 1. 10 -3 in all the simulations

	𝛽 1	Equal to 1 in all the simulations
	𝛽 2	Equal to 1 in all the simulations
	𝛽 3
	𝐸	Manager's welfare
	𝑥	Exploited proportion of the landscape
	𝑠	Ecosystem service
	Variables and parameters of the model (1)
	𝑝 𝑎	Proportion of occupied non-exploited patches
	𝑝 𝑏	Proportion of occupied exploited patches
	𝑐	Average rate of colonization events
	𝑒 𝑎	Average rate of extinction events in a non-exploited patch
	𝑒 𝑏	Average rate of extinction events in an exploited patch
	ℎ 𝑎	Proportion of non-exploited patches on the landscape
	ℎ 𝑏	Proportion of exploited patches on the landscape
	𝑠 0	Growth rate of the ecosystem service
		Decay rate of the ecosystem service in the absence of the
	𝑚
		engineer
		Relative influence of non-exploited patches compared with
	𝛼
		exploited patches in the ecosystem service
	Additional variables and parameters of the model (2)
	𝑁 𝑎	Number of individuals in the non-exploited patch
	𝑁 𝑏	Number of individuals in the exploited patch
	𝑟 𝑎	Intrinsic growth rate in the non-exploited patch
	𝑟 𝑏	Intrinsic growth rate in the exploited patch
	𝐾	Carrying capacity in the whole landscape
	𝐾 𝑎	Carrying capacity in the non-exploited patch
	𝐾 𝑏	Carrying capacity in the exploited patch
		Equals +1 if rb > 0; equals -1 otherwise (negative density-
	𝜀
		dependence)
	𝑑	Metapopulation dispersal rate
	Additional variables and parameters of the model (3)
	𝜅 𝑎	Competition within an non-exploited patch parameter
	𝜅 𝑏	Competition within an exploited patch parameter
	𝑑′, 𝑑	Dispersal rates


Ana Isabel Perez Jimenez

Lylian Challier

Margherita Di Pisa

Manon Guille-Collignon

Frédéric Lemaître

Solange Lavielle

Christelle Mansuy

Christian Amatore

Eric Labbé

Olivier Buriez

Margherita Di Pisa

Manon Guille- Collignon

Three-Electrode Analytical and Preparative Electrochemistry in Micro-Volume Hanging Droplets

Keywords: Electroanalysis, Micro-volume, Droplet, Platinization, Cell-penetrating peptide

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction.

Miniaturization of analytical techniques is an important issue for samples of low availability and for the development of portable devices, as well as for evident economic and ecological reasons (consumption of small amounts of precious and/or expensive reagents, waste production). Within this framework, electrochemistry brings a significant contribution to the analysis in small volumes, notably through the miniaturization of electrodes [1]. On the other hand, the development of microfabrication techniques such as photolithography [2][3][4][5][6][7], screen printing [8,9], plasma etching [10], and laser ablation [11,12] allowed the preparation of volume-limited structures. In this context, two main approaches have been developed aimed at running electrochemistry in small volumes. The first one, essentially developed in view of the detection of small numbers of molecules released by single cells, consists in the immersion of needle-type micro-electrodes in small chambers, or vials, prepared by micro-fabrication techniques [2,3,13,14]. Alternatively, the second route deals with the immersion of needle-type micro-electrodes within droplets, or meniscus, thus avoiding the preparation of containers [15,16,17]. This approach has been implemented with micro-, nano-, picoand femto-liter droplets deposited either on a substrate or onto the electrode itself.

Micro volume droplets are generally deposited onto the working electrode surface [18,19]. However, comparatively large liquid-solid reference electrode (e.g. SCE) cannot be immersed in such drops. Instead, pseudo-reference wire-shaped electrodes such as Ag, Ag/AgCl, or Pt are used though they suffer from potential instability with some species [20,21,22].

Herein, a new approach is proposed for real three-electrode operation with a commercial SCE in micro-volume droplets combined to smaller counter and working electrodes. This allows investigating both aqueous and organic solutions provided a chamber is used to prevent evaporation. Actually, the novelty of this device relies on the fact that droplets are deposited neither on a substrate nor on the working electrode surface, but suspended at the fritted glass of a bridge compartment in which the SCE is placed (Figure 1A). The bridge was aimed at preventing any contamination of the reference electrode from the solution and vice versa. This original electrochemical cell was shown to be perfectly adapted to electroanalysis of low quantities of highvalued molecules and also for controlled-potential electrolysis.

Figure 1.

(A) Experimental device used to perform cyclic voltammetry and controlled-potential electrolysis in 50 µL hanging droplets. Note that a U-shape bridge compartment could be also used in order to have the droplet deposited on the fritted glass surface rather than hanging from the latter.

(B) Voltammetry of ferrocenemethanol (0.5 mM) in H 2 O + KCl (0.1 M) at a 10 µm diameter carbonfiber electrode (v=20 mV/s). The same curve (1) was obtained in the droplet or in bulk solution or after cleaning the electrode surface. Curves (2) to (6) corresponded to successive scans run after 2, 3, 8, 13, and 15 minutes, respectively, in the micro-volume cell. Inset shows a voltammogram obtained at 20 V/s in the droplet.

(C) Voltammetry of ferrocene (1 mM) in DMF + nBu 4 PF 6 (0.1 M) at a platinum electrode after 15 minutes in the droplet (dotted line) and in a bulk solution (solid line). v= 20 mV/s. Inset shows a voltammogram obtained in the droplet at 20 V/s. (D) Voltammetry of (a) ferrocenemethanol (25 µM) and of (b) Fc-(Arg) 9 (50 µM) in PBS at a 10 µm diameter carbon-fiber electrode in the droplet (v=20 mV/s).

-

Cyclic Voltammetry in Water and in Organic solvent

Cyclic voltammograms of ferrocenemethanol were first recorded at a carbon fiber electrode immersed in a suspended 50 µL water droplet to test and validate the device under steady-state and transient conditions. Usual sigmoidal-and peak-shaped voltammograms were obtained at slow and high scan rates, respectively (Figure 1B). Identical voltammograms were obtained in a bulk solution (3 mL) of same composition upon using the same electrodes demonstrating that, in microliter volumes, voltammetric responses are not affected by depletion phenomena as observed in pico-liter cells [15]. Furthermore, it was calculated that 4.82 mC would be required to completely electrolyze 50 µL of 1 mM ferrocenemethanol. Considering a steady-state current of ca. 4 nA flowing during 25 seconds, i.e., ca. 0.1 µC passed, only about 2.10 -3 % of the analyte was oxidized during one slow scan. This confirms that depletion of ferrocenemethanol in the drop can be disregarded even at slow scan experiments.

Another important point to test was evaporation. It was already reported that evaporation from an air-exposed 1 µL aqueous droplet was significant after around 60 seconds (peak currents increased by more than 5%) [23]. Since relative evaporation rate is proportional to the surface-to-volume ratio of the sample, no evaporation, and therefore no peak current increase, was expected under our conditions within tens of minutes. Conversely, a very slight peak current decrease was observed with time (Figure 1B). This was due to the fouling of the electrode surface. Indeed, the initial peak current value was recovered after polishing the electrode (Figure 1B), and no evidence of evaporation was observed after 15 minutes exposure to air.

Organic solvent droplets could be suspended as well to the bridge fritted glass. Cyclic voltammetry of ferrocene was therefore performed in dimethylformamide or acetonitrile with a platinum microelectrode (the shaft of the above carbon-fiber was insulated by poly(oxyphenylene) which may dissolve in organic solvents). To prevent solvent evaporation, a tissue soaked with the same organic solvent was placed into the measurement chamber. An argon stream saturated with solvent was also gently blown before starting each experiment, for five minutes, to remove oxygen from the chamber. Transient and stationary voltammograms were recorded at high (20 V/s) and slow (20 mV/s) scan rates, respectively (Figure 1C). As observed in water, the electrochemical responses were identical to those obtained in a bulk solution (3 mL) and no solvent evaporation was noted during 15 minutes at least, in agreement with previous reports [23].

Controlled-potential electrolysis

The presence of a saturated calomel electrode in our device is particularly relevant to perform accurate controlled-potential electrolysis in micro volumes. Since the anode and cathode were in close proximity (few millimeters) any electrolysis of a reversible redox soluble couple would have led to redox cycling. For this reason, we considered the reduction of H 2 PtCl 6 or KAuCl 4 in aqueous solutions as test experiments. A series of controlled-potential electrolyses were run in suspended 50 µL water droplets, instead of usual milliliter volumes [24], as a function of the charge passed and analyzed by SEM (Figure 2). Platinization of the carbon-fiber microelectrode was successfully achieved. SEM images show that platinization occurred only on the active surface but not on the insulated one. Also, no interference of the counter electrochemical reaction was observed within the short time-scale electrodeposition (estimated maximum diffusion layer thickness c.a.: 0.5 mm). A typical variation of the amperometric current following the reduction of hydrogen hexachloroplatinate(IV), in the presence of lead acetate, on the carbon-fiber active surface is shown in Figure 2(IV)-left.

Similarly, electrodeposition of gold nanocrystals could be achieved onto carbon-fiber electrode surfaces. In Figure 2-right are shown the SEM images obtained after reduction of potassium tetrachloroaurate for 30 seconds at either + 400 or -500 mV / SCE corresponding to the plateau and the foot of KAuCl 4 reduction wave, respectively. As already demonstrated on glassy carbon electrodes [25], the size, structure, and distribution of gold nanocrystals depended on the overpotential, illustrating the benefit of the three-electrode cell operation.

Electrochemical detection of a high value-added molecule

Owing to its minute volume our method is perfectly suited for electroanalysis of precious/poorly available molecules. This was illustrated through the detection of a cationic cell-penetrating peptide (CPP) grafted to a ferrocenyl (Fc) moiety. Presently, CPPs are among the most popular and efficient shuttle to access the intracellular medium owing to their ability to penetrate cell membranes at micromolar concentrations without causing significant membrane damages while carrying biologically active cargoes (drugs, oligonucleotides, proteins, nanoparticles, etc.) [26,27]. For these reasons, they have opened new avenues for therapeutic applications, notably in cancer therapy [28,29]. Though fluorescence-based techniques are the most commonly exploited approaches to detect these peptides [30,31], electrochemical methods may also be of great analytical interest. In this context, we recently reported on the modification of a cationic CPP with a ferrocenyl unit [32]. However, these experiments required large CPP amounts (milliliter electrochemical cell and millimetric sized electrodes). Hence, the electrochemical detection of a CPP composed of nine arginines and tagged with a ferrocenyl unit (noted as Fc-(Arg) 9 ) at a carbon fiber electrode immersed in a micro-volume suspended water droplet was investigated. As shown in Figure 1D, a well-defined sigmoidal voltammogram was obtained at slow scan rate in the presence of micromolar concentration of the Fc-(Arg) 9 peptide. Also, the Fc-(Arg) 9 oxidation appeared at a more positive potential value than that of ferrocenemethanol illustrating an electron-withdrawing effect of the amino acids [33]. This again illustrated the analytical interest of the present three-electrode droplet cell approach.

Conclusion

This work established the validity of electrochemical investigations performed in micro-volume hanging droplets, in the presence of a conventional reference electrode. This original approach is simple, low-cost, easy to implement and versatile. Moreover, it does not require any construction of a miniaturized electrochemical cell, thus affording immediate electrochemical characterization.

Experimental Section

Chemicals. Tetra-n-butylammonium hexafluorophosphate (Acros Organics) and Potassium chloride (Sigma Aldrich) were used as supporting electrolytes in organic solvents and in highly purified water (HPW; resistivity = 18 M.cm; Milli-Q system; Millipore, Billerica, MA, USA), respectively. Acetonitrile (Sigma Aldrich), N,N-dimethylformamide (Acros Organics), ferrocene (Fluka), ferrocenemethanol (Aldrich), potassium tetrachloroaurate (Aldrich), sulfuric acid (Sigma-Aldrich), hydrogen hexachloroplatinate (Aldrich), and lead(II) acetate trihydrate (Sigma) were used as received. Phosphate-buffered saline (PBS) pH 7.4 (Gibo, Invitrogen Co.) was prepared with HPW.

Instrumentation.

Cyclic voltammetry and controlled-potential electrolysis experiments were performed in a Faraday cage using a Princeton Applied Research potentiostat (PARSTAT 2273). The quality of gold and platinum depositions were monitored by scanning electron microscopy (Eline model by RAITH).

Working Microelectrode Fabrication. Carbon fiber microelectrodes (10 µm diameter, Thornel Carbon Fibers, Cytec Engineered Materials, Greenville, SC, USA) were constructed, insulated and polished as described previously [34].

Platinum working disk-electrodes were made through cross-sections of platinum wires (25 µm diameter, Goodfellow) sealed in a Pyrex glass capillary (Harvard Apparatus, length: 10 cm; outer diameter: 1.2 mm; inner diameter: 0.6 mm). Electrodes were polished to a mirror finish with emery paper of decreasing grain size and alumina powder (Minitech 233 by PRESI), sonicated for 10 minutes in HPW, and dried before use.

Electrochemical Experiments. All measurements were made at room temperature. The reference electrode was a commercial saturated calomel electrode (SCE -Radiometer Analytical) equipped with a bridge compartment filled with the same solvent/supporting electrolyte solution as in the drop. The counter electrode was a 1 cm length platinum wire (0.5 mm in diameter, Goodfellow). The working and counter electrodes were held with the help of H-7 pipette holders (Narishige group). Their immersion in the drop was done manually without the help of any micromanipulators. Water and organic solvents droplets (50 µL) were dispensed with a Hamilton microsyringe directly on the fritted-glass of the bridge compartment. To avoid any solvent evaporation and/or when an inert atmosphere was required, the drop and electrodes were placed in a glass enclosed chamber in which an argon stream saturated with the corresponding solvent was gently passed during five minutes before starting the experiment. A piece of tissue (TX 604 cellulose/polyester, surface: 2 cm 2 ) soaked in the appropriate solvent was also placed in the chamber. [35,36].

Solid Phase Peptide Synthesis. All standard chemicals and protected amino acids were purchased from Sigma, Merck chemicals Ltd., Carlo Erba or Iris Biotech AG.

The Arg 9 sequence was synthesized on a 0.1 mmol scale by conventional tBoc/Tosyl SPPS starting from a MBHA-Resin (0.54 mmol/g loading) on an 431A ABI synthesizer (Applied Biosystem). Boc-Arg(Tos) (1mmol) was coupled with DCC/HOBt (1eq., 1:1). The dried Boc-peptidyl resin was transferred into a fritted syringe, swollen (NMP) and reacted with ferrocene carboxylic acid (0.5 mmol), HBTU (0.45 mmol), and DIPEA (1 mmol), in NMP for 30 min. The ferrocene peptidyl-resin was washed (NMP, DCM, MeOH), and dried under vacuum. Then, the peptide was cleaved from the resin with HF at 0°C, in the presence of anisole (1.5 μL/mg of resin), for 2 hours. Elution with a 10% acetic acid solution afforded the crude peptide, which was purified by RP-HPLC (

Figure 2 .

 2 Figure 2. Left: Scanning electron microscopy images of a carbon fiber (10 µm in diameter) electrode obtained before (I) and after a platinization charge of 20 µC (II) and 40 µC (III) within the droplet. In (IV) is shown a typical evolution of the amperometric current recorded during the reduction (E = -65 mV/SCE) of hydrogen hexachloroplatinate(IV) (28 mM) in the presence of lead acetate (0.57 mM).The electrodeposition started at t = 0 and was pursued until reaching a platinization deposit (30 µC in this case); inset shows the corresponding platinum deposit.Right: Scanning electron microscopy images of a 10 µm gold modified carbon fiber electrode by reduction of potassium tetrachloroaurate(III) (2.6 mM) dissolved in sulfuric acid (0.1 M). Images (I) and (II) were obtained after reduction at + 400 mV/SCE for 30 seconds. (III) and (IV) were obtained after reduction at -500 mV/SCE for 30 seconds.

Preparation of Gold-Coated and Platinized Carbon-Fiber Electrodes through Electrolysis

  SCE. Platinized carbon-fiber electrodes were obtained by reducing, at -65 mV vs. SCE, a solution prepared itself by mixing 1 mL of hydrogen hexachloroplatinate(IV) (8 wt.% solution in water) and lead(II) acetate trihydrate (99.8 %, 1.6 mg) in PBS(6.4 mL)

. All polished carbon-fiber surfaces were modified by controlled-potential electrolyses after immersion of the electrode in a 50 µL droplet of the appropriate solution. Deposition of gold metal at the surface of carbon-fiber electrodes was achieved by reduction of potassium tetrachloroaurate(III) (2.6 mM) dissolved in sulfuric acid (0.1 M), at + 400 or -500 mV vs.

  5-30% B in 30 minutes, VP 250 x 16 mm NUCLEODUR ® C18 HTec 5 μm column working at 14 mL/min, with the solvent systems: A (0.1% TFA in H 2 O) and B (0.1% TFA in CH 3 CN), and UV detection at 220 and 280 nm. Characterization of Fc-(Arg) 9 : i) Analytical RP-HPLC: 8.48min, purity over 95% (5-30% B in 10 minutes, Proto 200 C18-3 μm 100 x 4.6 mm column at 1 mL/min (Higgins Analytical)); ii) MALDI-TOF MS (DE-Pro, ABI) in positive ion reflector mode using the CHCA matrix (m/z): [M+H] + calcd 1635.73; found, 1635.09.

Acknowledgements. We thank J. Palomo for the SEM micrographs. This work was supported by the CNRS (UMRs 8640 and 7203), the Ecole Normale Supérieure, Université P. et M. Curie, and the Agence Nationale de la Recherche (grant number: ANR-12-BS08-0002-01 "ELIPTIC"). The authors also thank Conacyt for support of A.I.P.J.