domain: math info

text:
Hubert Baty
email: hubert.baty@unistra.fr

Léo Baty

Solving differential equations using physics informed deep learning: a hand-on tutorial with benchmark tests

We revisit the original approach of using deep learning and neural networks to solve differential equations by incorporating the knowledge of the equation. This is done by adding a dedicated term to the loss function during the optimization procedure in the training process. The so-called physics-informed neural networks (PINNs) are tested on a variety of academic ordinary differential equations in order to highlight the benefits and drawbacks of this approach with respect to standard integration methods. We focus on the possibility to use the least possible amount of data into the training process. The principles of PINNs for solving differential equations by enforcing physical laws via penalizing terms are reviewed. A tutorial on a simple equation model illustrates how to put into practice the method for ordinary differential equations. Benchmark tests show that a very small amount of training data is sufficient to predict the solution when the non linearity of the problem is weak. However, this is not the case in strongly non linear problems where a priori knowledge of training data over some partial or the whole time integration interval is necessary.

I. INTRODUCTION

Neural networks (NN) are widely used to solve problems in a variety of domains including computer vision, language processing, game theory, etc., as one can see in Le Cun et al. (2015) and references therein. The use of machine learning approaches in the field of scientific computing including differential equations is relatively recent. Indeed, the idea of leveraging prior knowledge of the physics in the learning process of a NN network was introduced by Raissi et al. (2017,2019).

Among other things, NN are a tool that can be used for supervised learning, one of the main machine learning settings. Supervised learning consists in finding a mapping function between given input objects and their associated output values. This is done by using knowledge about a dataset containing several input/output pairs. This dataset is used to parameterize the NN such that it minimizes the error between solutions predicted by the NN and true known solutions in the dataset. The convergence is achieved by minimizing a loss function which expression is based on the mean squared error. Finding "good" parameters is achieved by solving an optimization problem using a gradient algorithm that relies on automatic differentiation to back-propagate gradients through the network (Baydin et al. 2018).

In the case of differential equations, we can apply the supervised learning setting. Indeed, solving a given differential equation comes down to finding a mapping function between some physical input variable values (position, time, . . . ) and a corresponding unknown physical quantity which is the solution of the equation. Hence, by training a NN we can obtain a non-linear approximation of the solution, which can be used to instantaneously predict the equation's solution at any given input point. However, a first strong limitation comes from the impossibility to extrapolate the desired solution for input variable values situated outside the range of the training data. In other words, the NN is a bad extrapolation function.

Second, a minimum amount of training data is required, otherwise, wrong solutions or even absence of convergence can occur during the training process.

In order to tackle these limitations, classical NN can be enhanced by giving it additional information corresponding to the physics. These approaches are generally called physicsinformed neural networks (PINNs), in the context of simulating physical and engineering systems modeled by differential equations. The method consists in evaluating the solution at some other set of data points (called collocation points) at which the estimated solution must ensure the equations. A new loss function corresponding to the physics is thus defined and added to the previous one in the learning process. In other words, the training is penalized by this additional constraint. The space of available solutions is thus restricted, being partly driven by the original data and also partly driven by the physics. When solving partial differential equations (PDEs), one is particularly interested in using only a very minimal known data set, as for example the solution on the boundary and at the initial time. Then, in this sense, the PINNs approach can be said to be mostly physically-driven, as opposed to data-driven. In this work, we focus on such motivation, even if PINNs can also be used for many other aims like inverse or physics discovery problems (see the discussion and conclusion in the last section). One can also refer to Cuomo et al. (2022) and Karniadakis et al. (2021) for reviews.

The paper is organized as follows. We first review the basics of PINNs for PDEs in Section 2. Section 3 is devoted to a tutorial of the method to solve a simple first order ordinary equation. The results of benchmark tests performed on a series of different academic ordinary differential equations (with increasingly non-linearity) are presented in Section 4. Finally, a discussion and conclusions are drawn in Section 5, with a particular attention on highlighting the advantages and drawbacks of the PINN approach versus standard integration schemes.

II. PHYSICS-INFORMED NEURAL NETWORKS

A. The basics of PINNs for PDE

We consider a partial differential equation (PDE) written in the following residual form

F(x, t, u x , u t , ...) = 0, x ∈ Ω, t ∈ [t 0 , T ] , (1)
with the imposed initial condition u(x, t 0 ) = u 0 (x). A boundary condition must be also specified as, u(x, t) = u Γ (t) for x ∈ ∂Ω (a Dirichlet-like condition is chosen for simplicity). Ω and ∂Ω represent the spatial domain and associated contour, respectively. The variables x ∈ R d (d being a spatial dimension) and t ∈ R denote the spatial coordinate and time, respectively. Note that higher order differential operators can be also included in correspondence with multidimensional system of equations instead of Equation 1.

We introduce a neural network approximating the desired solution u(x, t) of the PDE with ûθ u, where θ is a set of model parameters. An artificial feed forward neural network is taken (see Figure 1), with neurons organized in different layers in order to perform calculations in a sequential way. A single input layer containing the input variables (x, t) is connected to a few hidden layers (two layers with four neurons in the schematic example of Figure 1), and finally to an output layer for the solution ûθ . Neurons are only connected in adjacent layers, and are not linked inside each layer. The neural network of L + 1 layers is a non linear approximation function, N L (z), that can be expressed recursively as follows.

The input vector is denoted by z ∈ R d i , with d i = d + 1, as it includes the spatio-temporal coordinates (x, t). Thus,

N 0 (z) = z. (2)
For the hidden layers (1 ≤ l ≤ L -1), we have

N l (z) = σ(W l N l-1 (z) + b l ), (3)
where we denote the weight matrix and bias vector in the l-th layer by W l ∈ R d l-1 ×d l and b l ∈ R d l (d l being the dimension of the input vector for the l-th layer). σ(.) is a non linear activation function, which is applied element-wisely. In this work, me choose the most commonly used hyperbolic tangent tanh function. For the (final) output layer, we get

N L (z) = W L N L-1 (z) + b L , (4)
and finally ûθ (x, t) = N L (z). The network can be also written as a sequence of non linear

functions ûθ (x, t) = (N L • N L-1 ... N 0 )(x, t), (5)
where the operator • denotes the composition and θ = {W l , b l } l=1,L represents the trainable parameters (weight matrices and bias vectors) of the network. The goal is to calibrate its parameters θ such that ûθ approximates the target solution u(x, t).

B. The training of PINNs for PDE

The resolution of the PDE is reduced to an optimization problem as follows, and as schematized in Figure 1. We first assume that a set of N data data is available for the known solution at different times, i.e. {t i data ,

x i data , u i data } N data i=1
that are the training data, which Two hidden layers with four neurons per layer are connected with the input and the output layer,

where the latter has a single variable (one neuron) representing the desired solution ûθ .

include the initial condition. A corresponding loss function L data (using the mean square error formulation) can be deduced from the residual as

L data (θ) = 1 N data N data i=1 ûθ (z i ) -u i data 2 . ( 6
)
In a similar way, defining a loss function L b corresponding to the knowledge of the boundary condition, we have

L b (θ) = 1 N b N b i=1 ûθ (z i ) -u i b 2 , (7)
where a set of N b known data is imposed via

{t i b , x i b , u i b } N b i=1 .
Finally, another loss function for the equation itself can be also obtained as,

L F (θ) = 1 N c Nc i=1 F( ûθ (z i )) 2 , (8)
that must be evaluated on a set of N c data points (generally called collocation points) as explained below. Indeed, one advantage of the neural network approach is given by the possibility to evaluate exactly the differential operators at the collocation points in L F and F by using automatic differentiation. The automated differentiation is also used to compute derivatives with respect to the network weights (i.e. θ), that is necessary to implement the optimization procedure (see below). Note that in this way, the derivatives can be obtained at machine precision, contrary to the use of some standard numerical scheme. Moreover, the latter operations are greatly facilitated by Python open source software libraries like Tensorflow or Pytorch.

A composite loss function is generally formed as

L(θ) = ω data L data (θ) + ω b L b (θ) + ω F L F (θ), (9)
where an optimal choice of values for hyper-parameters (ω data , ω b , ω F ) allow to ameliorate the eventual unbalance between the partial losses during the training process. These weights can be user-specified or automatically tuned. In the present work, for simplicity we fix the ω data value to be constant and equal to unity, and the other weight parameters are determined with values varying from case to case.

A gradient descent algorithm is used until convergence towards the minimum is obtained for a predefined accuracy (or a given maximum iteration number) as

θ i+1 = θ i -η∇ θ L(θ i ), (10)

C. The ODE case

In this study, we focus on ordinary differential equations (ODE). Thus, the spatial dependence is ignored in Equation 1, and we are left with a desired solution u(t) and its approximation ûθ (t). The input first layer of the neural network is supplied with N data values at different times t i ∈ [t 0 , T ] corresponding to u data (t i ) = u i data (with i = 1, N data ). The initial condition corresponds to the first point, u data (t 0 ) = u 1 data . In this work, for simplicity, we assume a uniform distribution of the N data points within a subinterval of the full time interval domain. As a consequence of the ODE particular case, we are not concerned with the boundary condition, thus ω b = 0. However, the central part of the PINNs concept concerns the loss function L F (θ), that is evaluated at N c collocation points which are not necessarily coinciding in time values with N data . The distribution of N c is also taken to be uniform in the full time domain, or on a subinterval of it. Some of the Pytorch Python-based codes and data-sets accompanying this manuscript are available on the GitHub repository at https://github.com/hubertbaty/PINNS-EDO. These have been inspired by the codes provided on https://benmoseley.blog and available on GitHub repository at https://github.com/benmoseley/harmonic-oscillator-pinn. We have chosen to use very simple deep feed-forward networks architectures with hyperbolic tangent activation functions. In this work, the optimal choice of detailed architecture of the network (number of hidden layers, number of neurons per layer) and of hyperparameters (learning rate, loss weights) is done manually. Although more systematic/automatic procedures could be used, this is a more complicated task not considered in this work.

Note that we use the notation y for the desired ODE solution in the following instead of u introduced in the previous section for a PDE case.

III. ILLUSTRATION OF THE METHOD ON A SIMPLE TUTORIAL EXAMPLE A. The differential equation: a tutorial example

Let us consider the following equation example, called tutorial equation below, to be solved for y(t), dy dt + 0.1t -sin(πt/2) = 0, (11) fo t ∈ [0, 30], using the initial condition y 0 = y(0) = 1. As can be seen below, the corresponding solution contains two time scales, a first one due to the sinusoidal forcing term, and a second one due to the linear term in t 10 that gives an exponentially decreasing envelope amplitude.

B. Solving with a normal neural network

We first consider a situation without any constraint coming from the differential equation, i.e. ω F = 0. As we are left with data coming only from the exact solution, we call it the normal neural network. The training procedure is illustrated on Figure 2 for two training steps (i.e. n t ), at n t = 4000 and n t = 24000. The solution is not fully converged for n t = 4000, contrary to n t = 24000. This is in agreement with the history of the loss function (see Figure 3), where the convergence is already roughly obtained when n t 15000. A learning rate of η = 3 × 10 -3 , with ω data = 1 is chosen. A choice of 3 hidden layers with 32 neurons per layer is also done. The exact solution that is drawn for comparison, is obtained using a classical Runge-Kutta method of order two. Note that the later is also useful to extract the training data corresponding to the N data values, y i data . Figure 2 clearly illustrates the ability of the normal NN to approximate the solution for a relatively high number of N data values, as N data = 101 is employed.

However, when the number of training points is not enough, the convergence towards the solution is bad, or it can completely fail. This is illustrated in two cases at the end of the convergence process stopped at n t = 48000. Indeed, the results for the first case obtained for N data = 26 values uniformly chosen within the full time interval, show bad convergence properties as illustrated in Figure 4 (left panel). Moreover, when the points are distributed only within a subinterval (with N data = 61), Figure 4 (right panel) shows that the method completely fails to obtain an acceptable solution in the other subinterval that is free of training data. In other words, the NN is not able to correctly extrapolate the solution in these two cases.

C. Solving with PINNs

In this sub-section, we consider now the possibility to add the constraint on the loss function with a non zero contribution coming from the differential equation at some chosen collocation points, i.e. ω F = 0. In other words, we minimize by adding a weighted partial loss function term ω F L F (θ), where F = dy dt + 0.1t -sin(πt/2). Typically, we define a uniform data set of N c = 50 points within the full time interval. The results of the training process is illustrated for the two cases previously studied using the normal network, for which the convergence is not satisfaying. The results plotted in left/right panels of Figures 5,display a spectacular amelioration after n t = 40000 iteration steps in both cases. Additionally, Figure 6 shows that taking collocation points only within a right subinterval (for the second above case) can also be sufficient. Note that, a minimum value for N c is required with an exact value that depends on the parameters of the network (i.e. number of layers, neurons, etc.).

The distribution of the collocation points can also influence the results, but the philosophy behind the PINNs technique remains. The learning rate together with the loss weight values can also influence the convergence of the gradient descent algorithm. Indeed, a too high value of η leads to strong oscillations in the loss function, whilst a too small value can induce a very small convergence speed. The combination of relative weights (ω data and ω F ) is also important, in order that the two partial losses converge at a similar rate. For the example studied above, we have taken optimal values ω data = 1, and ω F = 6 × 10 -2 for the weights of the two partial losses.

If we reduce the training data amount to the minimum possible, i.e. only one point corresponding to the initial value is taken, the solution obtained is also excellent as one can see in left panel of Figure 7. The corresponding loss function is also plotted in right panel of Figure 7. Note that, in this case, an optimal choice of 4 hidden layers with 32 neurons per layer is done. Now the question is, does the method work so nicely for any (ODE) differential equation for which a very reduced amount of data is known. In order to answer this question, we investigate a rather large number of academic cases in the following section, including the important class of second order differential equations.

IV. DIFFERENT EXAMPLES -BENCHMARK TESTS A. Harmonic oscillator

In this sub-section, we first consider the following harmonic oscillator equation,

d 2 y dt 2 + ω 2 0 y = 0, ( 12
)
where ω 0 is the normalized angular frequency, and where the time domain considered is t ∈ [0, 1]. We also take the initial conditions y(t = 0) = y 0 = 1 and dy dt (t = 0) = 0. Thus, the exact solution is a simple cosine function, i.e. y(t) = cos(ω 0 t). Our PINNs algorithm is first used to integrate the oscillator equation for ω 0 = 20, i.e. for a time interval slightly larger than 3 periods. Note that, in order to evaluate the corresponding ODE loss function L F , a second order automatic differentiation must be used for this example. The choice of the activation function is important to this respect (i.e. the hyperbolic tangent in this study).

The solution predicted when only one point corresponding to the initial value (N data = 1) is imposed, is plotted in left panel of Figure 8. It is clearly bad for the parameters taken in this case. The quality of the result slightly varies with this choice, but taking only one training data value is in general not sufficient to lead to an acceptable solution. This is not surprising, as a classical integration method (either analytical or numerical) requires two initial conditions for such second order differential equation.

The quality of the solution is greatly ameliorated by adding a second data value (i.e.

N data = 2). Another second improvement comes from adding another constraint, that is the conservation of the total energy E = 1 2 ( dy dt ) 2 + 1 2 ω 2 0 y 2 (up to addition by a constant which value is determined by the initial conditions). Indeed, we can add to the total loss another partial loss function defined as

L E (θ) = 1 N c Nc i=1 E( ûθ (z i )) -E 0 2 , (13)
with a corresponding weight ω E to be determined and E 0 being the initial constant total energy. More explicitly, we have L(θ) = ω data L data (θ) + ω E L E (θ) + ω F L F (θ). These new results are visible in right panel of Figure 8 for the predicted solution obtained at n t = 54000.

Moreover in Figure 9, one can compare the corresponding mean square error (M SE) and loss histories evaluated for 1000 points taken within the full time domain. One can clearly see the improvement with a lower minimum M SE value by more than two orders of magnitude obtained in the improved procedure. This is remarkable, as this is despite the fact that the loss function history (also visible on the same figure) converges in a similar way in the two cases.

Note also that, with these two improvements we have obtained that a relatively low minimum value for the number of collocation points with N c 24 is sufficient. Such very small minimum number of collocation points is a great advantage compared to a standard integration method for which the small time step restriction (because of stability and/or precision) requires a much larger number of points within the time interval. The choice of the other hyperparameters are η = 3 × 10 -4 , ω data = 1, and ω F = ω E = 3 × 10 -4 . Three , where the evaluation ŷθ (t i ) is done on N eval = 1000 points uniformly distributed within the whole time interval, and y i eval is the expected exact solution at t = t i .

hidden layers with 32 neurons per layer are taken for the neural network architecture. We have also investigated longer time cases with higher ω 0 values, typically up to ω 0 130 (not shown). The conservation of energy is fundamental to this respect, otherwise the algorithm doesn't converge towards the expected solution for such long time integration.

B. Non linear pendulum

Second, as a natural extension of the harmonic we consider the non linear pendulum example below,

d 2 y dt 2 + ω 2 0 sin(y) = 0, ( 14
)
where ω 0 is the normalized angular frequency (a value ω 0 = 25 is taken below), and where the time domain considered is t ∈ [0, 1]. We also choose the initial conditions y 0 = 0.1 and dy dt (t = 0) = 40. As for the harmonic oscillator, we investigate the possibility to use one single point for N data corresponding to the initial condition y 0 , the conservation of the total energy E being imposed via a non zero ω E L E (θ) term with E = 1 2 ( dy dt ) 2 -ω 2 0 cos(y) (up to addition by a constant which value is determined by the initial conditions). The hyperparameters used are η = 1 × 10 -3 , ω data = 1, ω F = 3 × 10 -6 , and ω E = 3 × 10 -7 .

When the NN architecture is composed of 4 hidden layers with 32 neurons per layer, the converged solution obtained for 72000 epochs is correct. However, taking 3 layers (instead of 4) leads to a wrong solution even if the loss function displays convergence. Indeed, in the latter case, the calculated solution is shifted with respect to the exact solution, and the M SE is dramatically high. This is illustrated in Figures 10 and11.

Of course, taking two data points (i.e. N data = 2) allows a nice convergence towards the expected exact solution (not shown) for 3 and 4 hidden layers. And, as for the harmonic oscillator equation, the minimum required value for the number of colocation points N c remains rather low, as it is of order 35 now.

We have also investigated the possibility to use an equivalent form of a system of two first order differential equations, as done in analytical or classical numerical integrations.

Indeed, we can consider the following equivalent system:

     dy 1 dt -ω 0 y 2 = 0, dy 2 dt + ω 0 sin(y 1 ) = 0, (15)
where y 1 represents the desired solution (i.e. the previous y parameter) and y 2 is its associated time derivative dy dt divided by ω 0 . The advantage of the latter normalisation is important, as it facilitates the use of the NN network because in this way y 1 and y 2 have values of the same order of magnitude. Otherwise, two networks (one per variable) must probably be employed. We have thus used our PINNs algorithm with an input layer containing one neuron for t, 4 hidden layers with 32 neurons per layer, and a final output layer containing two neurons for y 1 and y 2 . The partial loss function for the data, L data , is calculated using two points at t = 0 now, one value for y 1 (t = 0) = 0.1 and one for y 2 (t = 0) = 1.6.

The procedure is equivalent to a standard integration scheme using one initial condition for the solution and one for its derivative. The partial loss for the ODE,

L F = L F 1 + L F 2 ,
is now the sum of two terms corresponding to the two equations respectively. The results of the previous example with N c = 40 are successfully obtained using the following hyper parameters, η = 3 × 10 -3 , ω data = 1, ω F = 1 × 10 -1 , and ω E = 2 × 10 -6 , as plotted in

C. Anharmonic oscillators

If we consider an anharmonic potentiel of the form y 4 4 (instead of y 2 2 for the harmonic oscillator), we get the corresponding ODE with a non linear restoring force ∝ y 3 ,

d 2 y dt 2 + ω 2 0 y 3 = 0. ( 16
)
The PINNs integration for this problem leads to results and conclusions very similar to the non linear pendulum case. Indeed, the use of two training data points with energy conservation constraint considerably ameliorate the convergence. Thus we have E = 1 2 ( dy dt ) 2 + 1 4 ω 2 0 y 4 (up to addition by a constant which value is determined by the initial conditions). This is illustrated in Figure 13 for a case with y 0 = 1.5 and ω 0 = 15.5. We have chosen the following hyper parameters, η = 1.5 × 10 -3 , ω data = 1, ω F = 1 × 10 -5 , and ω E = 1 × 10 -6 .

The choice of five hidden layers for the neural network seems to be a better optimal choice for this equation. Note also that for the latter example, a minimum number of collocation points of N c = 37 seems to be necessary, that is slightly higher than for the non linear pendulum. However, as we have not investigated a large range of different initial parameters, conclusion is not firm and only gives a rough tendency.

It is also instructive to consider a more complex dynamical system corresponding to a double well potential of the form y 4 4 -y 2 2 , and according to the differential equation,

d 2 y dt 2 + ω 2 0 (y 3 -y) = 0. (17)
For this example, two families of solution exist (see below). In order to easier impose the initial conditions, we also take the equivalent system of two equations (as done for the non linear pendulum),

     dy 1 dt -ω 0 y 2 = 0, dy 2 dt + ω 0 (y 3 1 -y 1 ) = 0. (18)
Indeed, taking the initial condition y(0) = y 0 = 1.8 , or equivalently y 1 (0) = 1.8, together with zero derivative y 2 (0) = 0, a first solution is obtained as plotted in Figure 14. We use ω 0 = 12. The predicted solution is nicely reproduced when compared to the exact expected solution (obtained using a Runge-Kutta integration). For this example, two training data points seems to be necessary (i.e. N data = 4), with the use of the energy conservation as

E = ω 2 0 ( y 2 2 2 + y 4 1 4 - y 2 1
2 ). We have also chosen the following hyper parameters, η = 1 × 10 -3 , ω data = 1, ω F = 6 × 10 -3 , and ω E = 6 × 10 -5 , but they must be adjusted from case to case.

A second solution corresponding to y(0) = y 0 = 1.38 (with zero initial derivative condition) is also trained and nicely obtained, as one can see in Figure 15.

However, taking an initial condition closer to the critical value y 0 = √ 2, our PINN algorithm fails to converge to the exact solution. This is not completely surprising, as this corresponds to a threshold separating solutions having orbits in the phase space trapped in the well centered on (y 1 = 0, y 2 = 0) with solutions having orbits centered on (y 1 = 1, y 2 = 0).

In other words, two types of solutions coexist for y 0 = √ 2.

D. Van Der Pol oscillator

We now consider the Van Der Pol (VDP) oscillator equation given by, where ω 0 is a normalized angular velocity, and is a parameter having a value which determines the amplitude of a limit cycle in the phase space (see below). Note that the harmonic oscillator is recovered for = 0. The particularity of the system is the existence of a limit cycle as illustrated in Figure 16, obtained using a classical Runge-Kutta integration (of order 4) for ω 0 = 15, and = 5, for t ∈ [0, 3]. When the parameter has a lower value, the system exhibits a limit cycle with less distortion in the phase space, as illustrated in Figure 17.

d 2 y dt 2 + ω 2 0 y -ω 0 (1 -y 2 ) dy dt = 0, (19)
We have thus used our PINNs algorithm to explore its behavior for three cases corresponding to three values of , i.e. = 1 3 , 1, and 5. We have chosen ω 0 = 15 for t ∈ [0, 1.5]. We also take the initial condition y 0 = 1. The results are plotted in Figure 18, for a neural network having three hidden layers with 32 neurons per layer. The number of training data and collocation points, as well as the weight ω F associated to the equation (ω data = 1 being fixed) are varying from case to case (see the legend). The results clearly show that increasing the non linearity (via the parameter) require a higher number of points. This is the case of the number of collocation points, but also of the number of training data points.

for the = 5 case we need a collection of training data points distributed within the whole time domain, while for the smallest case, three points at early times are sufficient to obtain a convergence towards the exact solution. This is clearly a strong limitation of the PINN algorithm when one want to solve highly non linear problems.

V. DISCUSSION AND CONCLUSION

In this work, we have reviewed the basic concepts of using neural networks in order to integrate differential equations. More specifically, we have focused on the use of the equations knowledge to penalize the convergence of the training process, and generally referred as physics-informed neural networks in the literature. A tutorial example on a simple ODE is presented, with the aim to illustrate how adding a partial loss function associated to the physics information (i.e. via the differential terms ensuring the equation) can considerably ameliorate the results of a normal neural network.

Benchmark tests on different second order ODEs are used in order to highlight the benefits and drawbacks of this approach when compared to a traditional numerical integration method. When the problem equation displays weak non linearity, the training procedure is successful using known data representing only the initial conditions (as for a classical integration method). The first advantage of PINNs in this case, is the need to use a very low number of collocation data points. Indeed, for the problems illustrated in this work, between 20 and 40 points are sufficient. An integration using a Runge-Kutta method (of order two) for the same equations would require a number of points higher by at least one order of magnitude. The second advantage is that, once trained the solution for a given time (case of an ODE) is instantaneously predicted. This is not the case for a classical integration for which a new integration procedure must be realized. However, when the non linearity is increased (ses VDP oscillator tests with increasing parameter), the knowledge of a higher amount of training data is required with also a higher number of collocation data points.

The previously cited benefits of PINNs are consequently reduced. Nevertheless, the second benefit remains. The other drawbacks of the method concern the lack of general automatic differentiation to evaluate the equation terms at collocation data points. There is already a wealth of variants to do it differently in order to improve the efficiency of the optimization procedure. As a very incomplete list of NN-based ideas to solve PDEs, we have methods based on, learning the solution map, and variational formulation (Ritz and Galerkin). In the context of solving differential equations, PINNs approach can be also used for other different problems. For example, it is useful for inverse problems when data are known at some time different from the initial conditions that are desired (Raissi et al. 2019). Indeed, a neural network does not make difference between the two boundaries of the time interval. Second, PINNs technique can be also used when some terms of the differential equations are not completely known, and which can be added as additional inputs in the input layer. The aim of the latter problems concerns the data-driven discovery of the governing equations (Rudy et al. 2017).

FIG. 1 .

 1 FIG. 1. Schematic representation of the structure for a Physics-Informed Neural Network applied to the resolution of a differential equation. The input layer has two input variables (i.e. two neurons) noted x 1 and x 2 representing for example a space and a time coordinate respectively.

  for the i-th iteration also called epoch in the literature, leading to θ * = argmin θ L(θ), where η is known as the learning rate parameter. In this work, we choose the well known Adam optimizer. A standard automatic differentiation technique is necessary to compute derivatives (i.e. ∇ θ ) with respect to the NN parameters (e.g. weights and biases) of the model(Raissi et el. 2019).

FIG. 2 .

 2 FIG. 2. Tutorial solution ŷθ (t) (red solid line) predicted by the normal NN for n t = 4000 iterations (in left panel) and n t = 24000 iterations (in right panel) respectively, and compared to the exact solution (blue hatched line). The chosen training data set values are indicated using circles (with N data = 101). The physical information is not used, i.e. ω F = 0.

FIG. 4 .

 4 FIG. 4. Tutorial solution obtained in two cases using the normal NN, for a low number of data (N data = 26) uniformly taken within the whole time domain (left panel), and for N data = 61 taken within a left subinterval (right panel). The convergence is stopped at n t = 48000.

FIG. 5 .

 5 FIG. 5. Tutorial solution obtained in the two cases of the previous figure with PINNs, where a set of N c = 50 collocation points is used to calculate an ODE loss function L F . The cases with N data = 26 data points taken within the whole interval, and with N data = 61 data values within a left subinterval, are plotted in left/right panel respectively. The time values (for collocation points) at which the physical loss function is evaluated are indicated with the small green circle on t axis.

(

  FIG. 9. Histories of the total loss function L(θ) and M SE corresponding to the two panels of the previous figure respectively. The M SE is evaluated using the standard expression, M SE =

FIG. 10 .

 10 FIG. 10. PINNs solution obtained for the non linear pendulum with N data = 1 and N c = 40 collocation points, using a NN architecture of 3 hidden layers (left panel) and 4 hidden layers (right panel).

FIG. 11 .

 11 FIG. 11. History of the total loss function L(θ) and M SE corresponding to the two cases of the previous figure (left and right panel for N data = 1 and N data = 2 respectively). The M SE is evaluated using the standard expression, M SE = 1 N eval

FIG. 12 .

 12 FIG. 12. PINNs solution obtained for the non linear pendulum using an equivalent formulation with two equations for the two variables y 1 and y 2 (see text) with N data = 2 (one point for y 1 and one for y 2 , i.e. for the two initial conditions t = 0) and N c = 40 collocation points. The energy conservation is also applied. Solutions as functions of time and in phase space are plotted in left and right panels, respectively.

Figure 12 . 2 2 2 -

 1222 Figure 12. The corresponding M SE is similar to the M SE previously obtained for the PINN solution of the single second order equation using N data = 2, with the constraint on the energy conservation via a L E partial loss function. Note that, the total energy with the system formulation is, E = ω 2 0

FIG. 13 .

 13 FIG. 13. PINN solution obtained for the anharmonic oscillator with N data = 2 and N c = 40.

FIG. 14 .

 14 FIG. 14. PINNs solution obtained for the double well oscillator system for the initial condition y 0 = 1.8. We use N data = 4 with 2 points for each variable, and N c = 40 collocation points. A NN architecture with five hidden layers is taken. The energy conservation also applied. Solutions are plotted as functions of time (left panel), and in phase space (right panel).

FIG. 16 .

 16 FIG. 16. Exact solution of the VDP oscillator equation obtained for = 5 and the initial conditions y(t = 0) = 1, dy dt (t = 0) = 0. The solution y(t) as a function of time is plotted in the left panel, and in phase space in the right panel.

ACKNOWLEDGMENTS

Hubert Baty thanks Emmanuel Franck, Victor Michel-Dansac, and Vincent Vigon (IRMA, Strasbourg), for associating him to the supervision of the Master2 internship of Vincent Italiano in February-July 2022, which also gave him want to learn the PINNs technique.

procedure for a fine tuning of the hyperparameters in order to have an optimal convergence during the training. Nevertheless, the most interesting point of the PINNs is its meshfree property, contrary to traditional integration methods. This approach remains relatively recent, and many ameliorations are expected in the next years. In this work, we have focused on the original PINNs method based on the automatic


Anaïs Carlin

Sébastien Verel

Philippe Collard

Modeling Luxury Consumption: An Inter-Income Classes Study of Demand Dynamics and Social Behaviors

Keywords: Preference formation, Consumption, Luxury good, Agent-based model, Income classes. JEL classification: B52, D01, D11

HAL is

Motivations and related literature

The study of luxury consumption is relatively rare in economic analysis. Actually, the standard economic theory of demand, based on the simplified concept of homothetic preferences, does not allow for such a study. The purpose of this article is to fill this gap by building a structure of preferences based on sociological and psychological factors, such as the need for social distinction.

Especially, the paper investigates whether a specialization of consumption by income classes emerges from simple social behaviors such as imitation, differentiation and innovation. Our purpose is to build an agent-based model of luxury consumption, which could easily be included in macro-models in order to study the impact of demand for luxury at an aggregate level.

In 1899, Thorstein Veblen analyzed consumption as a class-phenomenon, devoting part of his study to the concept of conspicuous consumption. Then, Pierre Bourdieu investigated the links between consumption and social status, emphasizing that "the world of luxury goods [...] seems predisposed to express social differences" (Bourdieu 1979, p.249). The 1857's empirical analysis of Ernst Engel about the way people allocate their budget shares among different categories of * anais.carlin@gredeg.cnrs.fr † verel@i3s.unice.fr ‡ philippe.collard@unice.fr 1 expenditures with respect to their socioeconomic class, and the subsequent concept of hierarchic preferences both revealed links between income level, wants and consumption. Finally, [START_REF] Aversi | Demand dynamics with socially evolving preferences[END_REF] investigated the main features of socially based individual preferences using the agent-based computational economics methodology. In that sense our model may be seen as a simplification of this latter. Indeed, our model requires less parameters and for this reason, it is easier to deal with, to replicate and to include in a macro-model. However, it may also be seen as a deepening of the [START_REF] Aversi | Demand dynamics with socially evolving preferences[END_REF] model because of its focus on the analysis of preference for luxury.

Methodology

We build an agent-based model of luxury consumption. In the model n agents are distributed among different income classes, denoted by c i ∈ {1..m}, from the poorest to the richest. Two types of goods are available in the economy: basic goods (B) and luxury goods (L). At each time t, each individual i spends all his income on consumption; a share θ i t ∈ [0, 1] of his income is devoted to good (L) and a share 1 -θ i t to good (B). θ i t is called consumption profile of individual i at time t. We set that three main factors are involved in the process of demand formation in a social context: imitation, differentiation and innovation. We first study analytically each process, separately, and then simultaneously. We then use simulation (run on NetLogo) to implement the whole model where all individuals act simultaneously. This first stage of the model takes place in a non-monetary framework.

Summary of the results

We first conduct a two classes analysis (i.e. m = 2) where we primarily look at the impact of each process: innovation, imitation and differentiation, taking place alone. The central role of innovation is to include inter-individual variety into the model. The innovation process is a random walk. When analyzing the imitation and differentiation process, we find that they act in a similar way on the consumption profile. The role of both imitation intensity, controlled for the size of the class, and differentiation intensity are symmetric ones. Both processes act together on the strength and the speed of evolution of the profiles. While imitation tends to bring closer the consumption profiles in an exponential way, differentiation tends to drive them exponentially away one from another. We are thus in the presence of two antagonistic dynamics, so that a fixed point can only result from an equilibrium state between the two forces.

We also study the joint effect of imitation and differentiation. We observe a pursuit racing phenomenon, in the sense that the consumption profile of the lower class is following the shape of the one of the upper class. Indeed, the pursuit racing schema is such that, by the differentiation process, class 2 profile keeps staying ahead of class 1 profile, and this despite the imitation dynamics of this latter. We find that, outside of boundary cases, the distance separating both consumption profiles converges exponentially toward a limit distance. therefore, this primary analysis confirms our hypothesis that there actually exists a differentiation of the consumption profile by income classes. Note that this limit distance is proportional to the threshold δ. Moreover, the time of convergence toward a boundary is a decreasing function of one of the imitation and differentiation intensities and is decreasing with δ.

Concerning the whole model, we are seeking some sets of parameters' values for which each class c displays a unique and specific consumption profile. By doing so, we are looking for a stable structure of consumption, which fits the fact that the consumption profile of rich people embodies luxury goods, while poorer people are only purchasing basic goods. The strength of our model leans on the fact that we do not need to enforce special characteristics on luxury goods, such as being an innovative or a traditional good. An item is defined as luxurious by its ability to exhibit social distinction. One emerging property of the simulations is the appearance of two distinct consumption profiles within the upper class in the two-classes analysis, within the middle and the upper class in the three classes analysis, and mainly within the middle class in the four classes analysis. Note that these emerging profiles are also clearly distinct from all the consumption profiles of the other classes and are persistent in time. This splitting phenomenon can therefore be interpreted as the emergence of additional consumption classes, or more precisely as the emergence of consumption sub-classes. It appears that the splitting phenomenon is related to both the relative size of each class and the intensity of the imitation and differentiation parameters. We observe that low levels of the imitation-intensity parameter tend to foster the appearance of splits, while the differentiation-intensity and the threshold parameter δ act on the distance separating the two profiles resulting from the split. Note that we never observe any persistent split within the lowest class. Interestingly, this latter class is the only one which only imitates and never differentiates itself. This observation tends to suggest that it is the differentiation process which is in action when splits occur. Indeed when the degree of imitation is low, the differentiation process has a larger influence in the determination of the consumption profile. If we now try to interpret these results, it is interesting to see that they are quite in tune with sociological evidence. Indeed, [START_REF] Warner | Social Life of a Modern Community[END_REF] propose a social categorization based on the individual perception of one's social space. According to this empirical study society appears to be divided into three main classes: the upper, middle and lower class, themselves divided in two sub-classes, the upper-upper, the lower-upper, the uppermiddle, the lower-middle, the upper-lower and the lower-lower class. Apart from the sub-division of the lower class on which we will come back to later on, our splitting phenomenon seems to fit this sociological framework quite well. Indeed, this sub-class consciousness explains both the splits that occur within the middle and the upper classes in the three classes analysis, but also the fact that, when we already deal with four classes (the middle class being defined as the second lowest class c = 2), we observe a split only within the middle class since the distinction between lower-upper and upper-upper class is already done. For what concerns the lowest class, the absence of split, due to our previous hypothesis that individuals only differentiate themselves with respect to individuals belonging to strictly lower classes, fits the idea of Bourdieu, according to which intra-class distinction does not appear within the more modest class, but only within the upper classes, and especially within the bourgeoisie, which corresponds to our middle class.

Further developments would consist in including the model in a monetary framework. This second stage of the model would connect luxury goods with high price levels. The introduction of price would also allow investigating the link between luxury consumption and income inequality.

In particular we seek to test the existence of Veblen effects and to confront our model with the stylized facts on the elasticity of demand with respect to price for each categories of commodities.

  People aim at imitating their peers and individuals from the upper classes. The imitation process follows a bayesian-updating rule. Then, during the differentiation process individuals seek to distinguish their consumption profile from individuals belonging to strictly lower classes. A threshold parameter δ defines the maximal proximity tolerated by each individual with respect to individuals from lower classes. Differentiation occurs if individuals belonging to distinct groups are distant by less than δ. Finally, what we call innovation is the production of a new consumption profile by an agent, independently from other individuals' profile. The innovation process follows a truncated normal law between 0 and 1.


Gladys Loranger-Merciris

Gaëlle Damour

Deloné-Louis Jeune

Harry Ozier- Lafontaine

Marc Dorel

Jorge Sierra

Jean-Louis Diman

Patrick Lavelle

Management practices and incidence of pests in Plantain (Musa paradisiaca AAB) crops. Consequences on the sustainability of the cropping systems

Keywords: Agroecological diagnosis, Pest control, Plant health, Soil biodiversity, Tropical agroecosystems

Plantain (Musa paradisiaca AAB) is one of the most important staple crops in the tropics, particularly in the Caribbean. Pests are the main constraint to plantain production and yield increases may be possible by improving pest management. However, there is a lack of data on these cropping systems and a need to identify new elements of improved production systems that can control pests in line with the principles of agroecological transition. In this study, we test the hypothesis that crops grown in good quality soils are less susceptible to pests. To this end, an agroecological diagnosis of the biological, physical, morphological and chemical soil conditions and the occurrence of pests, in particular plant-parasitic nematodes, was carried out on 21 plantain plots with contrasting management in Guadeloupe (French West Indies). After classifying these plots according to their management practices, from conventional to agroecological, we searched for relationships between cropping systems, soil quality and pests using a series of synthetic quality indices. Four different cropping systems were identified and compared: conventional intensive, intermediate, low-input and agroecological, according to the type and amount of agrochemicals used. Our data showed that the Chemical Index was significantly improved in the agroecological systems due to increases in pH, CEC, base saturation and total C and N content. The Morphological Index increased regularly from intensive to agroecological systems and was significantly higher in agroecological systems with higher proportions of biogenic aggregates. Soil macroinvertebrate communities were significantly influenced by production systems, with the highest values of the Macrofauna Index in low input and agroecological production systems.

We also showed that plant parasitic nematodes were better regulated in agroecological production systems. Furthermore, the agroecological system achieved similar levels of crop production to conventional systems, without the use of pesticides. In addition, these agroecological cropping systems extend the duration of plantain plots, which is an indicator of the good health of these systems.

Introduction

Plantain (Musa paradisiaca, AAB) is a major food resource in the tropics. This crop plays an important socio-economic role for developing countries in tropical and subtropical zones. Indeed, plantain provided food security and income for small growers who represent the majority of producers in the tropics [START_REF] Kahane | Agrobiodiversity for food security, health and income[END_REF]. It generates an important source of income for producing countries [START_REF] Tchango Tchango | Plantain: post-harvest operations[END_REF] and constitutes the staple diet of more than 400 million people living in many developing countries, particularly in Africa, Latin America and the Caribbean [START_REF] Lassoudière | Le bananier, un enjeu mondial[END_REF]. In 2020, plantain was grown in 54 countries and more than 43 million tons were produced worldwide (FAOSTAT, 2021). Latin America and the Caribbean produced 18% of the global production of plantain (FAOSTAT, 2021). Therefore, sustainable production of plantain in this region is crucial to provide food security and income to millions of people [START_REF] Picq | Bananas and Food Security. Les productions bananières : un enjeu économique majeur pour la sécurité alimentaire[END_REF]. In spite of this, plantain is confronted to numerous sanitary constraints, which affect its sustainability. Accordingly, production and yield gains can be increased through pest control improving, indeed the major constraint for plantain production [START_REF] Godefroid | Associations of soil type and previous crop with plant-feeding nematode communities in plantain agrosystems[END_REF].

Soil organisms provide essential ecological functions in agricultural and natural ecosystems.

These functions are the basis for the provision of ecosystem services [START_REF] Bardgett | Belowground biodiversity and ecosystem functioning[END_REF][START_REF] Adhikari | Linking soils to ecosystem services. A global review[END_REF][START_REF] Brussaard | Biodiversity and ecosystem functioning in soil: The dark side of nature and the bright side of life[END_REF]. Agriculture is one of the human activities that relies on the soil as a support, but also as a resource. However, artificializing practices of conventional agriculture have often had negative impacts on the environment, such as groundwater contamination, biodiversity loss and soil degradation [START_REF] Tamburini | Agricultural diversification promotes multiple ecosystem services without compromising yield[END_REF][START_REF] Lavelle | Soil macroinvertebrate communities: a worldwide assessment[END_REF]. These environmental disturbances result in the degradation of some of the ecosystem services provided by soils (Millennium Ecosystem Assessment, 2005).

The major underground pests in plantain cropping systems are nematodes (mainly Radopholus, Pratylenchus and Helicotylenchus) and black weevils [START_REF] Gold | Banana Weevil, Cosmopolites sordidus (Germar) (Coleoptera: Curculionidae)[END_REF][START_REF] Quénéhervé | Integrated management of banana nematodes[END_REF]. Plant parasitic nematodes attack roots, inducing plant toppling and yield reduction. The banana black weevil, Cosmopolites sordidus, attacks the plant underground corm, causing stem breakage. At the present time, these pests are mainly controlled by synthetic pesticides which has generated major side effects on environment (soil and water pollution, disappearance of domestic biodiversity and auxiliary fauna, ...) but also on human health (Millennium Ecosystem Assessment, 2005). As an example, in Guadeloupe, the use of chlordecone to control C. sordidus in dessert banana plantations has resulted in a severe long-term pollution of soils and water and has reinforced the anxiety of consumers about the quality and safety of food production [START_REF] Cabidoche | Long term pollution by chlordecone of tropical volcanic soils in the French West Indies: A simple leaching model accounts for current residues[END_REF][START_REF] Levillain | Analysis of environmental and farming factors of soil contamination by a persistent organic pollutant, chlordecone, in a banana production area of French West Indies[END_REF]. There is a critical need for alternative plantain cropping systems that better conserve environmental quality while sustaining crop yield and regulating pests [START_REF] Loranger-Merciris | Fast improvement of macrofauna communities and soil quality in plantain crops converted to agroecological practices[END_REF]. Reduction of agrochemical inputs is the target proposed, especially in Guadeloupe where ecosystems have attained high levels of contamination (Ozier-Lafontaine and Lesueur-Jannoyer, 2014).

Soil organisms but also physical and chemical properties of soil affect plant vigour. Low soil biodiversity is actually associated with greater vulnerability to pests, due to altered natural regulation of these pests [START_REF] Nielsen | Soil Biodiversity and the Environment[END_REF]. Moreover, some studies showed an unexpected positive relationship between the frequency of pesticide applications and the occurrence of pests and diseases [START_REF] Altieri | Agroecología: teoría y práctica para una agricultura sustentable[END_REF]. On the other hand, only few studies have actually reported a positive relationship between soil quality, plant health and pest regulation (e.g., [START_REF] D'hose | The positive relationship between soil quality and crop production: a case study on the effect of farm compost application[END_REF], and none to our knowledge, has been done in plantain plantations.

This study aimed at identifying elements of cropping systems that would reduce the impact of pests and diseases and, thus, decrease the need for pesticides. For this, we realized an extensive comparative diagnosis of biological, morphological, physical and chemical soil conditions and plant disease occurrence in 21 farms with contrasted management of plantain in Guadeloupe (French West Indies). Plots were classified according to their management practices. We then calculated synthetic indexes to summarize the information collected for soil quality and plant health at each farm, and we linked theses indexes to the types of management practices. Throughout this study, we tested the hypothesis that plants grown in good quality soils are less prone to pests.

Material and Methods

Site information and study design

The study was carried out in the western island (Basse-Terre) of the Guadeloupe archipelago (15°59'-16°9' N, 61°34'-61°43' W) in the Caribbean, where most plantain plantations are found. We selected a plot in twenty-one farms. Farms were selected on the basis of a preliminary survey of 100 farms to have a representation of all production systems. Plots were all located between 0 m and 150 m altitude, with annual rainfall ranging from 2 to 3 m and mean annual temperature around 23°C. Soils are Haplic Nitisols (FAO-UNESCO classification, [START_REF] Driessen | Lecture notes on the major soils of the world[END_REF] developed on volcanic ashes. They are clayey and are mainly composed of halloysite [START_REF] Clermont-Dauphin | Effects of intensive monocropping of bananas on properties of volcanic soils in the uplands of the French west indies[END_REF]. Nitisols are similar to Ferralsols, but they are younger and therefore less evolved, less desaturated and generally less acidic. The cation contents (Ca, Mg and K) are typically higher than in Ferralsols.

For each plot, cropping systems were described with seven variables: the origin of plants (plants issued for tissue culture vs. suckers); the previous crop (fallow, plantain or pineapple); the longevity of the plot (perennial vs. annual); the number of pesticide applications (0, 1 or 2 applications of Nemathorin® per year); herbicide applications (0 to 2 or 3 to 5 applications of Basta® per year); fertilization (organic manure, low input chemical -2 to 3 applications of N20-P20-K20 fertilizer per year-, high input chemical -4 to 8 applications of N20-P20-K20 fertilizer per year-); and soil tillage before plantain plantation (with or without).

All plots were planted with French cultivar.

In each plot, soil chemical, morphological and physical parameters, soil macroinvertebrate and plant-root parasitic nematode communities were assessed next to the stem of five plants selected at random at the flowering stage. Data collected within each plot were averaged to provide mean values for each of the 21 selected plots. Thus, for each measured variable within each plot values represent means of five replicates.

Soil chemical characteristics

Soil chemical properties were measured 30 cm apart from each selected plantain plant. For each plant, a sample was taken from a soil block and analysed for total soil C and N using an auto-analyser (NF ISO 13878, 1998). Available P was measured using the Olsen method [START_REF] Dabin | Sur une méthode d'analyse du phosphore dans les sols tropicaux[END_REF]. Soil mineral N (N03 -and NH4 + ) was measured by colorimetry after extraction with a 0.5M KCl solution. Exchangeable Ca 2+ , Mg 2+ , K + and Na + were determined after extraction with ammonium acetate (NF X 31108). Cation-exchange capacity (CEC) was measured with IF07-10D (NFX 31-130) method and pH-H2O and pH-KCl with NF ISO 1770, 3696 and 11464.

Soil physical and morphological characteristics

Soil physical properties were measured 50 cm apart from each selected plantain plant. Three undisturbed 100 cm 3 soil cores (50 mm in diameter × 50 mm in length) were collected in order to measure bulk density, water content, useful water reserve and water retention at two different potentials, pF 2 (-0.01 MPa) and pF 4.2 (-1.5 MPa).

A morphological assessment of soil macroaggregates and other elements (Velasquez et al., 2007b) provided an integrative measure of accumulated soil biological activity [START_REF] Lavelle | Soil ecosystem services and land use in the rapidly changing Orinoco River Basin of Colombia[END_REF]. An undisturbed soil sample taken with an 8 cm diameter  8.5 cm height cylinder was collected 50 cm apart from each selected plantain plant. In the laboratory, blocks were gently separated and the samples were air-dried then sieved at 2 mm. The soil retained by the sieve was placed on a filter paper and the various elements were separated according to their origin and shape. The elements were sorted according to five categories: i) biogenic aggregates of rounded forms, created by macroinvertebrates, ii) root aggregates, iii) physical aggregates of angular forms, produced by the physical processes of the environment (especially alternating dry and wet periods), iv) non-macro aggregated soil (soil passing through the 2 mm sieve), and v) other components (roots, stones, litter, wood…). The separated elements were put in an oven at 60 °C for 15 days and weighed (Velasquez et al., 2007b).

Macroinvertebrate communities

Soil macro-invertebrates were hand sorted from soil monoliths taken 30 cm apart from the 5 selected plants of each plot, following the ISO 23611-S methodology [START_REF] Iso | Qualité du sol -Prélèvement des invertébrés du sol -Partie 5: Prélèvement et extraction des macro-invertébrés du sol[END_REF]. A central monolith 25 cm  25 cm  20 cm was completed with two lateral 25 cm  25 cm  10 cm blocks located 1 m N and S, respectively, apart from the central point. Macro-invertebrates were classified into fifteen groups: Diplopoda, Ants, Termites, Earthworms, Coleoptera, Chilopoda, Isopoda, Dermaptera, Blattodea, Araneidae, Heteroptera, Gasteropoda, terrestrial Turbellaria, Homoptera and Orthoptera. Macroinvertebrate density at each point was taken as the sum of the three monoliths. Individuals were also classified in three functional groups: ecosystem engineers, litter transformers and predators [START_REF] Turbé | Soil biodiversity: functions, threats and tools for policy makers[END_REF].

Pests and plant characteristics

Primary roots were removed from a 25 cm  25 cm  30 cm soil block located 30 cm apart from the 5 sampled plants. Plant-feeding nematodes were extracted from an aliquot of the roots (100 g) using a centrifugal-flotation method [START_REF] Coolen | A method for the quantitative extraction of nematodes from plant tissue[END_REF].

Populations were counted in aliquots and expressed as number of nematodes per 100 grams of root fresh biomass. The main plant-parasitic species were identified and counted under a light microscope.

To assess the severity of root damage, a 25 cm  25 cm  30 cm soil block located 30 cm apart from the 5 sampled plants was removed. Roots were removed in the whole volume of soil, carefully washed and cut longitudinally. Necrosis rates on the respective external and internal surfaces of these roots were scanned independently using the WinRHIZO 2009a Software (Regent Instruments Canada, Inc.). The colour of root necrosis due to soil-borne pathogens varied from reddish to black while healthy roots were white. The necrosis index was expressed as necrosed surface area over total root surface area.

Weevil C. sordidus populations were estimated using pheromone-baited pitfall traps [START_REF] Reddy | A pheromone-based trapping system for monitoring the population of Cosmopolites sordidus (germar) (Coleoptera: Curculionidae)[END_REF]. One trap was put in the middle of each plot and was collected after 30 days (one sample by plot). Weevils were trapped in soapy water. Plant yield was estimated by bunch weight for each selected plantain plant. The average harvest time is 12 months from planting.

Data analyses

All the data sets used for statistical treatments are grouped in Supplementary Table S1. All analyses were performed with the R software (R Core Team, 2013, R 3.1.2).

A multiple correspondence analysis (MCA) was run on the 7 variables describing the cropping systems (function "dudi.pca", package "ade4", [START_REF] Chessel | The ade4 package -I-One-table methods[END_REF]. A hierarchical cluster analysis using the Ward linkage method (minimum variance criterion) and Euclidian distances was run on MCA-axes 1 and 2 scores, to analyse distances between plots and group plots according to their cropping systems. Differences in MCA-axes 1 and 2 scores among groups were further tested by one-way ANOVAs (normality was assessed with Shapiro-Wilk tests), followed by Tukey post hoc tests (function "Tukey HSD", package "stats").

Principal Component Analysis (PCA) of each of the 5 sets of data (chemical soil properties, physical soil properties, soil morphological characteristics, macroinvertebrate communities, and pests) were performed. Monte Carlo tests on the coordinates of plots on the factorial axes allowed identifying significant differences among cropping systems. These analyses were done with the ade4 library [START_REF] Chessel | The ade4 package -I-One-table methods[END_REF].

We then calculated synthetic indexes, defined as a set of statistically selected soil attributes that incorporate key functions, which for our study corresponded to: soil physical properties, soil morphological characteristics, soil chemical properties, macroinvertebrate communities, and plant pests. Synthetic indicators were designed and calculated for each point using a methodology proposed by Velasquez et al. (2007a). Briefly, for each PCA, we selected, for the two first axes, the variables with highest contributions (with a weight at least higher than 50% of the maximum contribution calculated for F1 and F2 of the PCA). Then, the values of each variable set, adjusted to a range of 0.1 to 1.0 by a homothetic transformation, were multiplied by their corresponding weight factors and summed, giving i) Physical Index, ii) Morphological Index, iii) Chemical Index, iv) Macrofauna Index, and v) Pest Index (for the calculation details, see [START_REF] Rodriguez | Agroforestry systems in the Colombian Amazon improve the provision of soil ecosystem services[END_REF].

In a last step, a PCA analysis of all plots characterized by the set of the 5 synthetic indexes generated allowed to analyse co-variations among the indicators and localize each plot in the factorial space thus defined. Finally, a Monte Carlo test on plot coordinates allowed to test significant differences among cropping systems. These analyses were done with the ade4 library [START_REF] Chessel | The ade4 package -I-One-table methods[END_REF].

Results

Cropping systems

The first two axes of the MCA run on the seven variables related to the cropping system (Supplementary Table S1) explained 46.9% of the total variance of the observations (data not shown). The first axis alone explained 33.6% of this variance and opposed sites with highly intensive practices (tillage, 5 to 6 applications of chemical fertilizer, repeated pesticide applications and annual cropping systems with local seeds), to sites with agroecological practices (organic fertilization, zero application of chemical fertilizer, limited use of herbicides and planting with vitroplant seeds in perennial crops), (Supplementary Fig. S1). Plots were divided in 4 clearly separated groups (Fig. 1; Supplementary Fig. S2).

Groups 1, 2 and 3 are annual plots planted with suckers separated along axis 2 of the analysis according to their relative degree of intensification that increases from 3 to 1 (Table 1). Group 1, called "intensive" cropping system, is characterized by high levels of pesticides (2 applications during the crop cycle), herbicides (3 to 5 applications) and mineral fertilizers use (4 to 8 applications). Group 2, called "intermediate" cropping system, is characterized by an intensive use of herbicides (3 to 5 applications during the plantation cycle), moderate use of pesticides (1 application at the plantation) and moderate use of mineral fertilizers (2 to 3 applications during the plantation cycle). Group 3, referred as "low input" cropping system, is characterized by moderate levels of pesticides and mineral fertilizer applications (respectively 1, and 2 to 3 applications during the plantation cycle) and no herbicides. Group 4 (called "agroecological" system) is separated from the other 3 along axis 1 of the MCA (Fig. 1). It is comprised of perennial plots, planted with seeds issued from tissue culture conducted with organic fertilization (manure). Pesticides were not used and 0 to 2 applications of herbicides had occurred in the last 12 months. In this group, unlike the others, there was no tillage.

Soil chemical characteristics

PCA analysis performed on the matrix of chemical characteristics significantly separated the four production systems (p<0.01). Factor 1 (48.61% explained variance) separated the agroecological system, with better chemical conditions (Fig. 2a,2b). Analyses of variance showed that soil pH-H2O and pH-KCl were significantly higher in agroecological systems compared to other systems. Soil N-NH4, N-NO3, total N, total C, CEC, K, Mg, P and Na content did not significantly differ among the four cropping systems. Despite this, the highest total C and N contents were observed in the agroecological system. Moreover, Ca content was significantly higher in agroecological systems compared to intensive systems. Therefore, the Chemical Index was significantly enhanced in agroecological systems compared to other systems, due to high pH, CEC, base saturation and total C and N contents (Fig. 2c).

Soil physical and morphological characteristics

PCA analysis performed on the matrix of physical characteristics did not separate the four production systems (p>0.01). Soils generally had rather low bulk densities around 0.9, especially in the conventional cropping systems (Supplementary Table S1). Soil in the agroecological systems were comparatively more compact, although with bulk densities of less than 1.2, and a remarkably low value of water retention at pF 4.2 which indicate of a very low proportion of pores of ≤ 0.15 µm. The Physical Index did not exhibit significant differences among cropping systems.

PCA analysis performed on the matrix of morphological characteristics significantly separated the agroecological system from the others (p<0.01, Fig. 3a,3b). On average, a rather large proportion (49.4%) of the soil volume was macro-aggregated with large variations from 14.1 to 81.3% depending on treatments and sites. Biogenic aggregates comprised 31.4% of the total soil mass and physical aggregates 5.9%. There was a clear trend for soils to have a larger proportion of biogenic aggregates along the sequence from high inputs (23.8% of the total soil weight on average), to intermediate and low and agroecological cropping systems (46.1% of the total soil weight). Agroecological systems had significantly higher proportions of biogenic aggregates than the other three, and lower proportions of physical aggregates and non-aggregated soils than intensive and low input systems. Subsequently, the non-aggregated soil volume was higher in other systems (54.8%) compared to agroecological systems (33.9%). The proportion of physical aggregates was also higher in intensive systems (8.7%) as compared to low input systems (4.7%). The Morphological Index regularly increased from intensive to agroecological systems; it was significantly higher in agroecological than in conventional systems (Fig. 3c).

Soil macroinvertebrate communities

Soil macroinvertebrate communities varied among cropping systems, with respective densities ranging from a minimum value of 357 ± 48 ind.m -2 to a maximum of 886 ± 112 ind.m -2 and an overall average of 618 ± 135 ind.m -2 . Regardless the cropping system type, ants were the most abundant group (26%) followed by earthworms (25%), diplopods (23%) and termites (7%). Communities were significantly impacted by production systems (global PCA randtest p<0.01, Fig. 4a,4b) with highest values of the Macrofauna Index in low input

(3) and agroecological (4) production systems (Fig. 4c). Taxonomic richness did not vary significantly among cropping systems.

Pests and plant production

PCA analysis performed on the matrix of pests did not separate the four production systems (p>0.01). However, plant parasitic nematode communities varied among cropping systems, ranging from 60809 ± 39413 to 180166 ± 53667 ind.100 g roots -1 Total population of plant parasitic nematodes was not significantly different among cropping systems. Nevertheless, Radopholus similis population was significantly reduced in agroecological and intensive systems and was higher in intermediate and low input systems (Table 2). Helicotylenchus multicinctus population, on the contrary, was significantly higher in agroecological cropping systems (Table 2).

The severity of root damage induced by parasites also varied within cropping systems. The root necrosis index was lower, although not significantly, in intensive and agroecological systems (respectively, 37 ± 3 % and 36 ± 5 %) compared to intermediate and low input systems (respectively, 54 ± 4 % and 43 ± 3 %).

Black weevil captures varied among cropping systems, ranging from 16 ± 2 ind.month -1 to 64 ± 3 ind.month -1 .The abundance of black weevils was not significantly different among cropping systems.

The Pest Index exhibited rather important variations among cropping systems, with highest values in the intermediate system and lowest on average in the intensive system despite a large variability.

Plant yield (expressed by bunch weight) averaged 15.6 ± 0.8 kg in low input systems, 16.6 ± 1.3 kg in intermediate systems, 17.2 ± 1.5 kg in agroecological systems and 19.7 ± 0.8 kg in intensive systems. Differences, however, were not significant.

Covariations

PCA analysis performed on the matrix of farm general Chemical, Physical, Morphological, Macrofauna and Pest indexes significantly separated the four production systems (p<0.01).

Factor 1 (40% explained variance) separated the agroecological system, with better morphological and chemical conditions from the other 3 ones that had higher values of the physical indicator. Factor 2 showed a clear opposition between sites with high macrofaunal index (mostly found in low input (3) and agroecological (4) systems) and sites of production systems 1 and 2 with higher values of the Pest Index (Fig. 5a,5b).

Discussion

Cropping systems and soil health

Around the world, plantain is mostly cultivated in traditional, extensive systems, and consumed locally [START_REF] Tchango Tchango | Plantain: post-harvest operations[END_REF]. Plantain cropping systems in Guadeloupe are comparatively more intensive with high inputs of mineral fertilizers, pesticides and herbicides. Only 14% of the plots have been classified as agroecological and 28% as low input, characterized by the predominant use of organic fertilizers, and low applications of pesticides and herbicides. In Guadeloupe, plantain plantations are mainly found in the banana crop region. They coexist with intensive banana dessert monocultures that use high amounts of pesticide and mineral fertilizers [START_REF] Sierra | Observed and predicted changes in soil carbon stocks under export and diversified agriculture in the Caribbean. The case study of Guadeloupe[END_REF]. Most plantain producers indeed use similar farming practices as those used in banana dessert crops with 8-12 applications of chemical fertilizer during the growth season in monocultures [START_REF] Sierra | Observed and predicted changes in soil carbon stocks under export and diversified agriculture in the Caribbean. The case study of Guadeloupe[END_REF].

Our study showed that soil morphological, macrofauna and chemical characteristics were improved under agroecological plantain cropping systems. We observed in these systems a significant accumulation of biogenic aggregates. The high concentration of biogenic aggregates in agroecological plots can be partly explained by the reduction of tillage, which results in improved stability of the soil structure within a few years [START_REF] Haynes | Labile organic matter as an indicator of organic matter quality in arable and pastoral soils in New Zealand[END_REF]. This can be reinforced by the maintenance of a permanent litter of plant residues on the soil surface.

A previous study carried out on Guadeloupean andosols showed that deep tillage caused a decrease in soil organic matter content [START_REF] Clermont-Dauphin | Effects of intensive monocropping of bananas on properties of volcanic soils in the uplands of the French west indies[END_REF]. Indeed, our results suggest that in agroecological systems, favourable soil structure could be maintained by macrofaunal activity, as previously reported by [START_REF] Velasquez | Soil macrofauna as an indicator for evaluating soil based ecosystem services in agricultural landscapes[END_REF] Macrofauna characteristics were also improved in low input and agroecological cropping systems. In a recent study, [START_REF] Lavelle | Soil macroinvertebrate communities: a worldwide assessment[END_REF] showed that the worldwide abundance of soil macroinvertebrate communities drastically decreases in conventional cropping systems. In our study, the increase in macrofauna abundance within the ecological friendly cropping systems may be due to moderate levels of pesticides, herbicides and mineral fertilizer applications [START_REF] Sánchez-Moreno | Multivariate analysis of toxicological and environmental properties of soil nematicides[END_REF][START_REF] Niedobová | Effects of glufosinate-ammonium herbicide and pod sealant on spider Pardosa agrestis[END_REF][START_REF] Loranger-Merciris | Fast improvement of macrofauna communities and soil quality in plantain crops converted to agroecological practices[END_REF]Zulu et al. 2022).

In addition, our results showed that the Chemical Index was significantly enhanced in agroecological systems compared to other systems, due to high pH, CEC, base saturation and C and N contents. [START_REF] Sihi | Evaluation of soil health in organic vs. conventional farming of basmati rice in North India[END_REF] found that long-term application of organic residues in basmati rice crop improved physical, chemical, and biological indicators of soil health.

These authors observed that organic residues addition favoured biologically available nutrients in organic systems, and stimulated soil microorganisms that produce enzymes responsible for the conversion of unavailable nutrients to plant available forms. Thus, this study indicates, as ours, that agroecological farming practices may promote biotic and abiotic interactions in the soil which may support sustainable agriculture.

Cropping systems, pests and plant health

The plant-parasitic nematode communities sampled in our investigation (H. multicinctus, Hoplolaimus sp, Meloidogyne sp, Pratylenchus coffeae, R. similis, Rotylenchus reniformis) was similar to previous studies carried out on plantain cropping systems in Côte d'Ivoire [START_REF] Adiko | Evolution de la nématofaune du bananier plantain (Musa AAB) en Côte d'Ivoire[END_REF], Rwanda [START_REF] Sebasigari | les charançons et les nématodes du bananier dans les pays C[END_REF], and Cuba [START_REF] Rodriguez | Aspectos cocio económicos del cultivo del plátano en Colombia[END_REF]. An interesting result of our study was that R. similis, the most destructive species on Musa sp., was absent in the agroecological perennial cropping systems. This can be explained by the fact that these plots were planted with vitroplants and thus were healthy at the beginning. This result highlights that particular attention should be paid to the sanitary state of the plant material used during plantain plantation for ensuring plantation durability. In this sense, it has been already observed that the use of healthy plants, obtained through the PIF technique [START_REF] Lefranc | Macropropagation as an innovative technology: Lessons and observations from projects in Cameroon[END_REF][START_REF] Sadom | Comparaison de l'efficacité de deux méthodes de multiplication rapide de plants de bananier à partir de l'étude des caractéristiques agronomiques d'un hybride de bananier plantain (Musa spp.)[END_REF], significantly reduced pest infestations at the beginning of the plantation [START_REF] Loranger-Merciris | Fast improvement of macrofauna communities and soil quality in plantain crops converted to agroecological practices[END_REF]. In contrast, the higher density of H. multicinctus observed in agroecological perennial banana plantations showed that these ectoparasite nematodes are more persistent in more organic cropping systems. This same result was obtained by [START_REF] Pierrot | Déterminants de la production de plantain dans le sud-ouest du Cameroun : intérêt d'un observatoire[END_REF] in plantain plantations in

Cameroon. An important fact is that H. multicinctus is not the most damaging on banana (effect on bunch weight) as shown by [START_REF] Moens | Reproduction and pathogenicity of Helicotylenchus multicinctus, Meloidogyne incognita and Pratylenchus coffeae, and their interaction with Radopholus similis on Musa[END_REF]. The dynamics of plantain roots can also alter the ratio of H. multicinctus to R. similis [START_REF] Quénéhervé | Effect of banana crop mixtures on the plantfeeding nematode community[END_REF]. Indeed, the proportion of old roots (potential resource for H. multicinctus) increased in agroecological perennial plantain plots compared to that of young roots (potential resource for R. similis), [START_REF] Tixier | Modelling pest dynamics of new crop cultivars: the FB920 banana with the Helicotylenchus multicinctus -Radopholus similis nematode complex in Martinique[END_REF].

We also showed that pest occurrence, in particular plant-parasitic nematodes communities, was minimal in agroecological cropping systems. Our study revealed that cropping systems may impact populations of the plant parasitic nematode R. similis and subsequent effects on plant health. In fact, R. similis population was significantly reduced in agroecological systems characterized by a low use of agrochemicals, which indicates that organic amendments can regulate plant parasitic nematode populations in plantain crops, as already observed for banana [START_REF] Pattison | Soil quality improvement and nematode management on banana farms in Australia[END_REF][START_REF] Tabarant | Effects of four organic amendments on banana parasitic nematodes and soil nematode communities[END_REF]. These authors suggested that the regulation of these plant parasitic nematodes was associated with an increase in the diversity and functional composition of soil nematode communities. In particular, the supply of organic matter boosted the populations of microbivorous nematodes, particularly bacterivores [START_REF] Ferris | Nematode indicators of organic enrichment[END_REF] and of carnivorous nematode populations [START_REF] Tabarant | Effects of four organic amendments on banana parasitic nematodes and soil nematode communities[END_REF] and not of plant parasitic nematodes. The equilibrium state of communities also limits pest development [START_REF] Altieri | The ecological role of biodiversity in agroecosystems[END_REF]. As an example, in Senegal, millet is attacked by two main nematode species, Tylenchorhynchus gladiolatus and Scutellonema cavenessi. These two species represent more than 95% of the nematode community of these agroecosystems.

When millet fields were left fallow, the number of these two nematodes decreased, while there was a marked increase in the numbers of other nematodes. Interestingly, the level of damage caused by the two main pest nematode species was reduced when they were associated with other plant-feeding nematodes such as one of those recovering during fallow periods, Helicotylenchus dihystera [START_REF] Cadet | Role of plant parasitic nematodes and abiotic soil factors in growth heterogeneity of sugarcane on a sandy soil in South Africa[END_REF]. R. similis populations were also reduced in intensive cropping systems presumably as a result of the application of chemical nematicides which protected plants, as least temporarily.

In fact, chemical nematicides only allow a partial and short-term control of nematodes [START_REF] Moens | Enhanced biodegradation of nematicides after repetitive applications and its effect on root and yield parameters in commercial banana plantations[END_REF]. This justifies, the local adoption of short-time duration of intensive plantain plots; e.g., 1 cycle in 85% of the cases, and 2 or more cycles in 15%. Nevertheless, these chemical nematicides can be damaging to human health and the environment [START_REF] Matthews | Pesticides: health, safety, and the environment[END_REF]. Overall, our study highlights the fact that the agroecological system has the same performance in terms of pest control as the conventional system, but without the application of pesticides.

Conclusion

Our study clearly showed that well managed agroecological systems can provide the same level of production and plant protection than intensive conventional systems, while improving soil quality and pest regulation. As a matter of fact, our data highlighted that there was a significant improvement in Chemical, Morphological and Macrofauna indexes in agroecological plantain cropping systems as well as a better regulation of plant parasitic nematodes. Moreover, these agroecological cropping systems lengthen the duration of plantain plots (in our study, agroecological systems are cultivated for the duration of 3 to 10 cycles). This longevity is an indicator of the good health of the agroecological cropping systems. Furthermore, in these agroecological systems, production costs may be significantly reduced by eliminating costly replanting and equally costly pesticide applications, while maintaining yield levels that are at least equivalent to conventional intensive cropping systems, characterized by a strong use of agrochemicals. However, the problems faced by plantain production, particularly a high level of telluric parasitism, led the majority of farmers to promote intensive cultural practices that aim at reducing pests without restoring soil health and the sustainability of the cropping systems. Yearly plant renewal and a massive use of agrochemicals in conventional systems enable to compensate for the lower general quality of the soil. Fast transitions from conventional to agroecological systems are possible. However, the main condition for rapid adoption of these eco-friendly cropping systems is the implementation of appropriate public policies. "intensive", G2: "intermediate", G3: "low input" and G4: "agroecological" cropping system.

Figure captions Fig. 1 .

 captions1 Figure captions

Fig. 2 .

 2 Fig. 2. Projection of (a) plots and groups and (b) variables in the factorial plane F1F2 of a

Fig. 3 .

 3 Fig. 3. Projection of (a) plots and groups and (b) variables in the factorial plane F1F2 of a

Fig. 4 .

 4 Fig. 4. Projection of (a) plots and groups and (b) variables in the factorial plane F1F2 of a PCA analysis of a matrix of macroinvertebrates communities in plantain cropping systems.

Fig. 5 .

 5 Fig. 5. Projection of (a) plots and (b) variables grouped into plantain cropping systems of

  Fig. 1

Table 2

 2 Plant parasitic nematode communities in 4 types of plantain cropping systems in Guadeloupe. Values are means of 5 samples (± SE). For each nematode species, means with different upper script letter are significantly different based on HSD test (P ≤ 0.05).

		Intensive	Intermediate	Low Input	Agroecological	P
	R. similis	2093 ± 2017 a 18294 ± 12630 ab 39758 ± 16102 b	0 ± 0 a	0.003
	P. coffeae	116769 ± 34910 a 137516 ± 80124 a 85987 ± 29217 a 37527 ± 32568 a 0.621
	H. multicinctus	320 ± 255 ab	21690 ± 21603 ab	19 ± 19 a	20776 ± 16780 b 0.031
	Meloidogyne sp.	3171 ± 1964 a	2666 ± 1611 a	1378 ± 1105 a	2471 ± 1444 a 0.556
	R. reniformis	80 ± 46 a	0 ± 0 a	135 ± 85 a	18 ± 10 a	0.271
	Hoplolaimus sp.	70 ± 53 a	0 ± 0 a	0 ± 0 a	17 ± 12 a	0.120

Acknowledgements

The authors thank all the farmers involved in this study. The authors also thank Fred Burner, Chantal Fléreau, Danielle Ramaël and Frédérique Razan (INRAE Antilles Guyane) for their technical help. The authors would also like to thank Lucienne Desfontaines (INRAE Antilles Guyane) for soil chemical analyses. This work was supported financially by the French Ministry of Ecology and Sustainable Development, France (ECOPHYTO program, project ALTERBIO).

Table 1

Description of the four groups of plantain plots identified in the hierarchical analysis and described as different cropping systems. Number into brackets following pesticide, herbicide and fertilizer levels are the mean number of applications during the plantation cycle observed in the plots of these groups. applications of Nemathorin® at the plantation and 6 months after the plantation. § application of Basta ® (active ingredient Glufosinate) # 1 application of Diammonium phosphate (DAP) at the plantation, then application of a complete N20-P20-K20 fertilizer