Antoine Brias
email: antoine.brias@irstea.fr

Jean-Denis Mathias

Guillaume Deuant

A Brias Irstea

Accelerating viability kernel computation with CUDA architecture : application to bycatch shery management

Keywords: Viability kernel, Dynamic programming, CUDA, GPU, Fishery management Mathematics Subject Classication (2000) 90B50 -Management Decision Making, including Multiple Objectives, 90C39 -Dynamic Programming

Accelerating viability kernel computation with CUDA architecture: application to bycatch fishery management

Introduction

Viability theory provides mathematical and numerical tools and concepts for maintaining a dynamical system within a set of states (called the constraint set). This theory has numerous potential applications in food processing,nance, economics, environment [START_REF] Sicard | A viability approach to control food processes: Application to a Camembert cheese ripening process[END_REF][START_REF] Béné | A viability analysis for a bio-economic model[END_REF][START_REF] Doyen | A stochastic viability approach to ecosystembased sheries management[END_REF][START_REF] Bernard | Comparing the sustainability of dierent action policy possibilities: application to the issue of both household survival and forest preservation in the corridor of Fianarantsoa[END_REF][START_REF] Andrés-Domenech | Sustainability of the Dry Forest in Androy: A Viability Analysis[END_REF][START_REF] Chapel | Dening yield policies in a viability approach[END_REF][START_REF] Chapel | Viability and resilience of languages in competition[END_REF][START_REF] Mathias | Using the Viability Theory to Assess the Flexibility of Forest Managers Under Ecological Intensication[END_REF]. A essential step in this approach lies in the computation of the viability kernel. This is the set of states for which there exists a control policy that keeps the system in the constraint set for some (nite or innite) time. Since the 90's, several algorithms have been developed to compute viability kernels in dierent application elds [START_REF] Saint-Pierre | Approximation of the Viability Kernel[END_REF][START_REF] Bokanowski | An anti-diusive scheme for viability problems[END_REF]. Some of them are adapted to nite horizon problems [START_REF] Djeridane | Approximate Viability using Quasi-Random Samples and a Neural Network Classier[END_REF], linear dynamics [START_REF] Kaynama | A Modied Riccati Transformation for Decentralized Computation of the Viability Kernel Under LTI Dynamics[END_REF] or particular kinds of problems such as single output nonlinear control systems ane in the control [START_REF] Turri | A method to construct viability kernels for nonlinear control systems[END_REF][START_REF] Mattioli | Noyau de viabilité : une contrainte globale pour la modélisation de systèmes dynamiques[END_REF]. The concepts have been extended to the stochastic case [START_REF] Doyen | Stochastic viability and dynamic programming[END_REF][START_REF] Rougé | Extending the viability theory framework of resilience to uncertain dynamics, and application to lake eutrophication[END_REF], in which dynamics includes random variables accounting for uncertainties. However, most of them are limited by their computational complexity. When the number of dimensions increases, the required computation time and memory increase exponentially because they are mainly based on a grid covering the state space. This curse of dimensionality is a moot point in a number of domains, like data mining, numerical analysis [START_REF] Donoho | High-Dimensional Data Analysis: The Curses and Blessings of Dimensionality[END_REF] and multiple approaches aim at mitigating it and at speeding up viability kernel computing. For instance, [START_REF] Bonneuil | Computing the viability kernel in large state dimension[END_REF] attempts to handle large dimensional state space with computing control policies on a given time horizon with simulated annealing. [START_REF] Deuant | Approximating Viability Kernels With Support Vector Machines[END_REF] use Support Vector Machines for approximating viability kernels. [START_REF] Maidens | Lagrangian methods for approximating the viability kernel in high-dimensional systems[END_REF] try to overcome this curse by using Lagrangian methods, which do not call for a grid. It oers particularly interesting performances for linear dynamics. Designing new algorithms is a crucial way to reduce the curse of dimensionality, but existing technological solutions can also lead to important gains of time. Some programs like Vikaasa, used in by-cash shery management, embed a multi-core implementation of viability kernel computation [START_REF] Krawczyk | Computation of viability kernels: a case study of by-catch sheries[END_REF]. However, in recent years, another technical improvement arose which seems well adapted for the viability kernel problematic : the GPU parallel computing technology. Historically, GPUs are electronic components used in computer graphics hardware. They are getting more widely used to accelerate computations in many elds [START_REF] Mametjanov | Autotuning Stencil-Based Computations on GPUs[END_REF][START_REF] Cekmez | Adapting the GA Approach to Solve Traveling Salesman Problems on CUDA Architecture[END_REF][START_REF] Sabo | Fast evaluation of a time-domain non-linear cochlear model on GPUs[END_REF]. Due to its architecture, the GPUs solution is advantageous for problems for which a parallel execution is possible [START_REF] Goldsworthy | A GPUCUDA based direct simulation Monte Carlo algorithm for real gas ows[END_REF]. It often exceeds the capacity of computing clusters in terms of performance and for a lower cost. Thus, NVIDA has created the CUDA platform to oer tools and to allow programmers to use the parallel architecture of graphic cards. The basic point is to separate the program in two parts : the serial code executed by the host (the Central Processing Unit (CPU)), and the parallel code, executed in the GPU threads across multiple processing elements.

In this paper, we use the standard viability algorithm [START_REF] Saint-Pierre | Approximation of the Viability Kernel[END_REF] which can be considered as a particular case of dynamic programming. To illustrate our work, the implemented algorithm is tested on a basic multidimensional model of population dynamics. In particular, we compare the computation time between the parallelized version and the sequential one. Finally, the algorithm is used to compute the viability kernel in a case of bycatch shery management.

2 Viability theory and dynamic programming

Viability kernel

We consider a discrete dynamics in which h is the function mapping the state and the control at time t with the state at the next time step t + 1:

x (t + 1) = h (x (t) , u (t)) (1)
with x (t) ∈ X the set of states, u (t) ∈ U (x (t)) which is a nite set of possible controls allowing the regulation of the dynamics and h (., .) : X × U → X , associating a state x (t) ∈ X and a control u (t) ∈ U (x (t)) with a successor

x (t + 1) ∈ X . For example, it associates the closest grid point of the successor obtained with an Eulerian scheme dened with an appropriate local parameter.

We want to assess the viability kernel V iab T (K); the set of initial states for which a control strategy exists to maintain the system inside K, the set of desirable states (the viability set), until T (a nite horizon). This is rewritten in Eq.2 :

V iab T (K) = {x 0 ∈ K|∃ (u (0) , u (1) , ..., u (T -1)) , ∀t ∈ {0, 1, ..., T } , x (t) ∈ K} (2)
with x (t) the state of the system at date t and u (t) the chosen control value at t, which maintains the system in the viability set K.

Dierent algorithms are available for computing this set. Because of the approximation of the dynamical system on a grid, these algorithms are submitted to the dimensionality curse: the necessary time and memory increase exponentially with the number of dimensions of the problem. We choose to focus on a dynamic programming algorithm which is very standard and well suited for a parallel approach.

Dynamic programming

We assume that the dynamics of the system fulls the condition ensuring the convergence of the approximation to the actual kernel, when the resolution of the discretization tends to 0 (Saint-Pierre, 1990). The dynamic programming algorithm is a backward recursive algorithm and determines a value function : V (t, x) at date t. The seed value at T is :

V (T, x) =        1 if x ∈ K 0 if x / ∈ K (3)
In this case, the value function directly determines the viability kernel. At T , being viable is the property of belonging to the viability kernel. The linear recurrence relation is given by :

∀t ∈ {0, T -1} , V (t, x) =        max u∈U (x) V (t + 1, h (x, u)) if x ∈ K 0 if x / ∈ K (4)
Thus, at t = 0:

V (0, x) =        1 if x ∈ V iab T (K) 0 if x / ∈ V iab T (K) (5)
The principles of dynamic programming are respected because, at each step time :

-each point of the grid is a subproblem (V has to be computed for each point);

-a value function has to be optimized independently for these points;

The value function V takes only two values : 0 and 1. We are looking for initial states, for which the value function at time horizon T is 1. Eq.4

is equivalent to say that the value function V (t, x) is equal to 1 if a control u ∈ U (x) exists, maintaining the system in the constraint set K. Otherwise, the value function V (t, x) equals to 0. The backward dynamic programming equation denes the value function V (x, t). By writing a max instead of a sup, we implicitly assume the existence of an optimal solution for each time t and state x. The algorithm involves two main steps:

For each point of the grid, we store all the successors (the point of the grid at the next time step) obtained for each choice of control (the possible set of actions at each point is supposed nite).

The second part performs iterations of kernel approximations by progressively excluding points of the grid until reaching a xed point or a predened number of time steps.

Thus we have the algorithm Alg.1 as follows :

Algorithm 1 Detailed Dynamic Programming Viability Kernel Algorithm

Storing successors :

         i = 0 K 0 = K store coordinates, ∀x ∈ K store h (x, u) , ∀x ∈ K, ∀u ∈U (x)
Iterating kernel approximation :

         repeat K i+1 = x ∈ K i |V (i + 1, x) = 1 i = i + 1 until K i+1 = K i OR K i+1 = Ø OR i = T V iab T (K) = K i

Parallelizing dynamic programing algorithm

Dynamic programming is well suited to parallelization in the two main parts.

During the storing of successors, all succcessors x (t + 1) of each state and each value of u, dened in Eq.1, can be computed separately. In the second part of the algorithm, it is possible to parallelize the update of the value function V (t, x) applied to each state (Eq.4). In both cases, computations are completely independent for one step of time. The parallelizing of all these calculations thus improves the performance in terms of computation time.

Since GPUs have hundreds of processor cores running in parallel, they can be more ecient than using traditional CPU architectures.

3 Using GPU for computing viability kernels 3.1 GPU architecture As they are widely used in many scientic computing elds, graphics hardware programming for general purpose computing is well documented, this includes ways to write and debug the code [START_REF] Langdon | Debugging CUDA[END_REF]. Providing highly eective data parallelism, GPU cores contain multiple threads running at the same time which can accelerate viability algorithms. The CUDA programming framework enables the use of all the cores belonging to the GPU for general purpose computing. Programs written in this paper use the C++ CUDA framework. A GPU computing view consists in splitting the input data between the dierent threads, and performing operations on these pieces of data, independently and thus achieving parallelization. This kind of computing encourages to rethink programs to make this parallelization possible. Fig. 1 illustrates the architecture of a GPU device which is seen as a grid housing blocks. Each block involves a number of threads performing the computations. The GPU device owns a global memory, each block has a memory shared by all of its threads, and each thread owns a smaller local memory. Here, we use the global memory to store big arrays of data with which all the threads can interact. In a block, a thread is indexed over 3 dimensions (a block can be shown as a cube of threads), thus the indexation has to be carefully done. Each thread and each block are indexed at their level. In our implementation, we want to get the ID threadID of each one-dimensional thread in a grid of 3D blocks. According to the CUDA framework writing convention , this ID is given by :

threadID =threadIdx.i + blockDim.i * blockIdx.i+ gridDim.i * blockIdx.j * blockDim.i+ (6) gridDim.i * gridDim.j * blockIdx.k * blockDim.i
Here, threadIdx.i means the coordinates i of threadIdx. gridDim contains the dimensions of the grid (the grid owns gridDim.i×gridDim.j ×gridDim.k blocks), blockIdx contains the block indexes within the grid (in 3 dimensions : blockIdx.i, blockIdx.j and blockIdx.k), blockDim contains the dimensions of the block (a block contains blockDim.i × blockDim.j × blockDim.k threads in 3 dimensions) and threadIdx contains the thread indexes within the block (in 3 dimensions too, but for simplicity, only one is used here : threadIdx.i).

These variables are mainly used in assigning the work per block and thread.

The maximum values in each dimension in both blockDim and gridDim are GPU-dependent. At least, the i component needs to be declared. By default, the value for the j and k components is 1.

Overall strategy

The organization of the algorithm is summed up in Fig. 2. The gure shows the storing of trajectories and the iterating kernel approximation steps of Alg.1.

The steps are :

-the initialization and setting of the GPU device (a); -the copy of results from the GPU memory to the CPU (f) using the CUDA instruction cudaM emcpy with the ag cudaM emcpyDeviceT oHost;

-The clear of the memory of the GPU and the nalization of the GPU device (g), using the CUDA instruction cudaF ree;

-the disallocation of the CPU memory (h).

Parallelizing the computing and the storing of successors

The input of the algorithm is the number of points by dimension and each dimension bound. The algorithm computes the coordinates of all the points and their index. It can be parallelized since each point is independent from the other. Thus, one thread on the GPU will compute the coordinates of one point.

Now, we have to compute the all the successors of each point of the grid, when applying the dierent values of the control (if the control can take continuous values then it is discretized). We can thus suppose that the control takes nbControls values point of the grid. Then, on each thread, the coordinates of one state X i are sent, and all the nbControls possible evolutions from this state are returned. Again, the number of threads is equal to the number of points of the grid. This treatment is summed up on Fig. 3. There are four steps in the GPU computation of the trajectories :

-the inputs are the values of the discretized control, U = (u 1 , ..., u nbControls )

and the coordinates of all the points X = (X 1 , ..., X nbP oints ). These coordinates were determined beforehand and remain in the GPU global memory;

-this data is sent on threads, each thread holding the computation of all the controlled evolutions of only one point;

-each thread then sends back the result in the array H. This array contains the index of the controlled evolutions. Moreover, it contains -1 if the evolution doesn't belong to K. Thus, we store the index of i-th controlled evolution of X in H (X, i). As the number of threads is the same for each block, some threads can be superuous on the last block;

-once we get the controlled evolutions, we just have to copy them from the GPU to the CPU and we apply the second part of the dynamic programming algorithm.

Parallelizing the iterating kernel approximation

The second part of the algorithm computes the value function V for each point until the stopping condition is reached. As shown on Fig. 4

) If ∀i ∈ [1, nbControls -1] , (H (X, i) = -1) OR (H (X, i) ≥ 0 AND v (H (X, i)) = 0) vT mp (X) = 0 CheckChange = CheckChange + 1 End If End If
If the point is viable and if a viable successor exists, nothing is done, the point is still viable. But if none of its successors is viable, i.e. if none of its successors belongs to K (i.e., H (X, i) = -1) or if the value function V of the successors is equal to 0, then the point is considered non-viable : 0 is put in the appropriate cell of the array vT mp. This treatment represents the value function updating in dynamic programming. In this case, a change is made, then CheckChange is incremented. At the end of the step, the value function is transferred of vT mp to v and we check CheckChange. If it is equal to 0, the algorithm can stop. It happens after a sucient number of iterations and we then get the approximation of the viability kernel.

Results

First, we test the GPU programming on a simple example of population dynamics, in order to compare the performance with the sequential algorithm.

Then, we use our parallel algorithm to treat a problem of bycatch shery management up to 6 dimensions.

A multi-dimensional population model

One toy model is used here to evaluate the benets of GPU programming. We consider the population model studied by Malthus and Verhulst. Malthusian ideas are the basis of the classical theory of population and growth [START_REF] Ehrlich | The problem of population and growth: a review of the literature from Malthus to contemporary models of endogenous population and endogenous growth[END_REF]. This model is a simple dynamical system of population growth on a bounded space without any predator. The system state includes two variables : the size of the population and its evolution rate. [START_REF] Aubin | An Introduction to Viability Theory and Management of Renewable Resources[END_REF] introduced a control on this evolution rate and adapted it to the management of renewable resources. The control is applied on the derivative of the evolution rate at each time step with the inertia bound c. The equations ruling the system in discrete time are : The program is implemented in C++, using CUDA 5. The computer CPU is an Intel R Xeon(R) CPU E5-2687W 0 @ 3.10 GHz x 16 and the GPU card is a Tesla K20C with a computing capability of 3.5 and 2496 CUDA cores.

x i (t + 1) = x i (t) + x i (t) y (t) , i ∈ {1, .., n} (7)
y (t + 1) = y (t) + u (t) (8)

Speed-up

Here we show a comparison between the two versions of the program applied to this model. The speed-up S is a measure evaluating the performance of the parallelization:

S = T seq T par (9)
where T seq is the sequential algorithm execution time and T par the execution time of the parallel one. The rst one uses the usual algorithm; the second one uses the adapted algorithm to GPU computing.

We gure the evolution of the speed-up, in function of the size of the grid.

In Fig. 5, we study the inuence of grid renement in the case of a 4-dimension problem.

First, according to Fig. 5a, we can see that the speed-up exceeds 5 in most cases of 4-dimension problems. It rises, using the whole GPU capacity but, then, the curve stagnates. The speed-up reaches a threshold when the GPU is fully loaded and where no more time gain is possible. On Fig. 5b, a comparison of elapsed times during the CPU and GPU version of the algorithm is shown.

We notice that when the renement of the grid is high, the computation time increases signicantly in the CPU version.

Fig. 6, shows speed-up isolines in function of the number of dimensions and the number of points by dimension. We observe that when we have a 2-or 3-dimension problem, the speed-up is between 0 and 4. This is not signicantly time-saving, and the initialization of the GPU can be slower than the classic array lling. However, when the number of dimensions grows, the speed-up increases up until 20 in 7 dimensions or more in 5 or 6 dimensions with a ner discretization. The speed-up isolevels follow the total number of points of the grid for high-dimension problems (4 dimensions or more). In these problems, the main part of the time is spent in the trajectories storing part (the iterating kernel approximation ends quickly because of the low accuracy of the grid).

Parallelized parts

Here, we compare the speed-up during the storing of successors and during the computation one between the CPU and the GPU versions of the algorithm. In the following gure (Fig. 7), we study two dierent sizes of problems (2 and 6 dimensions) and we return speed-ups in function of the number of points by dimension. Elapsed times during the two parallelized parts are also returned.

On Fig. 7a, we see that the elapsed time during the iterating kernel approximation is now higher than the one for storing the successors for 2-dimensions problems (for problems with more than 3000 points per dimension for the CPU version, or 1600 points per dimension for the GPU version). For problems in 3 dimensions or more (like 6 dimensions in Fig7c), storing successors time is still higher than iterating kernel approximation time. The time for storing successors exponentially increases with the number of dimensions and parallelizing this part of the algorithm is very ecient.

The speed-ups of the storing of trajectories and iterating kernel approximation are calculated separately here. As shown on Fig. 7b and Fig. 7d the speed-up increases up to 27 in 2 dimensions and 34 in 6 dimensions for the storing of trajectories, and up to 15 in 2 dimensions and 3.7 in 6 dimensions for the iterating kernel approximation.

For a small problem with a low discretization, the time taken by the two algorithms is negligible, but we observe that the speed-up is below one (up to 400 points by dimension in the 2D-problem, and 11 points by dimension in the 6D-problem). The mandatory initialization of the CUDA device explains the less successful results of the GPU version in smaller problems.

Moreover, we notice that the parallelization of the storing of successors is more ecient whatever the number of dimensions. Actually, memory exchanges (reading the value function array in the global memory from each thread at each step) in the second part alter the performance of the parallelization.

Thus, the speed-up of the iterating kernel approximation falls down to 3.7 in 6 dimensions. Since the time elapsed during this step is small compared to the time elapsed during the storing of trajectories, its impact is not signicant on the global speed-up.

Application to a bycatch shery

We illustrate our results on a multidimensional bycatch shery model. This model is an extension of the one studied by [START_REF] Krawczyk | Computation of viability kernels: a case study of by-catch sheries[END_REF] from 3 dimensions to 6 dimensions. We study the inuence of the number of bycatch species on the shery management.

The two-species dynamics

The model of Krawczyk is reminded here. There are two ecologically independent populations of sh x and y, harvested by a single eet. The shery focuses on harvesting the target stock x. We dene a catchability coecient q x determining the quantity of biomass that each unit of eort extracts, relative to the total size of the biomass at the time. Thus, the harvest rate at time t is:

h x (t) = q x e (t) x (t) (10)
The product of q x and e (t) is generally termed the level of shing mortality and it is expressed as a proportion of biomass.

A linear bycatch production function is used :

h y (t) = αh x (t) (11)
where the constant α, 0 < α < 1 is the bycatch to target harvest ratio.

The respective population dynamics are then described by :

x (t + 1) = x (t) + r x x (t) 1 - x (t) L x -q x e (t) x (t) (12)
y (t + 1) = y (t) + r y y (t) 1 - y (t) L y -αq x e (t) x (t) (13)
where r x , r y , L x and L y are all positive constants. By convention, r x and r y represent the intrinsic growth rates, and L x and L y are the environment carrying capacities, of x and y, respectively. According to Krawczyk, the numerical values for the bycatch species are chosen to have a bycatch stock less productive and less valuable than the target stock.

We suppose that the eort e can be controlled:

e (t + 1) = e (t) + u (t) , u (t) ∈ [u min , u max ]
The speed at which the regulator can change shing intensity is bounded by u min and u max . The shing eet's prot is given by :

π xy (t) = p x h x (t) + p y h y (t) -ce (t) -C (14)
where p y represents the unit prices of y, p x the unit price of x The marginal cost of the shing eort is c and C is a xed cost. We dene the viability constraints set K :

K = (x, y, e, u) : x (t) ≥ L x 10 , y (t) ≥ L y 10 , π xy (t) ≥ 0, e (t) ∈ [e min , e max ] (15)
The safe minimum biomass levels are set to one tenth of each sh population's unexploited level, as commonly implemented in sheries worldwide [START_REF] Krawczyk | Computation of viability kernels: a case study of by-catch sheries[END_REF].

The n-species dynamics

We propose to extend this model to treat n-species dynamics in order to illustrate the benets of the GPU computation of the viability kernel. The dynamics of the target species remains unchanged (Eq.2). We suppose here that we have n bycatch species y i . Moreover, we assume that these species are at the same trophic level, competitively sharing the same food resource [START_REF] Ekerhovd | Optimization in the `Pelagic Complex': A Multi-Species Competition Model of North East Atlantic Fisheries[END_REF][START_REF] Rice | Does functional redundancy stabilized sh communities[END_REF]. The bycatch species growth is limited by the common carrying capacity L y . The discrete dynamics becomes:

y i (t + 1) = y i (t) + r yi y i (t)     1 - n-1 k=1 y k (t) L y     -α i q x e (t) x (t) , ∀i ∈ [1, n -1] (16)
r yi is the intrinsic growth rate. α i denotes the bycatch to target harvest ratio, a measure of how highly coupled the production relationship is. Each bycatch stock y i is a by-product of the production process. A linear bycatch production function is used :

h yi (t) = α i h x (t) (17)
We assume here that the sum of all bycatch harvest ratios α i is equal to α dened in the previous section :

n-1 i=1 α i (t) = α (18)
For reasons of simplicity, we suppose that all the α i are equal :

α i = α n (19)
As before, economic sustainability requires that the activity remains profitable. The shing eet's prot is given by :

π xy (t) = p x h x (t) + n-1 i=1 p yi h yi (t) -ce (t) -C (20)
Then, we have the viability constraint set K :

K = (x, y 1 , ..., y n-1 , e, u) : x (t) ≥ L x 10 , ∀i ∈ [1, n -1] , y i (t) ≥ L y 10n , π xy (t) ≥ 0, e (t) ∈ [e min , e max ] (21)
In the following, we study dynamics with 1, 2, 3 or 4 bycatch species, which means dynamics with 3, 4, 5 or 6 dimensions respectively (including the shing eort and the target species biomass).

Results

We use the parameters calibrated by Krawczyk [START_REF] Krawczyk | Computation of viability kernels: a case study of by-catch sheries[END_REF] and summed up in the Table 1.

With only one bycatch species, a similar shape of the viability kernel is obtained than the one in the previous study [START_REF] Krawczyk | Computation of viability kernels: a case study of by-catch sheries[END_REF]. This viability kernel is presented in Fig. 8. For the following gures, the biomasses are scaled between 0.1 and 1.

The bounds of the gure are the constraint domain K. From each state contained inside the viability kernel, it exists at least a control strategy keeping the shery in K. The states outside of the viability kernel are not viable, there is no control strategy satisfying the bio-economic constraints. For these states, the shery is facing a "crisis" situation. Now, bycatch species are added to the environment. For one bycatch species, 41.9% of states are inside the viability kernel. It decreases until less than 7%

for four bycatch species.

The Fig. 9 shows 3D slices of the 4D kernel for dierent shing eort values. It represents the biomass of the target species and the biomass of the bycatch species required to be viable to satisfy the ecological constraints while maintaining the protability of the shery. A small shing eort needs big available stock of the target biomass in order to conserve the system profitability. A mean shing eort leads to a large slice of the viability kernel according to the target biomass. Finally, keeping a big shing eort requires a mean target species stock to not deplete the bycatch populations. The results are displayed for 4D problems, but they can be extended to higher dimension problems.

Finally, we show on Fig. 10 the number of states belonging to the viability kernel for a 6D problem with a grid of 20 points by dimension. The biggest viability kernel is obtained when the starting population of bycatch species are equal to

Ly n (the oset is set by the grid step which is nearly 15.3 here). The environment is then completely lled with balanced increasing the possibility of sustainability. For larger initial populations, the environment is overloaded, causing the depletion of bycatch species. When the initial populations are smaller, the ecological sustainability requirement restricts the number of viable states.

Discussion and conclusion

The GPU parallelization provides signicantly faster viability kernel computations and tackles problems with a higher number of dimensions or more precise discretization than a sequential algorithm.

However, a huge storage space is required to save all the coordinates and all the successors for high-dimension dynamical problems. For example, in C++, storing one coordinate in double-precision oating-point format required 8 bytes (64 bits). Then, storing the coordinates of 100 millions points in a 8dimension problem (10 points per dimension) requires 6.4Gb. Moreover the successors multiplied this value by the number of possible values of the control.

Because of this size, it is impossible to send all the data at once. Instead, we cut the problem in smaller parts : computing the coordinates for some grid points and then computing the successors g (x, u) for these points, the result is then saved. The process then continues with other grid points until all the whole grid has been processed. This reduces the need for memory.

This parallelized version may be improved with other tools (like STXXL in C++ to use large arrays) and using the mapped memory, which is based on some pointers to the RAM. Further improvements could be achieved by associating multi-CPUs and multi-GPUs architectures, the rst ones cutting the problem in small parts and the second ones solving the subproblems.

The parallelization approach proposed here opens up perspectives in terms of parallelization for other viability problems. Then, revisiting other viability approaches such as an extension to stochastic dynamics [START_REF] Doyen | Stochastic viability and dynamic programming[END_REF] or algorithms such as support vector machine algorithm [START_REF] Deuant | Approximating Viability Kernels With Support Vector Machines[END_REF] could also lead to signicant improvements of eciency. Finally, the parallelization of the viability kernel algorithm provides a tool for the management of socio-ecosystems as illustrated with the management of the multi-species bycatch shery. v ( H (X threadID , u i ))

-

  the initialization of the model parameters (b). The user inputs these parameters; -the initialization of CUDA parameters (number of blocks, number of threads by blocks in each dimension)(c); -the creation of empty arrays, which will contain controlled evolutions and the value function of each point using the CUDA instruction cudaM alloc (d); -the copy of parameters from host to the GPU device (e) using the CUDA instruction cudaM emcpy with the ag cudaM emcpyHostT oDevice; -the computation of successors of each point i.e. the storing of next point of the trajectory for each possible value of the control Alg.1. This operation is parallelized on the GPU and it is done in the rst kernel function. Section 3.3 detailed this part. -the iterating kernel approximation step of Alg.1. The value function of each point is calculated and ran until the stop condition is reached. It is paral-lelized in the second kernel function on the GPU. Section 3.4 detailed this part.

  4.1.1 The model As we want to analyze the performances of the parallelized algorithm on multidimensional problems, we articially increase dimensions of the system in adding some independent variables in order to study the behaviour of the program in dierent dimensions. The new model is as follows : we considered n states x i representing the size of n populations, which all have the same evolution rate y (t) ∈ [d, e]. The size of the populations must remain in K = [a, b].

  with -c ≤ u (t) ≤ +c. The set K = [a, b] n × [d, e] is the viability constraint set. For the rest of this work, bounds for K are set as a = 0.2, b = 3, c = 0.5, d = -2, e = 2 . By default, the control is discretized in 5 values. Since the parallelization simply breaks down processes that are already independent in the usual algorithm, the discretization does not impact the accuracy of the results.

Fig. 1

 1 Fig.1GPU device architecture. The device is composed of a grid of threads blocks. Each level posseses its own memory. The parallelized tasks, like computation of controlled evolutions and computations of the value function, are made separately on the threads. A grid of 3x2x2 blocks of 4x2x1 threads is shown.

Fig. 3

 3 Fig. 3 Storing of trajectories : the parallelized controlled evolution computation from coordinates. The dierent steps of the parallelization are shown. The computation of all the controlled evolution for one point is done on one thread.

Fig. 4

 4 Fig. 4 Iterating kernel approximation: the parallelized viability computation. The dierent steps of the parallelization are shown.

Fig. 6

 6 Fig. 5 Eect of grid renement on the speed-up in a 4-dimensions problem. The elapsed time during the two versions of the algorithm are shown in (b).

  Fig. 7 Comparison of the storing of trajectories and the iterating kernel approximation elapsed times and speed-up between the GPU and CPU versions of the algorithm. Two sizes of problems are shown : 2-dimension problem (a) and (b), and 6-dimension problem (c) and (d).

Fig. 8

 8 Fig. 8 Viability kernel for the 3D shery model (with one bycatch species).Viability kernel for the 3D shery model (with one bycatch species). The blue surface symbolizes the zeroprot surface. Having a positive prot (being above the blue surface) is a property of the viability constraint set K.

  Fig. 9 3D slices of the 4D kernel for dierent initial values of the shing eort e.

Table 1

 1 Parameters of the bycatch shery. The parameters for all the bycatch species i are the same.

	Name	Variable Value
	Target carrying capacity	L x	600
	Target growth rate	r x	0.4
	Target unit price	p x	4
	Target catchability coecient	q x	0.5
	Bycatch common carrying capacity	L y	300
	Bycatch growth rate	r yi	0.2
	Bycatch unit price	p yi	1.9
	Bycatch harvest ratio	α yi	0.2
	Marginal cost	c	10
	Fixed cost	C	150
	Maximum eort	e max	1
	Minimum eort	e min	0.1
	Maximum eort variation	u max	0.01
	Minimum eort variation	u min	-0.01

Acknowledgments

This work was supported by grants from Irstea and Région Auvergne. This support is gratefully acknowledged.


Franck Jabot

Théophile Lohier

Non-random correlation of species dynamics in tropical tree communities

The importance of environmental stochasticity for tropical tree dynamics has been recently stressed by several studies. This has spurred the development of a "time-averaged neutral model" of community dynamics by Kalyuzhny and colleagues that extends the neutral model by incorporating environmental stochasticity. We here show that this framework can be used to assess the presence of non-random correlations between species dynamics. Indeed, the time-averaged neutral model makes the simplifying assumption that species responses to environmental variation are uncorrelated. We therefore propose to use this model as a null hypothesis against which observed community dynamics can be compared. This study makes five contributions. First, we describe a novel time-averaged neutral model of community dynamics that is close to, but more flexible than the one previously proposed by Kalyuzhny and colleagues. Second, we develop an inference method based on approximate Bayesian computation (ABC) and demonstrate the identifiability of the model parameters from community time series data. Third, we develop a test of the significance of environmental stochasticity, and a method to quantify its contribution to population variance. Fourth, we develop a test of non-random correlation between species dynamics.

Fifth, we apply these developments to three datasets of tropical tree dynamics. We evidence both a strong contribution of environmental stochasticity to population variance in the three datasets, and a non-random correlation of species dynamics in one of them. We finally discuss the implications of these results for the modelling of tropical tree community dynamics.

Introduction

The study of the assembly and dynamics of tropical tree communities has been largely influenced in the past fifteen years by neutral theory [START_REF] Hubbell | The unified neutral theory of biodiversity and biogeography[END_REF][START_REF] Rosindell | The unified neutral theory of biodiversity and biogeography at age ten[END_REF].

This theory considers as pivotal the role of demographic stochasticity in community dynamics, that is the random birth and death of tree individuals irrespective of the species they belong to. This theory has been shown to reproduce many community patterns including species abundance distributions [START_REF] Volkov | Neutral theory and relative species abundance in ecology[END_REF], beta-diversity [START_REF] Condit | Beta-diversity in tropical forest trees[END_REF] and phylogenetic imbalance [START_REF] Jabot | Inferring the parameters of the neutral theory of biodiversity using phylogenetic information and implications for tropical forests[END_REF]. In contrast, when looking at community dynamics in community time series, several patterns cannot be reproduced by neutral models. In particular, the scaling of population size variance with respect to population size has been shown to be steeper in real communities than in neutral ones [START_REF] Leigh | The average lifetime of a population in a varying environment[END_REF][START_REF] Chisholm | Temporal variability of forest communities: empirical estimates of population change in 4000 tree species[END_REF][START_REF] Kalyuzhny | A neutral theory with environmental stochasticity explains static and dynamic properties of ecological communities[END_REF].

The fact that populations of abundant species have a larger temporal variability than expected under demographic stochasticity indicates that the dynamics of conspecific individuals is correlated. One straightforward way to explain this correlated dynamics of conspecific individuals is to consider that they respond in the same manner to environmental temporal variability. This process has been named environmental stochasticity in the literature [START_REF] Lande | Stochastic population dynamics in ecology and conservation[END_REF]. [START_REF] Kalyuzhny | A neutral theory with environmental stochasticity explains static and dynamic properties of ecological communities[END_REF] made a critical contribution to the modelling of tropical tree community dynamics by showing that a community model incorporating both demographic and environmental stochasticity was able to reproduce a number of both static and dynamical patterns of the tropical tree community in Barro Colorado Island, including species abundance distribution, the scaling of population variance with population size, population fluctuation magnitude, and community temporal turn-over [START_REF] Kalyuzhny | A neutral theory with environmental stochasticity explains static and dynamic properties of ecological communities[END_REF]. They called this model "the time-averaged neutral model", since each individual has the same over-time-average prospect of birth and death irrespective of the species it

Author-produced version of the article published in Oikos, 2016, 125 (12), 1733-1742.

The original publication is available at http://www.oikosjournal.org/ doi : 10.1111/oik.03103 belongs to, although in contrast with a standard neutral model, individual birth rates differ among species at each time step.

In order to keep the model simple, they used a rather strong simplifying assumption: that species individual responses to environmental variability were uncorrelated (see methods).

This simplifying assumption may appear doubtful. For instance, tree species differ in their sensitivity to drought [START_REF] Engelbrecht | Drought sensitivity shapes species distribution patterns in tropical forests[END_REF], so that dry periods should benefit similarly to drought resistant species. And the same reasoning applies to any environmental component of species fitness [START_REF] Hutchinson | Concluding remarks. -Cold spring harbor symposium on quantitative biology[END_REF]. Still, an « effective » independence across species could result from the multi-dimensionality of environmental variation if species that respond similarly on one environmental axis do not co-vary on other axes. Recent analyses on tropical tree functional traits substantiate this possibility [START_REF] Baraloto | Decoupled leaf and stem economics in rain forest trees[END_REF]).

In addition, [START_REF] Kalyuzhny | A neutral theory with environmental stochasticity explains static and dynamic properties of ecological communities[END_REF] also considered in their model that there was a temporal autocorrelation in environmental conditions (modelled through a within-species temporal autocorrelation in fecundity rates), making species fitness temporally autocorrelated. This within-species temporal autocorrelation indirectly produces some level of between-species correlation in their dynamics. It could therefore contribute to alleviating the potential model error due to the assumption of between-species independence.

The aim of this study is thus to apply a rigorous statistical framework to analyze tree community dynamics with environmental stochasticity, and to use this framework to assess both the contribution of environmental stochasticity to population variance and the presence of non-random correlations between species dynamics. In particular, we assess whether the model parameters can be reasonably well estimated from community time series data, notably the variance A and the temporal autocorrelation ρ of environmental stochasticity. We then develop a simple test of whether environmental stochasticity and within-species temporal autocorrelation are at play in a community using a model selection approach, and quantify their contributions to population variance using virtual communities.

Author-produced version of the article published in Oikos, 2016, 125 (12), 1733-1742. The original publication is available at http://www.oikosjournal.org/ doi : 10.1111/oik.03103

Finally, we assess the presence of non-random correlations between species dynamics.

These developments are applied to three datasets of tropical tree dynamics [START_REF] Condit | The importance of demographic niches to tree diversity[END_REF].

Material and methods

Modelling community dynamics with demographic and environmental stochasticity

We consider a model very similar to the one of [START_REF] Kalyuzhny | A neutral theory with environmental stochasticity explains static and dynamic properties of ecological communities[END_REF] with two main differences. First, birth, immigration and death events are happening independently of each other at specified rates, so that there is no zero-sum rule and the community size can vary through time. Second, the temporal autocorrelation in environmental conditions is modelled for each species as a yearly draw from an autocorrelated lognormal distribution, rather than by independent draws every τ years as in [START_REF] Kalyuzhny | A neutral theory with environmental stochasticity explains static and dynamic properties of ecological communities[END_REF]. These two modifications do not change the spirit of the model of [START_REF] Kalyuzhny | A neutral theory with environmental stochasticity explains static and dynamic properties of ecological communities[END_REF]. They were made to increase the generality of the model and to make it more flexible to accommodate observed temporal variations in community size (Table 1).

More precisely, in the model considered, during year t, each individual of species i present in the community 1) dies with a constant probability d independent of the species considered, and 2) produces descendants at a rate (1 -m) × d × f i (t), where m is the expected proportion of immigrants in the recruits, and the fecondity f i (t) of species i during year t is drawn from a temporally autocorrelated lognormal distribution with mean 1, variance A, and withinspecies autocorrelation ρ on a log scale (i.e., cor[ ln(f i (t+1)) , ln(f i (t)) ] = ρ). Concretely, the number of descendants produced by an individual of species i during year t is drawn from a

Poisson distribution with mean equal to (1 -m) × d × f i (t). This model of local community dynamics is thus neutral in the sense that every individual has the same average over time prospect of birth and death [START_REF] Chave | Neutral theory and community ecology[END_REF]), although birth rates differ among individuals at each time step. This is why [START_REF] Kalyuzhny | A neutral theory with environmental stochasticity explains static and dynamic properties of ecological communities[END_REF] proposed to call this type of models, timeaveraged neutral models, a terminology that we adopt here. The choice of the lognormal distribution to model among year and species variabilities in birth rates is a standard choice

Author-produced version of the article published in Oikos, 2016, 125 (12), 1733-1742. The original publication is available at http://www.oikosjournal.org/ doi : 10.1111/oik.03103

for modelling environmental stochasticity [START_REF] Lande | Stochastic population dynamics in ecology and conservation[END_REF]. It is supported by previous analyses of tree demographic rates in tropical forests [START_REF] Condit | The importance of demographic niches to tree diversity[END_REF]. We nonetheless compared this model choice to an alternative one based on a Weibull distribution that was shown to better fit bird dynamical properties in North America (Kalyuzhny et al. 2014b) and showed that the lognormal distribution provided a better fit to tropical tree data (Table S6).

Importantly, alternative choices for the modelling of environmental stochasticity could be easily plugged into the approach detailed in this study.

On top of this local dynamics, immigration into the local community occurs at rate m × d, so that the total number of immigrating individuals during year t is drawn from a Poisson distribution with mean equal to J(t) × m × d, where J(t) is the community size at time t. These immigrants are drawn from a neutral regional pool with regional diversity θ, as in the neutral model of [START_REF] Hubbell | The unified neutral theory of biodiversity and biogeography[END_REF]. Here, for computing efficiency, we model the species abundance distribution in the regional pool with a logseries distribution.

We additionally considered two degraded versions of this model: a model without temporal autocorrelation, in which ρ is equal to 0, and a model without environmental stochasticity in which A and ρ are equal to 0. We also considered an augmented version of the model incorporating a trend in community size. In this augmented version, we additionally consider that the community as a whole has a carrying capacity K, so that death rates during year t are no longer considered constant and equal to d, but they are equal to d × J(t)/K. All the model versions considered are thus time-averaged neutral, except the one that does not include environmental stochasticity and is therefore neutral.

A commented R code to simulate this model of community dynamics with demographic and environmental stochasticity is provided in the supplementary material (Appendix S7).

Tropical tree data

Author-produced version of the article published in Oikos, 2016, 125 (12), 1733-1742. The original publication is available at http://www.oikosjournal.org/ doi : 10.1111/oik.03103

We make use of three datasets of tropical tree communities gathered by the Center for Tropical Forest Science, and published in [START_REF] Condit | The importance of demographic niches to tree diversity[END_REF]. The three 50 ha forest plots are located in Barro Colorado Island (BCI) in Panama, in Mudumalai in India, and in Pasoh in Peninsular Malaysia (Table 1). These three tropical tree communities have been censused more than three times, which is a minimum number of censuses to estimate both the variance and the temporal autocorrelation of environmental variability (BCI: 5 censuses, Mudumalai and Pasoh: 4 censuses, see Table 1). In these three forests, we will use all the tree individuals above 10 cm of diameter at breast height, which is a common arbitrary size threshold to distinguish adult trees from saplings [START_REF] Condit | The importance of demographic niches to tree diversity[END_REF].

Inference by ABC and a test of the presence of environmental stochasticity and temporal autocorrelation

We developed an inference method for these models based on Approximate Bayesian Computation (ABC, [START_REF] Beaumont | Approximate Bayesian computation in population genetics[END_REF]. ABC contains four main steps: 1) simulating a large number of times the model that is under study, with parameters drawn from prior distributions, 2) computing for each model simulation a set of summary statistics, 3) comparing these simulated summary statistics to the same set of statistics computed on the data, and 4) retaining the best fit simulations if their associated summary statistics are closer to the data summary statistics than an arbitrarily defined tolerance threshold. The parameter values used for these simulations form the posterior distribution of the parameters. A more detailed description of the ABC procedure for the present case study is given below.

Preliminary analyses revealed that successful inference could be performed with a set of 6 summary statistics: the average number of death events between two successive censuses N death that is correlated with model parameter d (Fig. S1a), the average community size J that is correlated with K (Fig. S1f), the average species richness S across censuses that is correlated with m (Fig. S1b), the total number of species across censuses S tot that is correlated with θ (Fig. S1c), the average across species initially present of their standard

Author-produced version of the article published in Oikos, 2016, 125 (12), 1733-1742. The original publication is available at http://www.oikosjournal.org/ doi : 10.1111/oik.03103 deviation of population sizes SD pop that is correlated with A (Fig. S1d), and the average across species initially present of their temporal autocorrelation coefficients R(1) pop that is correlated with ρ (Fig. S1e). This choice of summary statistics was largely inspired by the study of [START_REF] Kalyuzhny | A neutral theory with environmental stochasticity explains static and dynamic properties of ecological communities[END_REF] and by previous works on the inference by ABC of neutral and non-neutral models [START_REF] Jabot | Analyzing tropical forest tree species abundance distributions using a nonneutral model and through approximate Bayesian inference[END_REF]. More precisely, some parameters jointly influence several summary statistics, notably m and θ that both influence S ans S tot ., and A and ρ that both influence SD pop and R(1) pop . This will create correlations in the posterior distributions of these parameters, as we will see below.

We first performed a model selection through ABC, to assess whether all model parameters could be inferred, or whether there was some level of over-parameterization. Distinct uniform prior distributions were used for parameter K in the three datasets: [10000 ; 30000] for BCI, [1000 ; 15000] for Mudumalai, and [25000 ;50000] for Pasoh. These prior distributions were chosen to encompass observed summary statistics in the three datasets.

For each simulation, the model parameters were randomly drawn from the corresponding prior distributions. The community was initialized using the first census, and the community dynamics was simulated starting from this initial condition for the number of years separating the first from the last census of the studied dataset (20 years for the BCI dataset, 12 for Mudumalai and 13 for Pasoh). Note that this initialization procedure implies that parameter inference will be based solely on the community trajectory from this first census, as is classically done in community time-series analyses (e.g., [START_REF] Ives | Estimating community stability and ecological interactions from timeseries data[END_REF], and not on the likelihood of obtaining this initial community composition based on the model parameters evaluated. At the end of each simulation, the six summary statistics mentioned above were computed. They were normalized by their median absolute deviations across simulations [START_REF] Csilléry | abc: an R package for approximate Bayesian computation (ABC)[END_REF]. And the distance between the simulation and the data was

Author-produced version of the article published in Oikos, 2016, 125 (12), 1733-1742. The original publication is available at http://www.oikosjournal.org/ doi : 10.1111/oik.03103 computed as the Euclidean distance in the 6-dimensional space of the summary statistics between the simulated and observed normalized summary statistics [START_REF] Beaumont | Approximate Bayesian computation in population genetics[END_REF].

Based on these four sets of ABC simulations for the four model versions, model selection was performed using the method of [START_REF] Beaumont | Joint determination of topology, divergence time, and immigration in population trees[END_REF] with a tolerance rate of 0.01 (meaning that the best 0.01 × 80000 = 800 ABC simulations were retained), implemented in the R package "abc" [START_REF] Csilléry | abc: an R package for approximate Bayesian computation (ABC)[END_REF]. This method fits a multinomial logistic regression on the summary statistics based on the set of 800 retained simulations, and evaluates the four model posterior probabilities as the multinomial logistic prediction at the value of the summary statistics observed in the data. The idea of using model selection to assess the number of parameters that can be efficiently inferred within an ABC procedure was initially proposed in [START_REF] Lagarrigues | Approximate Bayesian computation to recalibrate individualbased models with population data: Illustration with a forest simulation model[END_REF]. Here, this model selection approach is also a way to test for the presence of environmental stochasticity and within-species temporal autocorrelation: their presence will be supported if the models without these components are rejected.

We also assessed the quality of the fit of the selected models with a model checking procedure [START_REF] Rubin | Bayesianly justifiable and relevant frequency calculations for the applied statistician[END_REF][START_REF] Csilléry | Approximate Bayesian computation (ABC) in practice[END_REF]. This consists in assessing a posteriori that the summary statistics of the best-fit simulations are effectively close to the values observed in the data. For each summary statistics, we computed a p-value of model-data mismatch by computing the proportion of the 200 best-fit simulations that is equal or more extreme that the observed value. Since there is no a priori expectation of the direction of a potential misfit, this p-value must be lower than 0.025 to indicate a significant deviation at the 0.05 level. However, a correction for multiple testing must be used. We used here the Holm-Bonferroni correction.

Assessing the contribution of environmental stochasticity to population variance

We quantified the respective influences of demographic and environmental stochasticities on observed mean population variances in the three datasets. For the BCI dataset, for which there was no trend in community size in the selected model (see results), we proceeded as

Author-produced version of the article published in Oikos, 2016, 125 (12), 1733-1742. The original publication is available at http://www.oikosjournal.org/ doi : 10.1111/oik.03103

follows: we simulated community dynamics using model parameter values from the posterior distributions except that environmental variance A was fixed at 0. We then compared, for three species abundance classes (1-19, 20-199, >200), the median population variances VAR pop_demo obtained without environmental stochasticity to the observed median population variances VAR pop_obs . The contribution of demographic stochasticity to observed population variances was quantified as VAR pop_demo / VAR pop_obs and the contribution of environmental stochasticity was quantified as 1 -VAR pop_demo / VAR pop_obs .

For the two other datasets for which there was a trend in community size in the selected model, we proceeded slightly differently. Indeed in these datasets, the trend in community size also contributes to the overall population variances. We also simulated community dynamics using model parameter values from the posterior distributions, but we performed two sets of simulations: one with environmental variance and the second without environmental variance (A fixed at 0). In these two sets of simulations, the fitted trend in community size was removed, using death rates independent of community size, with values drawn from the posterior distribution of parameter d. The contribution of demographic stochasticity to population variance was quantified as above as VAR pop_demo / VAR pop_obs while the contribution of environmental stochasticity to population variance was quantified as (VAR pop_env -VAR pop_demo )/ VAR pop_obs , where VAR pop_env is the median population variance obtained in the simulation set with both demographic and environmental variances. The remaining part 1-(VAR pop_env / VAR pop_obs ) quantifies the contribution of the trend in community size to population variance.

Assessing non-random correlations between species dynamics

Due to the short length of the time series considered (5 censuses in BCI, 4 in the two other datasets), direct tests of between-species correlations in their dynamics would have a low power. Besides, the within-species temporal autocorrelation incorporated in two of the four models considered can potentially generate some level of between-species correlations in their dynamics. Consequently, one needs to take this model component into account when

Author-produced version of the article published in Oikos, 2016, 125 (12), 1733-1742. The original publication is available at http://www.oikosjournal.org/ doi : 10.1111/oik.03103

assessing the assumption of species independence in their responses to environmental conditions. We thus developed a test of the presence of non-random correlations between species dynamics, using a model checking procedure [START_REF] Rubin | Bayesianly justifiable and relevant frequency calculations for the applied statistician[END_REF][START_REF] Csilléry | Approximate Bayesian computation (ABC) in practice[END_REF]. We proceeded in two steps. First, a more precise parameter inference with the selected model was performed using the sequential ABC algorithm of [START_REF] Lenormand | Adaptive approximate Bayesian computation for complex models[END_REF] within the R package "EasyABC" [START_REF] Jabot | EasyABC: performing efficient approximate Bayesian computation sampling schemes using R[END_REF]. This procedure required 14,500 ABC simulations for the BCI dataset, 14,000 for Mudumalai, and 19,000 for Pasoh. Second, we sampled parameters from the obtained posterior distribution and computed, for each simulation below the tolerance threshold obtained by the sequential algorithm, the correlation matrix of the most abundant species with an initial abundance of at least 50 individuals. This correlation matrix was computed as follows: for two species i and j, the term (i,j) in the matrix was equal to the correlation coefficient Cor ij between the two trajectories of population sizes across censuses of species i and j. 200 such correlation matrices were computed and they form the posterior predictive checks that can be compared to the correlation matrix computed from the data [START_REF] Csilléry | Approximate Bayesian computation (ABC) in practice[END_REF]. Rather than comparing each correlation coefficient one by one, we compared a unique index: the average of the absolute values of the matrix coefficients |Cor pop |.

We finally replicated the model checking analysis for the BCI dataset, using either the first three censuses or the first four censuses, to assess whether considering shorter time series might change our understanding of the correlations between species trajectories.

Results

Parameter identifiability and a test for the presence of environmental stochasticity and within-species temporal autocorrelation

Models with temporally autocorrelated environmental stochasticity were selected against the simpler nested models in the three datasets (Table 2), and their fits of the summary statistics were overall satisfactory (Table S2). This indicates that the model studied is not over-parameterized and that both environmental stochasticity and within-species temporal

Author-produced version of the article published in Oikos, 2016, 125 (12), 1733-1742. The original publication is available at http://www.oikosjournal.org/ doi : 10.1111/oik.03103 autocorrelation are critical to explain the observed dynamics in the three tropical forests considered. In two of the three datasets (Mudumalai and Pasoh), a model further incorporating a trend in community size was selected. This result reminds that the standard zero-sum rule abundantly used since the work of [START_REF] Hubbell | The unified neutral theory of biodiversity and biogeography[END_REF] should not be seen as a panacea for modelling community dynamics, despite its attractive simplicity.

Selected models had an overall satisfactory fit of the six summary statistics used in the inference (Fig. 1, Table S2). For the BCI dataset, two p-values were below the threshold of 0.025 (those associated with the average S and total S tot species richness). However, when applying the Holm-Bonferroni correction for multiple testing, they were no longer significant (since the corrected thresholds become equal to 0.025/6 = 0.0041 and to 0.025/5 = 0.005 for the two lowest p-values respectively). Still, this model checking procedure highlights that the fitted model tends to underestimate the average species richness S, while it tends to overestimate the total among censuses species richness S tot . This model-data mismatch is due to a larger mortality rate of very rare species in the data compared to the one of more abundant species (data not shown), a pattern that could not be reproduced with the models studied. Interestingly, the selected model for BCI was also able to reproduce observed values of two additional summary statistics not included in the inferential procedure: the average Shannon's index H across censuses, and the average Bray-Curtis similarity BC between each census and the initial one (Fig. S3). It was also able to marginally reproduce the number of times N trans that a population increase between two censuses is followed by a population decrease or the reverse (Fig. S3). These complementary results strengthen the demonstration of the overall good fit quality of the selected model regarding communitylevel static and dynamical properties. Similar results were obtained for the Mudumalai dataset (Table S2). For the Pasoh dataset, even the full model had a small but significant model-data mismatch for the average number of species S (Table S2). It is worth noting that for both BCI and Mudumalai, several models had good quality fits according to this model checking procedure (Table S2). The Bayesian model selection procedure still succeeded in clearly ranking these competing models (Table 2).

Author-produced version of the article published in Oikos, 2016, 125 (12), 1733-1742. The original publication is available at http://www.oikosjournal.org/ doi : 10.1111/oik.03103

For the BCI dataset, estimations of the mortality rate d and of the temporal autocorrelation ρ had low uncertainty with narrow posterior distributions (Fig. 2a,c). The environmental variance A was also well estimated, although with a wider posterior distribution (Fig. 2c). The parameters of immigration m and of regional diversity θ had wider posterior distributions (Fig. 2b), which is in line with previous studies (Etienne et al. 2006, Jabot and[START_REF] Jabot | Inferring the parameters of the neutral theory of biodiversity using phylogenetic information and implications for tropical forests[END_REF].

Similar results were obtained with the two other datasets, with the additional result that the community carrying capacity parameter K also had a relatively narrow posterior distribution for these two datasets (Fig. S4-5).

We found larger environmental variances (A = 2.2, Fig. 1f) for BCI than in the study of [START_REF] Kalyuzhny | A neutral theory with environmental stochasticity explains static and dynamic properties of ecological communities[END_REF] in which it was estimated at 0.8. This large environmental variance was also recovered in the Mudumalai dataset (A = 3), but less so in Pasoh (A = 0.94, Fig. S4-5). We also found that parameters A and ρ were negatively correlated (Fig. 2c). This last result explains the discrepancy between our estimate of A and the one of [START_REF] Kalyuzhny | A neutral theory with environmental stochasticity explains static and dynamic properties of ecological communities[END_REF] who considered a strong temporal autocorrelation in their model in which environmental variables were assumed constant during 10 years. Our analysis thus demonstrates that a larger environmental variance associated with a less strong temporal autocorrelation is more likely (Fig. 2c).

Assessing the contribution of environmental stochasticity to population variance

Values of temporal autocorrelation were variable among the three datasets, varying from 0.60 at Pasoh to 0.68 at BCI and 0.92 at Mudumalai. The resulting contributions of environmental and demographic stochasticity to population variance were found to be consistent across the studied sites. In BCI, where no trend in community size was observed, demographic stochasticity was the main driver of population variance for rare species, and it was progressively surpassed by the influence of environmental stochasticity for more abundant species (Fig. 3a). In Mudumalai where the community size decreased by more than 16% in 12 years, the trend in community size explained the larger part of population variance (between 49 and 77% depending on the abundance class, Fig. 3b). The remaining

Author-produced version of the article published in Oikos, 2016, 125 (12), 1733-1742. The original publication is available at http://www.oikosjournal.org/ doi : 10.1111/oik.03103

part of the population variance was again largely explained by demographic stochasticity for rare species, while environmental stochasticity was the main driver for more abundant species (Fig. 3b). In Pasoh where the community size increased by 6% in 13 years, the trend in community size also explained a large part of population variance (Fig. 3c). The remaining part was again mainly due to demographic stochasticity for rare species, and to environmental stochasticity for abundant ones (Fig. 3c).

Non-random correlations between species dynamics

Observed correlation coefficients between population trajectories tended to be larger in absolute values than the ones computed on simulated community dynamics, with a larger number of strongly positively or negatively correlated population trajectories (Fig. 4a). This trend was significant for the BCI dataset (p<0.005, Fig. 4b). This indicates that observed levels of correlations between species dynamics cannot be explained by the time-averaged neutral model in this dataset. Non-random correlations between species in their responses to environmental variability thus need to be invoked. In the two other datasets, the same test led to the same trends of stronger between-species correlations in empirical data than awaited under the time-averaged neutral model (Fig. 4c,d). These trends were however nonsignificant (p=0.09 in Mudumalai, p=0.24 in Pasoh, Fig. 4d). One possible explanation for this lack of significance for these two datasets may be that these communities were studied during a period of directional trend in community size which seems to be an important driver of the dynamics of these two communities (Fig. 3). Another possible explanation is the shorter length of the time series for these two datasets. Indeed, we found that deviations from model predictions were increasing with the number of censuses considered in the analyses for the BCI dataset (Fig. 5).

Discussion

Several recent studies have proposed to model plant community dynamics with both demographic and environmental stochasticity [START_REF] De Mazancourt | Predicting ecosystem stability from community composition and biodiversity[END_REF][START_REF] Kalyuzhny | A neutral theory with environmental stochasticity explains static and dynamic properties of ecological communities[END_REF], following the seminal presentation of [START_REF] Lande | Stochastic population dynamics in ecology and conservation[END_REF]. To further advance our

Author-produced version of the article published in Oikos, 2016, 125 (12), 1733-1742. The original publication is available at http://www.oikosjournal.org/ doi : 10.1111/oik.03103

understanding of the consequences of environmental stochasticity on community assembly and dynamics, rigorous inference methods need to be developed for these models, as was done previously with neutral and non-neutral dynamical community models (Beeravolu et al. 2009, Jabot and[START_REF] Jabot | Analyzing tropical forest tree species abundance distributions using a nonneutral model and through approximate Bayesian inference[END_REF]. The present contribution provides a quantitative Bayesian inferential framework for slightly modified versions of the time-averaged neutral model recently proposed by [START_REF] Kalyuzhny | A neutral theory with environmental stochasticity explains static and dynamic properties of ecological communities[END_REF].

We showed that it is possible to estimate both the variance A and the temporal autocorrelation ρ of environmental stochasticity from community time series of marked individuals (Fig. 2), and we provided a simple way to assess the importance of environmental stochasticity through a model selection technique. We recovered the results of [START_REF] Chisholm | Temporal variability of forest communities: empirical estimates of population change in 4000 tree species[END_REF] and of [START_REF] Kalyuzhny | A neutral theory with environmental stochasticity explains static and dynamic properties of ecological communities[END_REF] that environmental stochasticity was critical to understand the dynamics of tropical tree communities (Table 2). We further showed that environmental stochasticity was mostly important for abundant species, while rare species were mostly impacted by demographic stochasticity (Fig. 3). This may contribute to the fact that hyper-diverse tropical forests with many rare species often show neutral-like patterns [START_REF] Jabot | Analyzing tropical forest tree species abundance distributions using a nonneutral model and through approximate Bayesian inference[END_REF][START_REF] Kubota | Non-neutrality in forest communities: evolutionary and ecological determinants of tree species abundance distributions[END_REF][START_REF] Qiao | A latitudinal gradient in tree community assembly processes evidenced in Chinese forests[END_REF].

Our analysis further revealed that environmental variance was larger and temporal autocorrelation in environmental conditions lower (Fig. 2) than previously found by [START_REF] Kalyuzhny | A neutral theory with environmental stochasticity explains static and dynamic properties of ecological communities[END_REF]. Still, the presence of relatively large year-to-year autocorrelation in species demographic rates, already evidenced by [START_REF] Feeley | Directional changes in the species composition of a tropical forest[END_REF] in BCI, calls for a biological explanation. Indeed, climatic variability does not show large year-to-year autocorrelation [START_REF] Lande | Stochastic population dynamics in ecology and conservation[END_REF], so that other sources of autocorrelation needs to be invoked. First, gap dynamics is likely to produce long standing (decadal) conditions favoring some species over others in tree communities, depending on species light requirements and their sensitivities to crowding [START_REF] Shugart | A theory of forest dynamics: the ecological implications of forest succession models[END_REF][START_REF] Uriarte | A spatially explicit model of sapling growth in a tropical forest: does the identity of neighbours matter?[END_REF]). Second, the long generation time of trees and their capacity to store resources (Chapin et al. 1990, Poorter and[START_REF] Poorter | Carbohydrate storage and light requirements of tropical moist and dry forest tree species[END_REF] is likely to dampen their response to yearly environmental fluctuations and thus to increase the autocorrelation of the stochastic environmental variables. Third,

Author-produced version of the article published in Oikos, 2016, 125 (12), 1733-1742. The original publication is available at http://www.oikosjournal.org/ doi : 10.1111/oik.03103

synchronizing processes such as mast fruiting [START_REF] Norden | Mast fruiting is a frequent strategy in woody species of Eastern South America[END_REF] or extreme climatic events [START_REF] Brando | Abrupt increases in Amazonian tree mortality due to drought-fire interactions[END_REF]) may also tend to increase temporal autocorrelation. The inferred strong autocorrelation in stochastic environmental variables is thus biologically plausible.

We must however recognize that our inference of autocorrelation strength present a nonnegligible level of uncertainty (Fig. 2c, Fig. S4c, Fig. S5c) and that it is necessarily dependent on our (reasonable) model assumptions. Reported autocorrelation levels should thus be interpreted and used with caution.

Previous studies have investigated the signature and magnitude of environmental stochasticity in various types of ecological communities, using various mathematical approaches. For instance, model-based analyses of community time-series with likelihood [START_REF] Ives | Estimating community stability and ecological interactions from timeseries data[END_REF] or Bayesian techniques [START_REF] Mutshinda | What drives community dynamics?[END_REF][START_REF] Mutshinda | A multispecies perspective on ecological impacts of climatic forcing[END_REF] have been proposed.

These parameter-rich approaches require relatively long time-series so as to obtain reliable inferences and were consequently applied to a limited number of systems that have been intensively studied for a long time, with typically several dozens of censuses. Using these techniques, [START_REF] Mutshinda | What drives community dynamics?[END_REF][START_REF] Mutshinda | A multispecies perspective on ecological impacts of climatic forcing[END_REF] demonstrated that environmental stochasticity exerted a much stronger influence on community dynamics than demographic stochasticity and interspecific interactions for a diverse range of animal communities. Similar conclusions were obtained using model-free analyses of dynamical patterns for both plant and animal communities (e.g., [START_REF] Houlahan | Compensatory dynamics are rare in natural ecological communities[END_REF][START_REF] Chisholm | Temporal variability of forest communities: empirical estimates of population change in 4000 tree species[END_REF], Kalyuzhny et al. 2014a,b). The main advantages of these last approaches are that 1) they are not restricted to massive datasets, and 2) they are more flexible, in that they do not rely on necessarily debatable model assumptions. They constitute in this sense attractive ways to design pertinent model assumptions that can be later tested with model-based approaches. Inferences based on time-averaged neutral models are thus complementary to these existing approaches, in that they are a model-based approach that is "parameter-economical", and can thus be used on less massive datasets. They corroborated in this study the previous findings of the prevailing role of environmental stochasticity in shaping community dynamics, especially for abundant species (Fig. 3).

Author-produced version of the article published in Oikos, 2016, 125 (12), 1733-1742. The original publication is available at http://www.oikosjournal.org/ doi : 10.1111/oik.03103

A critical simplifying assumption of the time-averaged neutral model is that species respond independently of each other to environmental variability. The second main goal of this study was thus to assess the validity of this strong assumption of species independence. We performed this assessment using a model checking procedure comparing the correlation matrix of species trajectories in both observed and best-fit simulated communities. This assumption was rejected for the BCI dataset (p<0.005), but not for the two other ones (Fig. 4). These results are unlikely to be biased by the (non-significant) model-data mismatch surrounding the dynamics of very rare species, since only abundant species (>49 individuals)

were considered in the correlation matrix. These results show that non-random correlations between species dynamics need to be invoked to explain observed community dynamics, at least in some study sites. They nevertheless show that deviations from the predictions of the time-averaged neutral model are not systematic, neither very large. Indeed, at BCI, between-species correlation was found to be equal to 0.72 on average, a number only modestly larger than the predictions of the time-averaged neutral model that were located around 0.6 (Fig. 4b).

We then tried to understand what could explain this relatively modest failure of the model predictions on between-species correlations. We evidenced that the deviations from model predictions build over very long time, thereby making short time series insufficient to detect statistically significant deviations (Fig. 5). A more positive interpretation of these last results is that the simplifying assumption of species independence proposed by [START_REF] Kalyuzhny | A neutral theory with environmental stochasticity explains static and dynamic properties of ecological communities[END_REF] seems to be sufficient for some communities, and up to twenty year-long temporal horizons.

An important question is whether one can use time-averaged neutral models in communities in which the assumption of independence is violated, as in BCI. The present work and the previous one of [START_REF] Kalyuzhny | A neutral theory with environmental stochasticity explains static and dynamic properties of ecological communities[END_REF] shows that the model is still valid if one wants to make general predictions on community-wide dynamical patterns such as average population variability and its scaling with population abundance, or the respective contributions of demographic and environmental stochasticity to the observed community [START_REF] May | Moving beyond abundance distributions: neutral theory and spatial patterns in a tropical forest[END_REF] and anthropogenic impacts on biodiversity [START_REF] Solé | Habitat fragmentation and biodiversity collapse in neutral communities[END_REF][START_REF] Gilbert | Can Neutral Theory Predict the Responses of Amazonian Tree Communities to Forest Fragmentation?[END_REF][START_REF] Hubbell | How many tree species are there in the Amazon and how many of them will go extinct?[END_REF], as well as how diversity indices may depart from null predictions incorporating environmental stochasticity [START_REF] Münkemüller | From diversity indices to community assembly processes: a test with simulated data[END_REF]. A decade of research did not deplete the study of neutral models based on demographic stochasticity [START_REF] Rosindell | The unified neutral theory of biodiversity and biogeography at age ten[END_REF]. This myriad of works now needs to be generalized to embrace the combined effect of demographic and environmental stochasticity.

As a concluding remark, it may sound abusive to use the terminology of "time-average neutrality" to designate a model including both demographic and environmental stochasticity, since what drives environmental stochasticity is classically referred to as niche processes (as opposed to neutral processes) in the ecological literature. Actually, the distinction between "strong" and "weak" versions of the equivalence assumption of the neutral theory was stressed early on during the maturation of this theory [START_REF] Chave | Neutral theory and community ecology[END_REF],

where "strong" neutrality refers to the strict fitness equivalence among individuals, while "weak" neutrality refers only to average fitness equivalence among individuals, as in timeaveraged neutral models. A good reason for calling neutral these "weakly" neutral models is that individuals still behave neutrally from an evolutionary standpoint. Beyond this semantic discussion, the interesting originality of time-averaged neutral models is their way of merging what is commonly referred to as "niche" and "neutral" processes. Indeed, previously proposed syntheses of these two types of processes are commonly relying on the addition of neutral-based noise to niche-based processes that generate average fitness

Author-produced version of the article published in Oikos, 2016, 125 (12), 1733-1742. The original publication is available at http://www.oikosjournal.org/ doi : 10.1111/oik.03103 differences among species (e.g., [START_REF] Tilman | Niche tradeoffs, neutrality, and community structure: a stochastic theory of resource competition, invasion and community assembly[END_REF][START_REF] Gravel | Reconciling niche and neutrality: the continuum hypothesis[END_REF][START_REF] Jabot | Analyzing tropical forest tree species abundance distributions using a nonneutral model and through approximate Bayesian inference[END_REF]. In contrast, time-averaged neutral models integrate niche-based differences among species without breaking the weak version of the equivalence assumption among individuals [START_REF] Alonso | Response to Benedetti-Cecchi: neutrality and environmental fluctuations[END_REF][START_REF] Kalyuzhny | A neutral theory with environmental stochasticity explains static and dynamic properties of ecological communities[END_REF]. respectively, highlighting the correlations between these parameters. In these panels, points represent posterior samples obtained by ABC, the grey levels stand for posterior density, and the black lines are iso-density curves.

Author-produced version of the article published in Oikos, 2016, 125 (12), 1733-1742. The original publication is available at http://www.oikosjournal.org/ doi : 10.1111/oik.03103 The proportion of population variance due to demographic stochasticity is depicted in black, the one due to environmental stochasticity is in dark grey, and the variance due to the trend in community size is in light grey.

Author-produced version of the article published in Oikos, 2016, 125 (12), 1733-1742. The original publication is available at http://www.oikosjournal.org/ doi : 10.1111/oik.03103 Author-produced version of the article published in Oikos, 2016, 125 (12), 1733-1742. The original publication is available at http://www.oikosjournal.org/ doi : 10.1111/oik.03103 Author-produced version of the article published in Oikos, 2016, 125 (12), 1733-1742. The original publication is available at http://www.oikosjournal.org/ doi : 10.1111/oik.03103

  For each dataset, we performed four sets of 20, 000 ABC simulations with the four versions of the model: the full model with six parameters, and the three versions of the model without trend in community size: the models with environmental stochasticity with (ρ > 0) and without (ρ = 0) within-species autocorrelation, and the model without environmental stochasticity (A = 0 and ρ = 0). Uniform prior distributions were used: [0.01 ; 0.025] for d, [0.001 ; 5] for A, [0 ; 1] for ρ, [ln(0.001) ; ln(0.1)] for ln(m) and [ln(1) ; ln(10000)] for ln(θ).

Figure 2 .

 2 Figure 2. Posterior distributions for the BCI dataset. Panel A: marginal posterior distribution

Figure 3 .

 3 Figure 3. Contributions of demographic and environmental stochasticities to population

Figure 4 .

 4 Figure 4. Non-random correlations between species dynamics. Panel A: Distribution of

Figure 5 .

 5 Figure 5. Lower model failure with shorter time series. Panels A, B: Distributions of the

  Author-produced version of the article published inOikos, 2016, 125 (12), 1733-1742.The original publication is available at http://www.oikosjournal.org/ doi : 10.1111/oik.03103 dynamics. In contrast, this model is unlikely to make good predictions in such cases regarding between-species correlations in their demography.

	More generally, the present contribution employs a quantitative statistical framework to
	analyze community time series with time-averaged neutral models. Being able to quantify
	environmental stochasticity in ecological communities from community time series opens
	exciting new questions for ecological research. Theoretical predictions could be developed
	on how environmental stochasticity might shape spatial and temporal patterns in
	metacommunities (Hubbell 2001, Chave et al. 2002, Ulrich 2004, Economo and Keitt 2010,

Table 1 .

 1 Characteristics of the tropical tree datasets.

			Number of	Census	Species
	Site	Site area		Community size
			censuses used	date	richness
	Barro
	Colorado	50 ha	5
	Island

Table 2 .

 2 Posterior model probabilities for the three datasets. For each dataset, the probability of the selected model is in bold.Author-produced version of the article published inOikos, 2016, 125 (12), 1733-1742.The original publication is available at http://www.oikosjournal.org/ doi : 10.1111/oik.03103 Figure 1. Model checking plots for the BCI dataset and the selected model (see Table 2). Histograms represent the simulated values of the various summary statistics for the 200 best-fit simulations. Observed values are depicted by the vertical grey lines, and the corresponding p-values are indicated.Author-produced version of the article published inOikos, 2016, 125 (12), 1733-1742. The original publication is available at http://www.oikosjournal.org/ doi : 10.1111/oik.03103

	Site	BCI	Mudumalai	Pasoh
	Dem. stochasticity (A = 0 and ρ = 0)	0	0	0
	Dem. and env. stochasticity (ρ = 0)	0.03	0	0
	Dem. and env. stoch. with autocorrelation	0.63	0.03	0.0002
	Full model with trend in community size (K)	0.35	0.97	0.9998

Acknowledgements

F.J. was funded by the Irstea Dynindic project. T.L. was funded by Irstea and the Regional Council of Auvergne. We thank the Subject Editor and three anonymous reviewers for their numerous suggestions that greatly improved this manuscript. We thank the volunteers, workers and data providers of the Center for Tropical Forest Science for having made their datasets publicly available.

Author-produced version of the article published in Oikos, 2016, 125 (12), 1733-1742.

The original publication is available at http://www.oikosjournal.org/ doi : 10.1111/oik.03103