Nicolas Parisey
email: nicolas.parisey@inrae.fr

Melen Leclerc
email: melen.leclerc@inrae.fr

Katarzyna Adamczyk-Chauvat
email: katarzyna.adamczyk@inrae.fr

K Adamczyk

Optimal spatial monitoring of populations described by reaction-diffusion models

Keywords:

Using spatialised population measurements and related geographic habitat data, it is feasible nowadays to derive parsimonious spatially explicit population models and to carry on their parameter estimation. To achieve such goal, reaction-diffusion models are common in conservation biology and agricultural plant health where they are used, for example, for landscape planning or epidemiological surveillance. Unfortunately, if the mathematical methods and computational power are readily available, biological measurements are not. Despite the high throughput of some habitat related remote sensors, the experimental cost of biological measurements are one of the worst bottleneck against a widespread usage of reaction-diffusion models. Hence we will recall some classical methods for optimal experimental design that we deem useful to spatial ecologist. Using two case studies, one in landscape ecology and one in conservation biology, we will show how to construct a priori experimental design minimizing variance of parameter estimates, enabling optimal experimental setup under constraints.

Introduction

Both empirical and theoretical studies have well established that population spread is an essential ecological process for understanding most of the observed population dynamics. Consequently, over the last two decades spatial ecology has become central either in theoretical ecology or in empirical approaches. Hence, population management for their control, capture or protection, has started to be considered at large spatial scales that match the inherent dispersal capacity of the considered species. Populations can be identified and monitored using various strategies (e.g. tracking, trapping). However, despite the development of new technologies and devices for ecological monitoring (e.g. imaging sensors for animal detection [START_REF] Weinstein | A computer vision for animal ecology[END_REF], autonomous Unmanned Aerial Vehicle [START_REF] Cliff | Robotic ecology: Tracking small dynamic animals with an autonomous aerial vehicle[END_REF]), the survey of many species on large areas remains challenging, costly and empirically guided.

The growing interest for space in ecology has been accompanied by a proliferation of models and statistical methods for analyzing and predicting the spatial distribution of populations. Nowadays, existing statistical methods allows one to infer the dynamics of spreading populations from noisy and sparse data. However, statistical inference of spatiotemporal models from common monitoring data can be subject to practical identifiability issues and is often associated with an important uncertainty on parameters estimates. Thus, given the inherent difficulty of collecting spatial data on most population, model-based forecasts are still associated with an important uncertainty, which could theoretically be reduced by more dense and efficient monitoring strategies.

When planning a survey of population over a large area, one has to make some choices giving financial and human constraints: which observational method (transects, traps), where and when to start the survey, the number and times for repeating the survey. These questions can either be addressed empirically or using model-based techniques that can be organized into three major groups of methods [START_REF] Ucinski | Optimal measurement methods for distributed parameter system identification[END_REF]. The first group transform the problem into a state-estimation one by augmenting the state vector [START_REF] Malebranche | Simultaneous state and parameter estimation and location of sensors for distributed systems[END_REF]. The second group is based on random fields theory [START_REF] Sun | Inverse Problems in Groundwater Modeling[END_REF] whereas the third group corresponds to the classical theory of the design of experiments [START_REF] Ucinski | Optimal measurement methods for distributed parameter system identification[END_REF]. We can also point out a fourth group that consists in adopting an heuristic to find the best sampling strategies among tested situations [START_REF] Bellot | Inferring ecological processes from population signatures: A simulation-based heuristic for the selection of sampling strategies[END_REF]. Based on either of these groups, methods for designing ecological sampling has already been addressed by several authors [START_REF] Hooten | Optimal spatio-temporal hybrid sampling designs for ecological monitoring[END_REF][START_REF] Williams | Monitoring dynamic spatio-temporal ecological processes optimal ly[END_REF] but it has not yet percolated in the modellers community and is still seldom considered in population ecology.

The use of optimal design of experiments is perhaps one of the lesser known [START_REF] Cook | Optimal observation times in experimental epidemic processes[END_REF] whereas it is well developed for optimal sensor placements in spatially-distributed systems with non-linear dynamics models [START_REF] Ucinski | Optimal measurement methods for distributed parameter system identification[END_REF]. It corresponds to an entire branch of statistics that provides criteria measuring the amount of information about unknown model parameters carried out by the observed data.

Here, "experiments" is taken sensus lacto and refers to controlled observations of populations within their living areas that are planned given the constraints of available resources and with subsequent statistical analysis in mind.

In this study we use the optimal design of experiments to tackle the question of ecological monitoring of populations described by reaction-diffusion models. This type of partial differential equations originally emerged in chemistry for analyzing the change in space and time of chemical substances before becoming one of the most important mathematical model for the study of spreading dynamical processes in biology and ecology. Here, we consider reactiondiffusion equations that provides a parsimonious description of a spreading population in a domain Ω included in R 2 :

( , ) = ( )Δ ( , ) + ( , , ) (1)
where ( , ) is the population density at time and location ∈ Ω, ( ) corresponds to the dispersion rate at , Δ ( ,

)
is the Laplace operator of evaluated at ( , ) that describes the random movement of individuals, and ( , , ) is a general reaction function accounting for local birth, death and interactions within the spreading population. To fully describe the system, one needs to define some boundary conditions on Ω as well as initial conditions 0 ( ) = (0, ).

In addition to their ability to capture the essential observed patterns of spreading populations, reaction-diffusion systems have been well studied by mathematicians for a long time. Thus, though the mathematical and the numerical analysis of some reaction-diffusion equations still challenge mathematicians, several systems are now well characterized and understood. The most popular one is the Fisher-KPP system, that applies in ecology, for which ( ) = and the reaction term accounts for a limiting carrying capacity of the environment ( , , ) = (1-). Parameter estimation of mechanistic PDE models from noisy observational data, obtained during common population monitoring strategies, can be performed using different statistical methods. In any case, one needs to link the mechanistic PDE model that describes the changes in space and time of the population density with a probabilistic model describing the observation process. This approach refers to the physical-statistical or mechanistic-statistical models that are specific instances of the more general hierarchical state-space model framework [START_REF] Clark | Population time series: process variability, observation errors, missing values, lags, and hidden states[END_REF][START_REF] Soubeyrand | Parameter estimation for reaction-diffusion models of biological invasions[END_REF][START_REF] Wikle | Hierarchical models in environmental science[END_REF]. Such approach has been successfully used for inferring spatio-temporal processes on large areas for invasive and beneficial insects [START_REF] Parisey | Rearranging agricultural landscapes towards habitat quality optimisation: In silico application to pest regulation[END_REF][START_REF] Roques | A statistical-reaction-diffusion approach for analyzing expansion processes[END_REF][START_REF] Roques | Modelling population dynamics in realistic landscapes with linear elements: a mechanistic-statistical reactiondiffusion approach[END_REF], plant diseases [START_REF] Abboud | Dating and localizing an invasion from post-introduction data and a coupled reaction-diffusion-absorption model[END_REF] or aquatic and terrestrial mammal species [START_REF] Louvrier | A mechanistic-statistical species distribution model to explain and forecast wolf (canis lupus) colonization in south-eastern france[END_REF]. As illustrated by [START_REF] Williams | Monitoring dynamic spatio-temporal ecological processes optimal ly[END_REF], the optimal design of experiments can be applied on mechanistic-statistical systems for ecological monitoring,

but is yet poorly considered by both modellers and ecologists collecting spatio-temporal data.

In this article we first introduce the main concepts of the optimal design of experiments applied on mechanisticstatistical systems for population dynamics. Then, we consider two example populations whose spatio-temporal dynamics is described by reaction-diffusion equations, but monitored with different strategies: beneficial insects for agriculture that are monitored using Barber pitfall traps at various locations within the landscape [START_REF] Parisey | Rearranging agricultural landscapes towards habitat quality optimisation: In silico application to pest regulation[END_REF], and invasive wild horses that are counted during aircraft transects [START_REF] Beeton | Modelling horse management in the Australian Alps[END_REF]. Assuming that the more realistic degree of freedom in the existing monitoring designs is the choice of the locations of observations we consider static designs, where locations of sampled populations are fixed before the start of the survey and don't change in time, and focus on the spatial aspect of monitoring. We use D-optimal designs that define optimal locations for population observation and show how this framework can help to design ecological monitoring strategies over large areas under identical constraints. We finish the paper by discussing the use of the optimal design of experiments framework to improve the monitoring of populations over large areas, especially to reduce the cost of sampling and support the link between modelling and empirical studies.

Mechanistic-statistical model for population dynamics

We consider a mechanistic-statistical approach that combines a reaction-diffusion system which describes essential spatio-temporal dynamics of a population with an observational equation describing the stochastic process leading to detection and enumeration of individuals at location ∈ Ω and time (as in [START_REF] Soubeyrand | Parameter estimation for reaction-diffusion models of biological invasions[END_REF]). However, following the framework used in control system analysis or optimal design for parameter identification [START_REF] Ucinski | Optimal measurement methods for distributed parameter system identification[END_REF], we decompose the observation process into two parts: 1) the measurement process which links the ecological process and factors that affect the measurement of population density (e.g. equipment features, operator behaviour, environmental effect), and 2) the data process that links population data with the measurement process. Then, given some boundary conditions on Ω the spatially distributed system is described by the following hierarchical system :

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ( , 0) = 0 ( ) ( , ) = Δ ( , ) + ( , , , ) ( , ) = ( ( , ), ) ∼( ( , ), ) (2)
where ( , ) is the measurement process determined from population density ( , ) with function (.), and the data process corresponds to the draw of population data from a general probability distribution . Depending on the considered system, parameters that needs to be estimated, or identified, from the monitoring data can occur in the population process , , the measurement process ( ) and the data process ( ). The set of unknown parameters is thus given by = , , , . Albeit initial conditions 0 ( ) of reaction-diffusion systems may also be estimated from noisy data [START_REF] Abboud | Dating and localizing an invasion from post-introduction data and a coupled reaction-diffusion-absorption model[END_REF][START_REF] Soubeyrand | Parameter estimation for reaction-diffusion models of biological invasions[END_REF], in this study we assume that 0 ( ) is known and fixed. Parameter estimation can be achieved using maximum likelihood estimation (see next section). It generally involves a numerical optimization that requires the simulation of the reaction-diffusion system with a suitable numerical method [START_REF] Fornberg | A review of pseudospectral methods for solving partial differential equations[END_REF][START_REF] Hundsdorfer | Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations[END_REF].

Optimal design

For monitoring a population, the questions that we should address at first are where to look at the population and when to observe the targeted population? Though the time scheduling of the monitoring is often constrained by numerous factors, the choice of spatial locations generally offers more degree of freedom. In the following, we assume that the dates of the monitoring design (or experiment sensus lato) are scheduled in advance and also that the sensors locations will not change in time. This simplifies the problem of defining optimal monitoring designs and thus an optimal spatial distribution of observational locations over the survey area.

This question can be addressed through experimental design theory, a branch of applied statistics introduced by Sir Ronald Fisher in his seminal book edited in 1935 [START_REF] Fisher | The design of experiments[END_REF]. The purpose of this section is to introduce some key-notions of experimental design theory, essential for understanding its application on the use cases presented in section 4. For a complete description of the statistical framework the reader can refer to [START_REF] Silvey | Optimal design: an introduction to the theory for parameter estimation volume 1[END_REF], [START_REF] Atkinson | Optimum Experimental Designs[END_REF] or [START_REF] Walter | Identification of parametric models from experimental data[END_REF]. We start by introducing the key notions of optimal design and then we give their interpretation for the mechanistic-statistical model presented in equation 2. We indifferently use the terminology inherited from the design of experiments theory. Therefore, the experiments will refer to the survey of the populations and sensors corresponds to any monitoring device or a visual inspection.

Let us consider a random variable with a probability distribution depending on: (1) a vector of real variables that can be chosen by the experimenter, (2) and a vector of parameters , supposed to be fixed and unknown for the experimenter. We assume that belongs to the set  ⊂ R and that belongs to a parametric space Θ ⊂ R . The variables composing are called control variables. Let us suppose that for a given and the distribution of is given by a probability density function ( | , ). The experimenter is allowed to take independent observations on at vectors 1 , … , chosen from the set . The set = { 1 … } will be referred as -observation design. The question, primarily, is how to select the design ?

The criterion of choice depends on the purpose of the experiment. Here our primary interest is in estimating the parameter from the experimental data. Let us denote by a vector ( 1 , … , ) and let be a vector of values of taken at . The log-likelihood function of is defined as:

log Λ( ; , ) = ∑ =1 log( ( | , )) (3)
The maximum likelihood estimate of maximizes log Λ over Θ:

̂ ( , ) = arg max ∈Θ log Λ( ; , )
Under some regularity conditions on the family of densities { ∶ ∈ Θ} the estimate ̂ is asymptotically normal as the sample size tends to infinity (see [START_REF] Lehmann | Theory of Point Estimation[END_REF] for instance). Moreover, the asymptotic variance of ̂ reaches its lower bound given by Cramer-Rao inequality and is equal to the inverse of the Fisher Information Matrix, defined as:

( ; ) = ∑ =1 ( , ) = ∑ =1 Var log ( | , ) (4)
where ( , ) is the information matrix for an observation on taken at . Roughly speaking an optimal design ⋆ makes the variance of ̂ "as small as possible" or, alternatively, makes the Fisher matrix "as large as possible". To be more precise we seek for an ⋆ that maximizes some real-valued function of ( ; ). We assume here that is homogeneous and concave.

We consider the design with distinct vectors 1 , … , replicated 1 , … , times, where

∑ =1 = .
We assign to a discrete probability distribution which puts the probabilities = ∕ at . Consequently, the Information Matrix (4) can be rewritten as :

( ; ) = ( ; ) where :

( ; ) = ∑ =1 ( , ) (5)
As a criterion function is homogeneous, optimizing ( ( ; )) amounts to optimizing ( ( ; )) over . An optimal design ⋆ is thus defined as:

⋆ ( ) = arg max ∈Ξ ( ( ; )) (6)
where the matrix is assumed to be non-singular.

D-optimal design

The choice of the optimality criterion relies on the final purpose of statistical analysis. The experimenter may wish to improve the precision of ̂ or to reduce the variance of the predicted values of or again to better discriminate between the candidate models. The main optimality criteria and the corresponding functions can be found for example in [START_REF] Walter | Identification of parametric models from experimental data[END_REF]. Because of its versatility of purpose, one of the mostly used criterion is the D criterion leading to a D-optimal design, maximizing the logarithm of the determinant of the matrix :

( ( ; )) = log det ( ; ) (7)
A D-optimal design aims at minimizing the volume of confidence ellipsoid for model parameters. Consequently the statistical treatment of experimental results become more efficient. For example, it helps to assess if a model's parameter, and its linked hypothesis, has a non negligible effect on a population dynamics e.g. whether a parameter is different from zero. The D-optimal design also helps strengthen parameter comparison. For example it's usually important to estimate which habitat is the best (or worst) for a given species e.g. to significatively rank their corresponding growth rates. Finally, this design can also lead to minimizing the maximum variance of the predicted values [START_REF] Kiefer | General equivalence theory for optimum designs (approximate theory)[END_REF].

Exact and approximate design. When the set of the matrices defining a domain of ( 7) is discrete, the design is said to be exact.The problem of finding an exact optimal design ⋆ could be numerically challenging, especially for the large values of . The solution proposed in the theory of optimal design is to calculate the optimum of over the extended domain and to look for the discrete N-observational design which is "close" to the optimum. The definition of the matrix can be extended by considering the set of all probability distributions on  in place of the discrete distributions :

( ; ) = ∫  ( ; ) ( ) (8)
where is a probability distribution on . Then the approximate solution of the initial optimization problem can be calculated. An optimal (continous) approximate design ⋆ is defined as:

⋆ ( ) = arg max ∈Ξ ( ( ; ))
for the matrix given by the formula (8).

Average and local optimal design. An optimal design for a non-linear model depends on the unknown parameter value and is referred to as a local design. Different methods are proposed to "remove" the dependence from . One possible approach is to assume that is a random variable following a known distribution and to maximize the expected value of the criterion, calculated with respect to . Consequently, an on-average exact design maximizes:

E ( ( ; )) = ∫ Θ ( ( ; )) ( ) (9)
Fisher Information Matrix for generalized regression model. According to the definition given in [START_REF] Atkinson | Elemental information matrices and optimal experimental design for generalized regression models[END_REF], a generalized regression model is specified by the probability density that depends on and through the dimensional regression function ( , ):

( | , ) = ( | ( , ))
Consequently, the Fisher Information Matrix ( , ) for a single observation on can be rewritten as :

( , ) = Var log ( | ( , ) = D ( , ) Var log ( | ) ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏟ ( ) D ( , ) (10)
where D ( , ) is a × matrix of the first derivatives of the regression function with respect to and ( ) is a × Fisher Information Matrix for the reparametrized model. The matrix ( ) is called the elemental information matrix and is known for many useful probability densities.

Design settings for spatial population monitoring. We assume that the targeted population is observed at locations , … , . For each location, the measurement is repeated -times. The time measurements are fixed in advance and common for all the locations whereas the set of the locations is solution to the optimal design problem.

The measurement taken at a location at a time moment is denoted by ( , ). Typically ( , ) is either a population count or its suitable transform. In terms of mechanico-statistical model introduced in section 2, ( , ) is the result of a data process, a random variable following a probability distribution . The distribution  depends on a measurement process ( , ) governed by a partial differential equation and on a vector of parameters , related to all the levels of the hierachical model (see eq. 2). The likelihood of , given the vector of observations ( ( , )) =1,…, =1,…, is given by equation 3. The optimal design problem adressed in this paper is to find an exact D-optimal on-average design * = ( , … , )

(implicitly 1 = … = = 1). According to the short introduction given in this section, * satisfies :

* = arg max ∫  log det( ( , )) ( ) (11)
where ( , ) is calculated according to equation 5 and 10. From this point on, we assume that the parameter follows a uniform distribution  ( , ). In the following chapters, we illustrate the solution of optimal design * for two spatial monitoring examples.

Spatial population monitoring case studies

Agroecological case study

As a first example we consider the dynamics of Poecilus cupreus, a carabid beetle known for its weed seed consumption, within agricultural landscapes. This beneficial insect often complete its life cycle in a couple of months and is commonly monitored using pitfall traps that are picked up and replaced weekly. This species exhibits a single peak during its activity season and is known to react differently to semi-natural habitats, grasslands and cereal fields [START_REF] Marrec | Crop succession and habitat preferences drive the distribution and abundance of carabid beetles in an agricultural landscape[END_REF]. This leads to a reaction-diffusion model introduced in Parisey et al. ( 2016) with a birth rate decays parameter ( ), to mimick the single pick, and a spatially heterogeneous growth rate ( ), to express the habitat dependencies. In [START_REF] Parisey | Rearranging agricultural landscapes towards habitat quality optimisation: In silico application to pest regulation[END_REF], the population dynamics of carabids has been studied within a landscape of a few kilometer squared size. The dynamics of carabids was described by the following mechanistic-statistical model:

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ( , 0) = 0 ( ) ( , ) = Δ ( , ) + ( ( ) --) ( , ) ( , ) = ∫ - ( , ) ∼ ( ( , )) (12)
where ( , ) is in this case the density of carabids, describes the exponential speed at which the birth rate decays during the activity season, is the death rate (i.e.

1 life expectancy ), ( ) is an habitat specific growth rate such that :

( ) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ if x in a crop field if x in a semi-natural habitat if x in a grassland 0 otherwise (roads, buildings, . . . )
, is a scaling factor that allows the integration of the carabids density over the surface of a non-attractive pitfall trap from which count data is described by a Poissonian data process (.) with intensity ( , ).

The population process is parametrized by ( , ) with = ( , , , , ) and the measurement process depends on = ( ), that we considered known for a given type of traps. Hence the vector of unknown parameters were = ( , , , , ℎ , ). The parameters of the model were already estimated in [START_REF] Parisey | Rearranging agricultural landscapes towards habitat quality optimisation: In silico application to pest regulation[END_REF] for an arbitrarily chosen set of pitfall traps. We used these results in order to fix the bounds of parameters in an on-average design. and were set so they were biologicaly relevant and coherent with previous estimations. All growth rates were explored so they give half to twice as many descendants, during a season of several months, as previously estimated. The diffusion coefficient varies between a population of individuals half to twice as fast. As is linked with life expectancy, it varies between a life twice as long and half as long. Finally, , measuring birth decays during the activity season, varies between 0.9 and 1.1 times its nominal value. Details for all parameters values can be found in appendix A.1.

In this example, as described in [START_REF] Parisey | Rearranging agricultural landscapes towards habitat quality optimisation: In silico application to pest regulation[END_REF], the experiment consisted in placing = 24 pitfall traps in an agricultural landscape, owned by several farms, and sample them over the course of a season, in that case = 9

times over the course of three months in 2010. The sampling design of this study was decided beforehand so that, first, an agricultural field was chosen then several traps (here three) were planted, spatially grouped, in the field. From now on, we will refer to this as a clustered design.

Conservation biology case study

As a second example we considered the dynamics of Equus caballus, an invasive feral horse that has detrimental effects to the ecosystem of the Australian Alps and whose management through aerial culling or trapping and mustering is under debate. This question has been addressed by [START_REF] Beeton | Modelling horse management in the Australian Alps[END_REF] who used a two years aerial survey of feral horse populations, consisting in line transect data within the Australia Alps region [START_REF] Cairns | Feral Horses in the Australian Alps: the Design and Analysis of Surveys Conducted in April-May[END_REF]. They derived a reaction-diffusion model assuming horse populations are limited by density dependence in births via a logistic model, and that their movement through the landscape depends on the local horse density and, spatially heterogeneous, carrying capacity. We supplemented their spatio-temporal dynamics with a measurement process and a data process to obtain a full mechanistic-statistical model. The measurement process model the use of aerial transects to assess horse populations while the data process assumes normally distributed residuals with zero mean and variance 2 . This lead to a mechanistic-statistical system given by :

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ( , 0) = 0 ( ) ( , ) = Δ ( , ) ( ) + ( (1 - ( , ) ( ) ) -)) ( , ) ( , ) = ∫ ( , ) ∼  ( ( , ), 2 ) (13)
where ( , ) is in this case the density of horses, ( ) the habitat-dependent carrying capacity, the birth rate and the mortality rate. Population densities are monitored over subdomains ⊂ Ω, each subdomain corresponding to a horizontal transect with lower-left corner position , over which we integrate to obtain observations that are drawn from a Gaussian distribution  (.). We assume all transects have the same width (from aerial line of sight) but that the length will depend on each potential transect. The habitat-dependent carrying capacity was defined as :

( ) = ⎧ ⎪ ⎨ ⎪ ⎩ if (x) in a highly favorable habitat field, , otherwise.
The population process is parametrized by ( , ) with = ( , , ,

) and the data process depends on = ( 2 ). In this second case study the question we addressed was to find the best transects, over . We chose to simulate a survey of the population performed 7 years after the original monitoring (i.e. 2021).In [START_REF] Beeton | Modelling horse management in the Australian Alps[END_REF], the mortality rate was fixed and the carrying capacity ( ) was estimated independently of the population process. Hence, the vector of unknown parameters, for the conservation biology case, was = ( , ). One can note that is classically not part of such optimal design as it can be estimated by residual standard deviation. There was no previous estimations for but there were lower and upper boundaries, expressed as 'low growth and dispersal' and 'high growth and dispersal' scenarios in the original paper. We used them to set and used in the on-average exact design whereas for local design, we considered . Details for all parameters values can be found in appendix A.1.

In this example, as described in [START_REF] Beeton | Modelling horse management in the Australian Alps[END_REF], the experiment consisted in surveying = 32 aerial transects, of different lengthes, one time, in the year 2014, to estimate wild horses populations in the autralian alps [START_REF] Cairns | Feral Horses in the Australian Alps: the Design and Analysis of Surveys Conducted in April-May[END_REF].

Criteria for comparing designs and numerical implementation

Performance assessment of designs

In optimal design, one can quantify the suboptimality of any given designs compared to the D-optimal design, as described by eq. 5 and 7, using the notion of the D-efficiency [START_REF] Ucinski | Time-Optimal Path Planning of Moving Sensors for Par ameter Estimation of Distributed Systems[END_REF] which is defined as follows :

( ) = det( ( ; )) det( ( ; * )) 1 (14)
with the number of parameters in . We can use this measure to compare our designs among themselves, or with any other designs, e.g. taken from the literature or even randomly drawn. It can also be of interest to visualize the design properties. Especially, one can review design positions on maps. One might also want to focus on the relationship between pair of parameters, using confidence ellipses [START_REF] Murdoch | A graphical display of large correlation matrices[END_REF] derived from the inverse of the Fisher matrix.

Evaluating designs according to a practical constraint

When planning an ecological monitoring, the statistical properties of the design is only part of the decision making.

One important factor for planning the monitoring in the two study cases considered here is the length of travelling that induces important costs. In order to integrate this component we consider the tour length of travelling salesman solutions [START_REF] Rosenkrantz | An Analysis of Several Heuristics for the Traveling Sa lesman Problem[END_REF] where one start at a relevant 'base camp', then tour a design before going back to camp. Formally, as defined in [START_REF] Dantzig | Linear Programming and Extensions[END_REF], we label the traps from 1, … , and define :

= ⎧ ⎪ ⎨ ⎪ ⎩
1 the path goes from trap to trap 0 otherwise Taking > 0 to be the Euclidean distance from trap to trap , the travelling salesman problem (TSP) can be written as the following integer linear programming problem:

min ∑ =1 ∑ ≠ , =1 ∶ (15) subject to ∶ ∑ =1, ≠ = 1 = 1, … , ; ∑ =1, ≠ = 1 = 1, … , ; ∑ ∈ ∑ ≠ , ∈ ≤ | | -1 ∀ ⊊ {1, … , }, | | ≥ 2
As one can see, the TSP minimize the tour length while ensuring each trap is arrived at from exactly one other trap, departed to exactly one other trap and that the solution returned is a single tour and not the union of smaller tours. For the second use case, involving lines instead of points, a heuristic is used instead of the problem 15. The tour length is roughly approximated by a weighted mean between the diagonal of the bounding box of the transects and their cumulated lengthes. Finally, the best compromises between statistically efficient designs and travelling lengths can be visualized on a Pareto front in the D-efficiency -tour length space [START_REF] Sheftel | The geometry of the pareto front in biological phenotype space[END_REF]. Of course, we can note that performance assessment of an experimental design is not limited to the above examples.

Numerical implementation

We implemented different strategies to solve numerically the computationally demanding problems we worked on.

As the computation of the on-average design (eq. 11) with a uniform distribution would be cumbersome, we relied on a sampling from a Latin Hypercube, a space-filing design [START_REF] Pronzato | Design of computer experiments: space filling and beyond[END_REF], to 'cover' at best a domain with a limited number of sampled points.

Moreover, we noted that one useful feature that come at no cost with our designs is the ability to filter the set of possible points (or transects) before searching for a solution. For example, in the agroecological case, we used a regular grid to select possible trap locations within fields ( >> ) while excluding road from pitfall trap's possible positions.

In addition, in order to solve the information matrix that requires the derivatives of the deterministic process models we chose to numerically extract them ahead of solving the design, i.e. to do some memoization [START_REF] Michie | Memo" Functions and Machine Learning[END_REF], which proved relevant for the task at hand. For performance, one could also rely on analytic, symbolic or automatic differentiation (adjoint code) if applicable.

Finally, we relied on an excursion algorithm [START_REF] Fedorov | Theory Of Optimal Experiments[END_REF], with multiple starts, for solving the designs. These algorithms take the particular form of the problem into account. Starting from a non-degenerated design (i.e. with positive determinant), at each iteration , exchange one support point by a better one, in the sense of the D-optimality.

Considering all * ( -) possible exchanges successively, retaining the best one among these each iteration, it converges to a local solution. Hence, we used multiple starts to get the best local solution as an approximation of a global solution.

Numerical calculations were performed in R (R Core Team, 2019). Using a Runge-Kutta solver combined with the method of lines to solve the reaction-diffusion models [START_REF] Soetaert | Solving differential equations in r: Package desolve[END_REF] and some naive parallelization. It took about an hour to solve an on-average design on a Intel ® Xeon ® E5 with 8 cores 8 Gb of RAM.

Results

We reduced the size of the original landscapes, and thus the number of sampling points, to reduce the computational cost. For the agroecological case study, we used a virtual agroecosystem, of one km 2 , inspired by [START_REF] Soubeyrand | Parameter estimation for reaction-diffusion models of biological invasions[END_REF] with just six pitfal traps sampled over a season. For the conservation biology case study, we used a simplified habitat map, of several thousand km 2 , of a subpart of the australian alps (see supplementary materials of [START_REF] Beeton | Modelling horse management in the Australian Alps[END_REF]), sampled one time, 7 years after the first experiment (hence in 2021). For this example, we reduced the number of surveyed transects to eight, chosen among the original thirty two, through D-optimality.

For each use case, we compared empirical designs (e.g. random and/or clustered) to local and on-average D-optimal ones, assuming they both outperform the empirical designs. Theoretically, locally optimal designs should be the most efficient ones regarding the information gain on the processes. However, they require a good knowledge on parameters before the survey which is, in practice, rarely the case. Thus, given the usual uncertainty on populations traits and considering on-average designs appears more relevant.

Agroecological case study

In most studies investigating insect populations in agricultural landscapes, traps are placed according to what we call a clustered design. In practice, when a field is chosen, three spatially grouped traps are put in this field as illustrated in Fig. 1b. For this case, we compared some properties of local D-optimal, on-average D-optimal, clustered and random As shown in Figure 1, the initial conditions were generated as a gradient from bottom left to top right, and the landscape contains 2 crop fields (orange), two grasslands (yellow), one region of semi-natural habitat (e.g. forest patch in green) and two crossing roads (gray).

The local D-optimal design displayed in Figure 1d concentrates the points on the bottom left corner of the landscape where the highest slope of initial conditions is. Moreover, it seems to cover all the habitats, even if unequally, except roads where the growth rate is set to 0. Figure 2 shows the position of the traps in an on-average design and 95% confidence ellipses, estimated by inverting the Fisher information matrix, for the pairs of growth rates ( , ), ( , ) and ( , ). Interestingly, compared to the local-design, the stencil shape seems less apparent here as we see four points align within the initial condition density but two points are farther away and distributed between two different habitats. As expected, the confidence ellipses surface are smaller for all pairs of parameters for the on-average D-optimal design compared to the cluster one (Fig. 2c-d &Appendix A.2 for those related to , and ). The uncertainty associated with semi-natural growth rate is severely reduced which is logical given the clustered design did not account for that habitat. In fact, the optimal design suggests that only one point within this habitat drastically reduces the variance of its estimates. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1.5 2.0 2.5 3.0 0.00 0.25 0.50 0.75 1.00 Efficiency Tour Length (a) Locally D-optimal q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1.5 2.0 The comparison of the efficiency of the designs in relation to the tour length is presented in figure 3 for both local and on-average designs. Random designs are represented by black circles (there are a hundred of them), the clustered design by a blue dot, the D-optimal design by a green dot and if at least one red dot exist, it is a random sampling that is pareto optimal for at least one dimension. One can see that the clustered design is, of course, less statistically efficient than the D-optimal but also that it has a longer tour. In fact, the D-optimal designs here are dominant both in statistical and tour length. Several simulations of random designs show that it is not always the case but the optimal designs are still among the shortest (data not shown). One can also note that most random designs are more D-efficient than the clustered design. This is probably because they either place traps within semi-natural habitat and/or simply cover more ground, with a longer tour length. As in this study we drew only one clustered empirical design, the solid comparison between random and clustered sampling strategies, which was not the scope of this work, would require more simulations.

Conservation biology case study

In this example the design problem consists in searching = 8 aerial transects, of different lengths, chosen among N=32 reported in [START_REF] Cairns | Feral Horses in the Australian Alps: the Design and Analysis of Surveys Conducted in April-May[END_REF]. The expected values and the range of parameters are given in Appendix A.1. As shown in figure 4a the monitored area has less high quality habitats (green) than average and low quality ones (grey). The original transect coverage was extremely dense (Fig. 4b) and it appears natural to reduce it, for instance by taking only 1 4 of the possible transects. As short transects are frequent in the original empirical distribution of transects, there are a few chances of randomly choosing the longer ones (Fig. 4c). When defining a local D-optimal design with a 'low growth and dispersal' scenario, the longer and more informative transects appear to be favored (Fig. 4d). (c) Confidence ellipse q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 3.0 The on-average D-optimal design (5b) appears to be similar to the local design (Fig. 4d). It seems to favor longer transects and, in that case, only one transect on the southern part of the map. The spatial extent of the optimal transects is smaller than most random transects but, as it favors longer transects, the tour length, forcibly passing by transects, is actually longer than most (Fig. 5d). As expected, D-optimality reduces the confidence ellipse of estimated parameters compared to a random design (Fig. 5c). The variance reduction appears to be greater for the diffusion coefficient than for the growth rate . Finally, the Pareto front displayed in Figure 5d to assess the compromise between design efficiency and tour length, shows that the on-average D-optimal transects dominate the D-efficiency but is here the worst choice in term of tour length. This means that there are multiple point on the Pareto front that could be candidate if a compromise has to be chosen between D-efficiency and tour length.

Discussion

In this study we addressed the question of where to make observations within a large spatial domain to efficiently capture the changes in space and time of a population described by a reaction-diffusion model. By considering two example systems, i.e. the dynamics of carabid beetles in agricultural landscapes and wild horses in the Australian Alps, we showed that the use of the optimal design of experiments framework can be useful to find optimal locations for trapping or counting individuals that maximize the information on the population for parameter estimation. In particular, we were able to produce spatial designs that outperform both those used in the initial studies and random designs. More precisely we first obtained local D-optimal designs that assume that the expected values of parameters are known before planning data collection. Then, as this situation is generally unlikely to occur in the population ecology, we investigated the situations when only bounds around population parameters are known considering on-average D-optimal designs. While local designs appeared to be obviously the most efficient for minimizing the variance of parameters estimates, our results points out the usefulness of optimal design without strong a priori to think spatial monitoring so it will enable efficient model fitting.

Although the question of optimized survey in ecology or epidemiology is still seldom considered, it has already been addressed in several studies using various frameworks. For instance, [START_REF] Williams | Monitoring dynamic spatio-temporal ecological processes optimal ly[END_REF] proposed a framework that produces designs that minimizes prediction uncertainty, [START_REF] Cook | Optimal observation times in experimental epidemic processes[END_REF] calculated optimal observation times for botanical epidemiology experiments by maximizing the Kullback-Leibler divergence from the posterior to the prior, or [START_REF] Bourhis | Artificial neural networks for monitoring network optimisation-a practical example using a national insect survey[END_REF] who proposed a strategy for choosing insect traps locations by minimizing uncertainty in neural networks outputs. Here, we considered the optimal design of experiments framework that is briefly introduced in part 3. While its use is common for solving experimental design issues in environmental sciences or systems biology [START_REF] Steiert | Experimental design for parameter estimation of gene regulatory networks[END_REF], it is poorly considered in ecology. The monitoring of population on large spatial domain is actually similar to the problem of optimal sensors placement for spatially distributed systems [START_REF] Ucinski | Optimal measurement methods for distributed parameter system identification[END_REF].

Despite the recent development of devices for increasing ecological monitoring, in practice the survey of population in large spatial domains remains challenging and a major bottleneck to the study of how populations change in space and time and the design of efficient management strategies [START_REF] Nichols | Monitoring for conservation[END_REF][START_REF] Lindenmayer | The science and application of ecological monitoring[END_REF].

As stated above, statistical frameworks developed to optimize how and where one should make observations can be used to design more efficient population surveys. This model-based design of experiments can also take into account several practical constraints. In our study we chose to consider the tour length between traps for collecting carabid bettles and aerial transects from which wild horses are enumerated. Then we used a Pareto front to asses how designs efficiency changes with the tour length. As illustrated in Fig. 5d, it is possible to find good candidate that maximizes design efficiency while minimizing tour length so one can choose those that fit the compromises that are acceptable. For the agroecology example, we found that the best design was also minimizing the tour length. This may be explained by the smaller spatial extent of the domain compared to the wild horse case. In fact, it is likely that increasing the spatial domain would create a more progressive Pareto front between design efficiency and tour length. Depending on the considered problem, other constraints could be included in the approach to find designs that are statistically efficient and minimized the constraints. For instance one could imagine several travellers and seek for the strategy that minimizes the cost of the survey, considering both human and transport costs. While, in our opinion, finding monitoring designs that provide maximum information on population spread is an important question on which the community of modellers in ecology should concentrate more, integrating realistic constraints in model-based design is likely to catch the interest of field ecologists and practitioners. Moreover, in our two study cases we focused on the spatial aspect of monitoring assuming that the only degree of freedom is the survey designs was the locations of observations whereas the time of observations was fixed. This assumption could be relaxed and both temporal and spatio-temporal D-optimal designs may be obtained in further studies. Optimal times at which observations should be done may be obtained by solving local or on-average D-optimality problems before the survey starts. In addition, if applicable, the problem can also be tackled sequentially as an alternative to the on-average design. In this case, optimal groups of observations are proposed, over time, after having performed a parameter estimation on already acquired data, this approach is know as a sequential design [START_REF] Chernoff | Sequential design of experiments[END_REF].

Apart the theoretical background, the main obstacle to a widespread usage of optimal design of experiments in spatially explicit systems is the computational cost. Using usual algorithms and computing strategies (i.e. memoization and parallelization), in this study we were able to solve local and an-average D-optimality problems on reactiondiffusion models with a reasonable amount of time. However, increasing the complexity of the problem, for instance by adding more practical or economic constraints or the temporal dimension of sampling, may require more advanced computing methods and grid computing. Structural and practical identifiability issues may also arise and make the problem even more complex. Optimal design is useful for maximizing the information gain from data collection, however the modeller has to ensure that enough information is available so the parameters can be determined. In this study, before seeking for an optimal design, we assessed the practical identifiability of the model by multiple start of a numerical local approach [START_REF] Walter | Identification of parametric models from experimental data[END_REF] within each relevant domain Θ.

The issue of initial conditions is well known by mathematicians and modellers working with dynamic models. In some conditions, it is possible to estimate initial conditions from monitoring data. Nevertheless, in most studies one has to make strong hypothesis and find a way to fix it before performing a statistical inference, or seeking an optimal design. Regarding spatial ecology, the most frequent approach is perhaps to use past (or early data) with a regression model to someway estimate the distribution of population density at the beginning of the study. In our study, we did not include the initial conditions in the designs problems. It is very likely that optimal designs are dependent on the initial conditions. This may be assessed by augmenting the problem and exploring how changes in the initial conditions impact optimal-designs. Yet, augmenting the problem for taking care of initial conditions has its limits.

To cite in extenso [START_REF] Uciński | Optimum experimental design for infinite dimensional inverse problems[END_REF], it is "[a] very prospective direction [where] the infinite dimensional nature of the resulting parameter space is inherently associated with the ill-posedness which means that even low noise in the data may make the estimates extremely unstable.". The reader could refer to [START_REF] Alexanderian | A-Optimal Design of Experiments for Infinite-Dimension al Bayesian Linear Inverse Problems with Regularized 0 -Sparsification[END_REF] and [START_REF] Alexanderian | A fast and scalable method for A-optimal design of experiment s for infinitedimensional Bayesian nonlinear inverse problems[END_REF] for further information.

In this article we focus on the D criterion but there exist different criteria (e.g. different functions (.)) that can be applied to the information matrix. For example, one can seek to minimize the average variance of the estimates by minimizing the negative of the trace of the inverse of the Fisher information matrix (i.e. A-optimality). Here, we focused on the D criterion because (i) through the equivalence theorem [START_REF] Kiefer | General equivalence theory for optimum designs (approximate theory)[END_REF], we know that main classical criteria are linked in some ways (e.g. and ), (ii) the D criterion is invariant by reparameterization (e.g. helpful to control parameters positivity) and (iii) there exist variants of the D criterion that are useful for a variety of purposes.

For instance, if one is interested in model selection, the criterion enable to focus the problem on the parameters of interest versus 'nuisance' parameters, hence being a canonical criterion in nested model selection [START_REF] Atkinson | Planning Experiments for Discriminating Between Models[END_REF]. For more general purpose model selection in case of gaussian observations, the generalized DT criterion is a criterion that can be described as a normalized arithmetic mean of criterions and 2 distance between predictions, handling both optimal parameters estimation and model selection [START_REF] Atkinson | DT-optimum designs for model discrimination and parameter estimation[END_REF].

One can also note that, quite often, the aim of ecological surveys is not only to estimate reaction-diffusion model parameters but rather to monitor the state of the population and provide abundances. For such goal, it would be better to use a criterion more directly related to predictions. For instance, the G criterion seeks to minimize the maximum entry in the diagonal of the hat matrix, minimizing the maximum variance of the predicted values i.e. ̂ ( , ). As [START_REF] Pázman | Optimum design accounting for the global nonlinear behavior of the model[END_REF] refined the G criterion for nonlinear homoscedastic or heteroscedastic regression models, we think an on-average G-optimal design would be more suitable than D-optimal designs to focus on population state monitoring.

Even if knowing model's parameters and knowing is related, it corresponds to different mathematical objectives (e.g.

G versus D optimal designs) that may produce distinct designs. One can note that some limitations might have to be worked through for the G-optimal design in case of the generalized regression model (i.e. with any kind of data model) and the numerical aspects of the problem.

To finish with, the purpose of this article was to push the discussion on the statistical efficiency of experimental designs linked with spatial ecological models, as question that is yet seldom expressed this way. The main goal of a A. Appendix

A.1. Use case's parameters

Table 1 Parameters used for the use cases when searching for local and on-average optimal designs. The symbol ′ -′ means 'not applicable' e.g. when used for explored min or max, it means those parameters were not part of the optimal design problem so they were kept at their original values. Reference [1] is [START_REF] Parisey | Rearranging agricultural landscapes towards habitat quality optimisation: In silico application to pest regulation[END_REF] and reference [2] is [START_REF] Beeton | Modelling horse management in the Australian Alps[END_REF].

Parameters

A.2. Supplementary confidence ellipses

Figure 1 :

 1 Figure 1: Pitfal trap placements comparisons : (a) initial population densities; (b) clustered design with 3 traps per field (pseudo-replicate) in different fields; (c) randomly positioned pitfal traps ; (d) local D-optimal design

Figure 2 :

 2 Figure 2: On-average exact D-optimal design : (a) position of pitfal traps in case of on-average D-optimal exact design ; (b-d) 95% confidence ellipses for pair of growth rates, estimated by inverting the Fisher Information Matrix for , for the on-average optimal design (red) and the clustered design (blue) seen in fig. 1.

Figure 3 :

 3 Figure 3: Linear approximation of the Pareto front : (a) for the local design; (b) for the on-average design. Each point represents the D-efficiency (x) and tour length (y) of a design. Black circles represent 100 random samplings (the same set is projected on each panel), blue dots represents the clustered sampling, green dots are optimal and red lines and dots form the linear approximation to the pareto fronts. All the designs D-efficiencies are evaluated locally (a) or on-average (b).

Figure 4 :

 4 Figure 4: Transects comparisons: In the original study, high initial densities were linked with a high carrying capacity while mid and low densities were linked with another carrying capacity (see also Appendix A.1). We can see here (a) the initial population densities (b) then all 32 original transects (in black) over the carrying capacities ; (c) some random subset of 8 of the transects (in blue) ; (d) finally a locally D-optimal design of 8 transects for (i.e. 'low growth and dispersal' scenario) (in red)

Figure 5 :

 5 Figure 5: On-average exact D-optimal design : (a) random transects; (b) on-average D-optimal ; (c) 95% confidence ellipse of b and D, blue is for the random design figured in (a) and red for the optimal design, both evaluated at and, finally, (d) a linear approximation of the Pareto front with D-efficiency evaluated at (color coding and shapes are the same as in figure 3)

Figure 5 :

 5 Figure 5: Confidence ellipses at 95% for the agroecological case study, where parameters , and are paired with the crop growth rate.

Acknowledgements.

The authors thank Jean-Sébastien Pierre, Jacques Baudry and Solenn Stoeckel for useful discussions about this work. The authors want to thank especially Professor Darius Uciński for a brief but very supportive discussion about an early version of this work.

Software and data availability

All our results can be reproduced, extended and adapted to other use cases by cloning the code and data available at https://github.com/nparisey-pro/LANDoE under GPL license. The software was developed using R (R Core Team, 2019) with a dozen package dependencies. Packages list, versions and tested operating systems are listed at the git repository.

field ecologist is to sample 'as much as possible' within the constraints of its budgets (monetary and manpower), but this strategy can either fail to capture the essential information on the population or be very expensive. In addition to the numerous studies that had demonstrated that formalizing ecological processes into mathematical models and using them to analyze empirical data offers a mean for improving our understanding of populations spread, we hope this work points out that the difference between optimal and non-optimal monitoring can be significant in term of information gained, and thus that model-based design of ecological survey is a promising path that could also contribute to reduce environmental costs of samples collection, storage and processing.

CRediT authorship contribution statement

Nicolas Parisey: Implemented the framework, ran the numerical experiments and provided focus on optimal design. Melen Leclerc: provided insight on ecological modeling. Katarzyna Adamczyk-Chauvat: provided fruitfull guidance on statistics.All authors designed the numerical experiments and the overall framework, drafted, managed, reviewed and approved the manuscript.

Competing interests.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Funding.

The authors thank the ANR project "Clonix 2D" (ANR-18-CE32-0001), and its leader Solenn Stoeckel, for its financial support of this work.