oaie Tudor-Bogdan Airimit

Christophe Collewet

Improving robust stability by increasing the number of controlled degrees of freedom

In this article we demonstrate the necessity of increasing the number of controlled degrees of freedom when controlling an uncertain system. We show this on an infinite dimensional system representative of reaction-diffusion problems. However, our final goal is concerned with fluid flow control. In both cases, fluid flow control or reaction-diffusion problems, a common approach to control is based on a spatial discretization followed by some linear control design. Unfortunately, this approach is sensitive to uncertainties introduced during linearization of the nonlinear system. In this paper, we show that robust stability is dependent on the number of controlled degrees of freedom and performance of a closed loop system can be improved by considering this aspect.

Notation

H l (0, 1) with l ∈ N * denotes the Sobolev space W l,2 (0, 1). H 0 (0, 1) is also denoted as L 2 (0, 1). We use the subscript notation for partial derivative, for example ξ x (x, y) represents the derivative of the dependent variable ξ(x, y) with respect to the independent variable x, ∂ξ /∂x, and similarly ξ xy (x, y) = ∂ 2 ξ /∂x∂y. ξ i,j denotes the value of the dependent variable in a point (i∆x, j∆y) on a grid of discretized position, where ∆x and ∆y represent the discretization steps in each of the two independent variables.

I. Introduction

The objective of this article is to point out some of the difficulties related to the control of systems described by partial differential equations (PDEs) as it is the case in flow control since a flow is governed by the Navier-Stokes equations. These represent systems where the state evolves continuously not only in time but also in space, in contrast to ordinary differential equations (ODE) which describe the evolution only in time.

In the early years of automatic control, finite dimensional systems described by ODE where studied due to the fact that only a finite number of actuators, sensors and states representing different variables of the system (velocity, pressure, temperature, humidity, etc.) at predefined locations where considered. While this is still appropriate for a large number of situations encountered in practice, there are still many physical, mechanical or chemical systems where it is not the case and the evolution in space has to be taken into account too.

In this paper, the position is denoted by x, where x ∈ R n with n = 1 (but generally n ∈ {1, 2, 3} depending on the size of the space) and the state of the system by ξ(x, t) which depends on both position x and the time t. Thus ξ(x, t) can be considered as an element of a functional space of infinite dimension. The term distributed parameters system, sometimes used in reference with systems described by PDEs, comes from the fact that the state is considered through the entire space, thus in an infinite number of positions. This is in contrast to ODE systems (called also lumped parameter systems) where the state evolution is considered only in a finite number of points from the domain.

Various discretization methods have been proposed, such as the integral, finite-difference, finite-volume, finite-element, or spectral methods. [START_REF] Kuzmin | A guide to Numerical methods for transport equations[END_REF]2,[START_REF] Rao | The finite element method in engineering[END_REF][START_REF] Tannehill | Computational fluid mechanics and heat transfer[END_REF] Control using discretized linearized models of infinite dimensional flow problems has also been proposed [START_REF] Aamo | Tutorial on Feedback Control of Flows, Part I: Stabilization of Fluid Flows in Channels and Pipes[END_REF][START_REF] Joshi | A system theory approach to the feedback stabilization of infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF][START_REF] Mckernan | A Linear State-Space Representation of Plane Poiseuille Flow for Control Design-A tutorial[END_REF][START_REF] Mckernan | Control of plane Poiseuille flow: a theoretical and computational investigation[END_REF] (see also [START_REF] Baramov | H∞ Control of Nonperiodic Two-dimensional Channel Flow[END_REF] where a finite-difference approach is used). However, in this paper we want to cope with the problem of robust stability. In fact, linearization and discretization introduces uncertainties and as soon as the distance between the state of the system and the equilibrium used for liniarization exceeds some robustness margin, the closed loop system becomes unstable. In this paper we demonstrate that robust stability can be greatly improved by increasing the number of controlled degrees of freedom, i.e. the number of independent actuators. This will be very helpful in the more complicated case of flow control. [START_REF] Dao | Minimization de l'énergie cinétique transitoire dans l'écoulement 3D plan de Poiseuille commandé[END_REF][START_REF] Dao | Fluid flow control by visual servoing[END_REF] In the next sections of this article, we describe the classical approach to controlling infinite dimensional systems (see [START_REF] Baramov | H∞ Control of Nonperiodic Two-dimensional Channel Flow[END_REF] ), which is spatial discretization using a carefully chosen set of grid points. A basic method of approximation is presented in section II and applied in section III to the heat equation with constant diffusivity to point out some of the difficulties that come from using a finite number of actuators in infinite dimensional systems. Then, in section V, by adding actuators it is shown the the robustness of a linear controller is improved.

II. Finite-difference approximation method

Most of the systems encountered in nature can be described by partial differential equations (PDEs). Therefore, in an attempt to control such systems, the need to better understand them arises. As such, numerical analysis has seen a lot of development in the last decades as the computational power of personal computers increased. Numerical analysis is based upon discretisation of PDEs and various techniques have been described.

In this section, a brief review of the finite-difference method (FDM) is presented (for a more in depth description see [START_REF] Anderson | Computational Fluid Mechanics and Heat Transfer[END_REF] ). This method is based on approximating derivatives by differences, resulting in an algebraic representation of the PDEs. The first step is to replace the continuous domain by a finite-difference mesh. For a two-dimensional (2-D) problem, ξ(x, y) is replaced by ξ(i∆x, j∆y) (written also ξ i,j ) for given distance increments ∆x and ∆y.

Recall next the definition of the derivative of a function ξ(x, y) at x = x 0 , y = y 0 :

∂ξ ∂x = lim ∆x→0 ξ(x 0 + ∆x, y 0 ) -ξ(x 0 , y 0 ) ∆x . ( 1
)
Another way to find this equation is by considering first a Taylor-series expansion for ξ(x 0 + ∆x, y 0 ) about ξ(x 0 , y 0 ) which is given by

ξ(x 0 + ∆x, y 0 ) = ξ(x 0 , y 0 ) + ∂ξ ∂x 0 ∆x + ∂ 2 ξ ∂x 2 0 (∆x) 2 2! + . . . ( 2
)
or in the (i, j) notation

ξ i+1,j = ξ i,j + ∂ξ ∂x i,j ∆x + ∂ 2 ξ ∂x 2 i,j (∆x) 2 2! + . . . (3)
which gives

∂ξ ∂x i,j = ξ i+1,j -ξ i,j ∆x + O(∆x), (4)
where (ξi+1,j -ξi,j ) /∆x represents the "forward" finite-difference approximation of the derivative and O(∆x) is the truncation error between the derivative and the proposed representation. In a similar manner, it is also possible to find a "backward" formula

∂ξ ∂x i,j = ξ i,j -ξ i-1,j ∆x + O(∆x) (5)
and a "central" formula

∂ξ ∂x i,j = ξ i+1,j -ξ i-1,j 2∆x + O(∆x) 2 . ( 6
)
As for the second derivative, this can be evaluated using the central difference formula

∂ 2 ξ ∂x 2 i,j = ξ i+1,j -2ξ i,j + ξ i-1,j (∆x) 2 + O(∆x) 2 . ( 7
)

III. Heat equation

Consider the following diffusion equation, described by the initial boundary value problem (IBVP):

ξ t (x, t) = θξ xx (x, t) + b T (x)u(t), x ∈ (0, 1) , t > 0, (8)
with boundary conditions (BCs)

ξ x (0, t) = ξ x (1, t) = 0, (9)
and initial condition (IC) ξ(x, 0) = ξ 0 (x) (10) where ξ 0 ∈ H 1 (0, 1) represents the initial system state (temperature), u(t) is a control signal, and θ represents the constant diffusivity. This is also known as the heat equation with constant thermal diffusivity.

The input influences the PDE (8) only at certain locations specified by the input shaping function

b T (x) = [b 1 (x) b 2 (x) b 3 (x) b 4 (x)] (11)
with

b i+1 (x) = 1, 0.175 + 0.2i ≤ x ≤ 0.225 + 0.2i 0, otherwise, (12)
∀i ∈ {0, 1, 2, 3}. The output is described by the equation

y(t) = 1 0 ξ(x, t)c(x)dx, t > 0. ( 13
)
where c(x) is given by

c T (x) = [c 1 (x) c 2 (x) . . . c 400 (x)] (14)
with

c i+1 (x) = 1, i 400 < x < i+1 400 0, otherwise (15)
The eigenvalues of the PDE (8) can be computed as

λ i = -θ(i -1) 2 π 2 , i ≥ 1. ( 16
)
The system has one zero eigenvalue and an infinity of negative ones. The numerical simulation of the IBVP can be done by using the method of lines. The semi-discretization of the PDE is done on the spatial domain by dividing the interval on which x is defined into 399 equal subintervals, obtaining thus N = 400 grid points.

The method of lines is used in initial boundary value problems and it consists in discretizing the system along all minus one of the variables. This results in a system of ordinary differential equations (ODE) which can be solved using various existing methods (see for example the Euler method and the various alternatives of the Runge-Kutta method described in [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF] ) efficiently implemented in simulation software such as MATLAB, Maple, Mathematica, or Scilab. After discretization and taking into account the boundary conditions, a lumped parameter system described by a state space representation of the form

χ(t) = θ ∆x 2 A L χ(t) + B L u(t), ( 17a
)
y L (t) = C L χ(t), ( 17b
)
is obtained, where χ(t) ∈ R N represents the discretized system state,

A L ∈ R N,N A L =               -2 2 0 0 0 . . . 0 1 -2 1 0 0 . . . 0 0 1 -2 1 0 . . . 0 . . . . . . . . . . . . . . . 0 . . . 0 1 -2 1 0 0 . . . 0 0 1 -2 1 0 . . . 0 0 0 2 -2               (18)
is the state matrix, B L ∈ R N,4 , and C L ∈ R N,N have zero values everywhere exception of a few elements of value 1 at grid positions corresponding to eqs. ( 12) and (15). Note that with the given output shaping function, C L = I N and therefore y L (t) = χ(t). The diffusivity constant has been chosen as θ = 0.1 and the spatial discretization step is ∆x = 1 N -1 . Using the continuous time system described by eq. ( 17), it is possible to design an LQR controller 14 that minimizes a quadratic cost function of the form

J = ∞ 0 (χ T Qχ + u T Ru)dt. ( 19
)
In the last equation, Q and R are state and input weighting matrices chosen in order to satisfy certain design restrictions. The solution of the state feedback gain matrix is given by

K = R -1 B T L P, ( 20
)
where P is found by solving a continuous time algebraic Riccati equation

A T L P + P A L -P B L R -1 B T L P + Q = 0. ( 21
)
Figure 1 shows the evolution of the temperature in the 1-D heat conductor over a 2 sec interval. This varies from an initial state ξ 0 (x) to a final one. The reference profile is marked as desired temperature profile in figure 2. In the same graphic, the obtained temperature profile is shown shown for comparison. One can observe that very good results are obtained at the point where collocated control and measurement exist while in the rest of the domain the results are not satisfactory. The main reason for this is the small number of point where measurement and control is done. Of course one can increase this number, but there will always be intermediate points where it is difficult to know exactly what the evolution of the system is.

IV. Nonlinear reaction-diffusion system

Consider now a reaction-diffusion system with nonlinear reaction given by

F : L 2 (0, 1) → L 2 (0, 1), F (ξ) = 10ξ(1 -10ξ 2 ). ( 22
)
The IBVP for this system is obtained by modifying eq. ( 8) of the IBVP described in Section III with for x ∈ (0, 1) , t > 0. The lumped parameters systems obtained after discretisation in the spatial coordinate x is given by

ξ t (x, t) = θξ xx (x, t) + F (ξ(x, t)) + b T (x)u(t), (23)
χ(t) = f (χ(t)) + B L u (24a) y L (t) = C L χ(t), (24b)
where

f (χ(t)) = θ ∆x 2 A L χ(t) + F L (χ(t)), (25)
with χ(t) representing the discretized system's state, A L , B L , and C L given previously, and F L defined as

F L : R N → R N , F L (χ(t)) = 10χ(t)(1 -10χ(t) 2 ). (26)
One approach to control system (24) is to linearize it and then build a controller based on the linear model. Linearizing around the zero temperature yields the state space system

χ(t) = θ ∆x 2 A L + 10I N χ(t) + B L u(t) (27a) y L (t) = C L χ(t) (27b)
which is unstable having four positive eigenvalues. Nevertheless, an LQR controller can be designed to stabilize the linearized system. Furthermore, a simulation run on the nonlinear system (24), shows that this is also stable, even though the set point value is not the one used for linearization. In figure 3, the temperature profiles obtained in stationary regime with the nonlinear system in closed loop with the LQR controller is compared to the desired profile. It can be observed that the tracking error is quite important in both cases and the controller doesn't manage to follow the desired profile anymore. Remark: if one would consider a nonlinear term given by +100ξ 3 instead of -100ξ 3 the closed loop system with the LQR controller is no longer stable and a different design technique has to be considered.

V. Increasing the number of controlled degrees of freedom

In this section, we will show using simulation that it is possible to control even the unstable case with nonlinear term given by +100ξ 3 using a linear LQR controller. The only change is in the number of controlled degrees of freedom. In other words, instead of having a length 4 input vector, we will use a length 10 input vector implying thus a constructive change on the physical actuators combined with an increase in the number of actuators. To be more precise, the input shaping function (11) and ( 12

for i ∈ {0, 1, . . . , 9}. The output shaping function remains unchanged.

With the increased number of actuators, the robustness of the closed loop system is improved. In figure 4, the desired profile and the final temperature profiles obtained in closed loop after the transitory period are shown. For the nonlinear system with negative nonlinearity "-", the stationary profile is very close to the desired one. Also for the nonlinear system with positive nonlinearity "+" the closed loop remains stable and the final profile is close to the desired one. In figure 5, the evolution of the head along the bar is shown.

VI. Concluding remarks

In this article, we have shown the weaknesses of conventional control, which is appropriate for lumped parameters systems, when applied to distributed parameters systems, as it is the case in flow control. We have shown the basic ideas used when trying to apply conventional control to infinite dimensional systems. The difficulties of applying linear control theory to distributed systems are highlighted. However, we demonstrate by a simulated example the practical advantages of increasing the number of controlled degrees of freedom in such systems, showing that increasing the number of actuators, when possible, can have a beneficial effect on closed loop stability. As already mentioned, this result will be of a great interest in closed-loop flow control. Recall that few degrees of freedom are used in practice, typically one or two in the case of the plane Poiseuille flow. [START_REF] Aamo | Tutorial on Feedback Control of Flows, Part I: Stabilization of Fluid Flows in Channels and Pipes[END_REF][START_REF] Joshi | A system theory approach to the feedback stabilization of infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF][START_REF] Mckernan | A Linear State-Space Representation of Plane Poiseuille Flow for Control Design-A tutorial[END_REF] Future research will deal with quantifying the robustness gain obtained by the increase in the number of controlled degrees of freedom and propose judicious actuator placement solutions.

Figura 1 .

 1 Figura 1. Temperature evolution of the controlled heat conductor.

Figura 2 .

 2 Figura 2. Comparison of desired and obtained temperature profiles (after 2 seconds).

Figura 3 .

 3 Figura 3. Comparison of desired and obtained temperature profiles (after 2 seconds).

  . sys. "-" meas. nonlin. sys. "+"

Figura 4 .

 4 Figura 4. Comparison of desired and obtained temperature profiles for the nonlinear system (after 2 seconds).

Figura 5 .

 5 Figura 5. Temperature evolution of the controlled nonlinear heat conductor with positive nonlinear term "+".

of 7 American Institute of Aeronautics and Astronautics

of 7 American Institute of Aeronautics and Astronautics

of 7 American Institute of Aeronautics and Astronautics

of

American Institute of Aeronautics and Astronautics


Professeur, sociologue Yankel Fijalkow

Yaneira Wilson

l'exigence journalistique Logement

Le gouvernement vient de proposer des mesures pour relancer la construction, la rénovation et l'accession à la propriété. Mais ces mesures prendront-elles en compte les pollutions que présentent certains habitats ? En effet on considère souvent -à tort -les logements comme des abris à l'écart des toxicités de la vie urbaine. Celles-ci sont par ailleurs élevées. Selon le The _Lancet_Planetary Health du 18 mai 2022, neuf millions de personnes meurent ainsi chaque année dans le monde à cause de la pollution, de l'air, de l'eau ou des sols soit trois fois plus que les morts cumulés des suites du sida, de la tuberculose et du paludisme. En Île-de-France, le nombre de décès évitable est estimé à 8000 par an.

Les Olympiades, Paris 13e, un des sites d'étude du programme SAPHIR (2023). SAPHIR, Author provided https://theconversation.com/logement-un-enjeu-de-sante-publique-207041 2/4

Or, l'Organisation mondiale de la Santé (OMS) a montré le rôle primordial de la qualité du logement dans la prise en charge de ces questions. Ces éléments correspondent à la notion de santé dans le sens de l'OMS (un sentiment de bien-être physique et mental) on peut étudier l'effet de la représentation par les habitants de leur habitat sur leur santé. La qualité des habitats peut-elle être évaluée à partir de la notion de santé ?

Impliquer les habitants

La qualité des habitats peut-elle être évaluée à partir de la notion de santé ?

4

 4 enjeu de santé publique https://theconversation.com/logement-un-enjeu-de-sante-publique-207041 3/Santé et confort de l'habitat ont une longue histoire Le XIX e siècle correspond à l'entrée du logement comme question politique. En France, la loi sur l'habitat insalubre, adoptée en 1850, vise les appartements loués dont les caractéristiques sont « susceptibles de nuire à la vie ou à la santé des habitants ». En 1891, le docteur Jacques Bertillon, qui publie le premier recensement des conditions de logement et cartographie le surpeuplement des ménages, compare cette carte avec celle des décès (Bertillon, 1894). Les architectes du mouvement moderne, guidé par Le Corbusier, mettront à profit ces constats. Au nom de l'hygiène, ils préconisent des logements aérés, lumineux et respectant des dimensions standardisées devant procurer à l'habitant un sentiment de bien-être et de confort. De son côté, la statistique publique des conditions de logement se fonde depuis l'après-guerre sur la présence d'équipements sanitaires comme indice de salubrité. À partir de 2000, observant que plus de 90 % des logements sont équipés, des chercheurs de l'Insee ont montré que d'autres défauts affectent l'habitat et la santé. Ainsi, la qualité de l'habitat ne se réduit pas à des dimensions techniques. Elle comprend la possibilité pour l'habitant de partir ou de déménager (principe de mobilité) ; le principe d'adaptabilité du logement à tous les âges et cycles de la vie ; le principe d'identité qui permet ou non à l'habitant de se reconnaître dans son lieu de vie et d'en parler (principe de narrativité).

  Ainsi les maladies respiratoires et cardiovasculaires, les dépressions et le stress sont aussi à mettre en rapport avec la qualité de l'air et des matériaux, l'humidité, l'isolation phonique et thermique, le manque d'espace de nos logements. Le rapport de l'Institut des Hautes Études pour l'Action dans le Logement IDHEAL (2021) souligne la méconnaissance des usages et modes de vie de la part des concepteurs. Sur la base d'une analyse de 200 plans d'appartements construits depuis vingt ans, il constate une détérioration des surfaces aménageables dans les appartements, notamment la surface des pièces, les cuisines ouvertes, et les espaces extérieurs.

	Depuis une dizaine d'années, plusieurs équipes d'experts réfléchissent à la mise au point d'un index
	des qualités de l'habitat en matière de santé.
	Cependant, les mesures établies par ces chercheurs, qu'il s'agisse de l'Indice de salubrité du logement
	(2014), Healthy Homes Barometer (2015), HEQI (2019), ou Domiscore (2020), considèrent peu la
	santé comme un état de bien-être physique et mental (définition de l'OMS de 1946), et ne
	s'intéressent pas aux immeubles ordinaires non classés comme insalubres.
	En 2021, le gouvernement a lancé un vaste appel d'idées et d'expérimentations sur « la qualité du

De plus les études nécessitent de tenir compte des évolutions démographiques comme le vieillissement de la population, ou l'augmentation des familles monoparentales qui montrent que les qualités attendues de l'habitat sont aussi plurielles que liées à la santé physique et mentale (chutes, accidents domestiques, stress…).

La qualité architecturale des logements… sans la santé !

Le problème de la qualité des logements est en effet général. Depuis une dizaine d'années, les architectes et les constructeurs réfléchissent à cette question. Après la loi Élan (2018) qui propose de « construire plus, mieux et moins cher », trois rapports officiels abordent le sujet. Le rapport Lemas/Badia (2020), met l'accent sur la question des surfaces des logements sociaux, leur modularité et leur adaptabilité, leur luminosité et leur rapport à l'extérieur. Le rapport Laurent Girometti et François Leclerq (2021) qui soutiennent « un référentiel destiné à améliorer la qualité d'usage des logements pour répondre aux nouveaux besoins des occupants et retrouver le désir d'habiter en ville ». logement de demain » et désigné une centaine de lauréats faisant l'objet d'un suivi scientifique et technique.

Cependant, dans ces consultations, la notion de santé n'apparaît pas ou peu, du moins de façon explicite. Ceci illustre les difficultés des deux domaines du logement et de la santé à dialoguer, ce qui n'a pas toujours été cas.

  Tel est le projet du programme SAPHIR, soutenu par L'Agence Régionale pour la Santé. Il consiste en une rechercheaction dans le cadre d'un partenariat entre le bailleur social Paris Habitat et IDHEAL. Docteur Jacques Bertillon, qui a parmi les premiers au XIXᵉ siècle conçu l'idée de logement sains Docteur Jacques Bertillon, qui a parmi les premiers au XIXᵉ siècle conçu l'idée de logement sains. Wikimedia Focus group 1 à Aubervilliers. SAPHIR 2022, Author provided (no reuse)

Il s'appuie sur une typologie de 12 immeubles représentatifs, par leur taille, l'époque de construction du Nord-est parisien. Ils reflètent une gamme diversifiée de logements sociaux et très sociaux. Auprès des habitants de ces immeubles, il s'agit de mesurer la capacité des habitants à se saisir des questions d'habitat par le biais de la santé, qu'il s'agisse du chauffage, de l'humidité, de la qualité de l'air, des nuisances sonores, etc.

https://theconversation.com/logement-un-enjeu-


Baptiste Régnier

Judith Legrand

François Rebaudo
email: francois.rebaudo@ird.fr

Modeling temperature-dependent development rate in insects and implications of experimental design

Keywords: development rate, laboratory, nonlinear models, temperature

Characterizing the temperature-dependent development rate requires empirical data acquired by rearing individuals at different temperatures. Many mathematical models can be fitted to empirical data, making model comparison a mandatory step, yet model selection practices widely vary. We present guidelines for model selection using statistical criteria and the assessment of biological relevance of fits, exemplified throughout a Lepidoptera pest dataset. We also used in silico experiments to explore how experimental design and species attributes impact estimation accuracy of biological traits. Our results suggested that the uncertainty in model predictions was mostly determined by the rearing effort and the variance in development times of individuals. We found that a higher number of tested temperatures instead of a higher sample size per temperature may lead to more accurate estimations of model parameters. Our simulations suggested that an inappropriate model choice can lead to biased estimated values of biological traits (defined as attributes of temperature dependent development rate, i.e. optimal temperature for development and critical thresholds), highlighting the need for standardized model selection methods. Therefore, our results have direct implications for future studies on the temperature-dependent development rate of insects.

Introduction

Arthropods are ectothermic organisms, so that ambient temperature is a determining factor for their life-history traits such as survival, fecundity, and development. Consequently, the study of arthropod responses to environmental factors has a long history in quantitative ecology, and particularly for insects, considering their role in agroecosystems (Uvarov, 1931). Insect life-history traits are essential components for modeling and predicting insect occurrence, phenology, and population dynamics [START_REF] Ewing | Modelling the effect of temperature on the seasonal population dynamics of temperate mosquitoes[END_REF][START_REF] Haridas | Effects of temporal variation in temperature and density dependence on insect population dynamics[END_REF]Nielsen et al., 2016). Insect phenological models have received much attention in entomology [START_REF] Chuine | Processbased models of phenology for plants and animals[END_REF]Régnière, 2017, Shi et al., 2017b). Used primarily for pest management (e.g. [START_REF] Crespo-Pérez | Modeling invasive species spread in complex landscapes: the case of potato moth in Ecuador[END_REF]Rebaudo et al., 2017;[START_REF] Roy | Relationship between temperature and developmental rate of Stethorus punctillum (Coleoptera: Coccinellidae) and its prey Tetranychus mcdanieli (Acarina: Tetranychidae)[END_REF], but also in forensic science (e.g. Martín-Vega et al., 2017) and vector borne diseases (e.g. Kreppel et al., 2016), phenological models have been developed to predict insect emergence time and duration, distribution, outbreaks, and voltinism [START_REF] Shi | Comparison of thermal performance equations in describing temperature-dependent developmental rates of insects[END_REF]. Additionally, phenological modeling has been an important tool in recent decades for studying the impact of climate change on the population dynamics of arthropod species (e.g. Mwalusepo et al., 2015;Pollard et al., 2020).

Phenological models depend on two elements: i) temperature time series, which can be used to represent past, present or future conditions, and ii) characterization of the relationship between development rate and temperature for each arthropod life stage. For the latter, the relationship is well documented in the literature as nonlinear and takes the form of a thermal performance curve (fig. 1; [START_REF] Huey | Integrating thermal physiology and ecology of ectotherms: a discussion of approaches[END_REF]Sinclair et al., 2016;Uvarov, 1931). The nonlinear variation of development rate according to temperature is usually described using mathematical models with parameters corresponding to different temperatures of biological meaning. Three parameters are typically measured: the critical thresholds for development (CT min and CT max ) and the optimal temperature for development (T opt ). These three parameters are used to compare species' responses to temperature (e.g. Stejskal et al., 2019) or to predict the potential impact of climate change on insect development (Sinclair et al., 2016).

The experimental data used to quantify the relationship between development rate and temperature are usually obtained by using rearing units for insects under controlled conditions. Insects are reared at different constant temperatures and monitored every day (or hour) to estimate the length of each life stage in each individual (or group of 2 Fig. 1: Theoretical performance curve for temperature-dependent development rate. Models generally assume that development rate is null under a critical temperature (CT min ) from which it increases slowly as temperature increases, reaches a temperature zone where the relationship can be considered as linear [START_REF] Campbell | Temperature requirements of some aphids and their parasites[END_REF] until the development rate reaches a maximum at an optimal temperature (T opt ). It then rapidly decreases and is considered null under a critical temperature (CT max ). individuals), despite known limitations on the use of constant versus fluctuating temperatures [START_REF] Colinet | Insects in fluctuating thermal environments[END_REF]. Such experiments have been carried out on a large variety of insect species with different purposes. For example, the characterization of temperature-dependent development rates of Copitarsia incommoda (Lepidoptera: Noctuidae) with nonlinear models has been used to build a risk map of C. incommoda outbreaks in quinoa crops in the Bolivian Andes (Rebaudo et al., 2017). Another example can be found for Phratora vulgatissima (Coleoptera: Chrysomelidae), a herbivorous coleoptera feeding on willow, cultivated in Europe and North America as a perennial biomass crop (Pollard et al., 2020). A phenological model based on the characterization of the relationship between development rates and temperature for the different life stages has been built and its analysis showed that the currently univoltine species could be bivoltine under the RCP8.5 climate change scenario (IPCC, 2014). Since Réaumur's first experiments on caterpillars in the 1730s, many efforts have been made to build mathematical models describing the temperature-dependent development rate (Rebaudo and Rabhi, 2018;[START_REF] Roy | Relationship between temperature and developmental rate of Stethorus punctillum (Coleoptera: Coccinellidae) and its prey Tetranychus mcdanieli (Acarina: Tetranychidae)[END_REF]. One approach toward characterizing the relationship between temperature and development is grounded on the degree-day theory using a simple linear regression model which has been successfully used to predict development rate curve [START_REF] Campbell | Temperature requirements of some aphids and their parasites[END_REF]. However, this characterization is limited to the range of temperatures where development rate can be considered linear and could lead to misinterpretations for insects experiencing temperatures outside this range, which are frequently encountered by insects in regions marked with seasonality or high diurnal thermal amplitude (Pincebourde and Casas, 2015). In addition, because global warming is characterized by an increasing frequency of extreme events (IPCC, 2014), it is essential to evaluate the response of insects to extreme temperatures (Ungerer et al. 1999;[START_REF] Shi | Influence of temperature on the northern distribution limits of Scirpophaga incertulas Walker (Lepidoptera: Pyralidae) in China[END_REF]. Thus, the characterization of the nonlinear parts of the development rate curve for insects has become an important challenge. For that purpose, several nonlinear equations have been proposed to characterize development rate according to temperature (e.g. [START_REF] Briere | A novel rate model of temperature-dependent development for arthropods[END_REF]Lactin et al., 1995;[START_REF] Sharpe | Reaction kinetics of poikilotherm development[END_REF], and more than 30 equations can be found in the literature to date. Comparing models is consequently an essential step in characterizing development rate in insects, but remains to this day an irregular practice (Quinn, 2017).

As described by Lamb et al. (1984), a modeling study should comply with several criteria: (1) the model should fit the data accurately (which can be estimated using statistical measures of goodness-of-fit); (2) the parameters of the model should have a biological meaning; and (3) the model should be parsimonious, i.e. have as few parameters as possible. However, a conclusion from the numerous discussions on model selection (see [START_REF] Damos | Temperature-driven models for insect development and vital thermal requirements[END_REF]Kontodimas et al., 2004;Quinn, 2017;[START_REF] Régnière | Effects of temperature on development, survival and reproduction of insects: experimental design, data analysis and modeling[END_REF][START_REF] Shi | A comparison of different thermal performance functions describing temperaturedependent development rates[END_REF][START_REF] Shi | Comparison of thermal performance equations in describing temperature-dependent developmental rates of insects: (I) empirical models[END_REF]Shi et al., , 2017a;;Zahiri et al., 2010) is that model comparison is mandatory, since it is not possible to find a consensus on one model. Goodness-of-fit assessment plays a central part in model selection, while model selection methods are heterogenous (e.g. using the Aikaike Information Criterion in Rebaudo et al., 2017, Residual Sum of Squares and AIC in Pollard et al., 2020, RSS, R 2 , and AIC in Noor-ul-Ane et al., 2018) and model comparison is not systematic in the literature (Quinn, 2017). The main motivation for this article is to present a guideline to assess and compare model adjustments from empirical datasets. The guidelines are exemplified with an empirical dataset obtained for Tecia solanivora (Lepidoptera: Gelechiidae). Also, because several factors other than the criteria used for goodness-of-fit assessment can play a role in model selection, we explore how experimental design influences parameter estimation accuracy. We studied the experimental design characteristics (number of experimental temperatures, number of individuals) and insect population characteristics (variance among individuals in development time) to identify and quantify the main determinants that optimize temperature-dependent development rate characterization.

The purpose of this work is i) to specify the experimental designs that favor an accurate characterization of the nonlinear relationship between development rate and temperature, ii) to propose a guideline for model assessment. Although these two parts are treated independently, they are complementary by proposing methods on two technical components of the characterization of development rate as a function of temperature. We first assess and compare model fits using simulated datasets corresponding to different scenarios of experimental design, under different biological assumptions. Then we describe methods for model goodness-of-fit assessment and comparison, together with the tools to implement them, using a dataset extracted from the literature. Through these two aspects of temperature-dependent development rate modeling, our goal is to specify the best practices for the characterization of development rate as a function of temperature.

Methods

The development rate model for each life stage can be written as a function of temperature r(T). This function can have several parameters that can be estimated by fitting the model to observed data. A typical experiment consists in rearing n insects at J different temperatures and observing the length of the stage for each individual. We denote as n j the number of insects reared at temperature T j (j=1, … J) and D ij the length of the stage for the ith individual reared at temperature T j . For each individual, the development rate is computed as the inverse of the length of the stage r ij =1/D ij . Below, we denote as r ij the observed values of the development rate of individual i at temperature T j and as r j the prediction from the model r(T j ).

Impact of experimental design, a simulation study

Estimating the model parameters and biological traits from the development rate models depends on two kinds of characteristics: those imposed by the species biology (i.e., the probability distribution of development times at each temperature), and those set by the experimental design (i.e., the number of individuals reared at each temperature, the number, range and distribution of experimental temperatures). The purpose of this study is to define the characteristics of experimental designs that allow an accurate (i.e. precise and unbiased) characterization of the nonlinear relationship between development rate and temperature for insects. We defined an in silico biological scenario with (1) a theoretical model used to describe temperature-dependent development rate, (2) a set of parameter values for this model, and (3) a probability distribution of the development times at each temperature, including the variance parameter value.

Using simulations of different in silico biological and experimental scenarios, we explored how the variance of the development times at each temperature (species attribute), the number of individuals at each temperature, the number of tested temperatures, and the distribution of tested temperatures (experimental design attributes) impact the precision and bias of model parameters and biological traits estimations.

In silico biological and experimental scenarios

All in silico biological scenarios were simulated with the Brière-2 model [START_REF] Briere | A novel rate model of temperature-dependent development for arthropods[END_REF], chosen for its wide acceptance in the literature (e.g. [START_REF] Aghdam | Temperaturedependent development and temperature thresholds of codling moth (Lepidoptera: Tortricidae) in Iran[END_REF][START_REF] Roy | Relationship between temperature and developmental rate of Stethorus punctillum (Coleoptera: Coccinellidae) and its prey Tetranychus mcdanieli (Acarina: Tetranychidae)[END_REF]Tochen et al., 2014).

r (T )=aT ( T -CT min )( CT max -T ) (1/ b ) (1)
The model presents the advantages of having two parameters corresponding to CT min and CT max and T opt can be computed with:

T opt = (2b CT max + (b+1)CT min )+√4 b 2 CT max 2 + (b+1) 2 CT min 2 -4 b 2 CT min CT max 4 b+2
(2)

To mimic real-world values for the theoretical development rate curve, we extracted the set of parameter values of the Brière-2 model fitted on the Cydia pomonella pupal stage from [START_REF] Aghdam | Temperaturedependent development and temperature thresholds of codling moth (Lepidoptera: Tortricidae) in Iran[END_REF] (CT min = 8.2, CT max = 34.4, a = 0.00009, b = 3.9, T opt computed = 31.03°C). For all simulations, we considered that the duration of an individual's life stage for a given temperature follows a log-normal distribution to represent the unimodal asymmetrical distribution typical of development times distributions in insects [START_REF] Régnière | Effects of temperature on development, survival and reproduction of insects: experimental design, data analysis and modeling[END_REF][START_REF] Régnière | A method of describing and using variability in development rate for the simulation of insect phenology[END_REF]Stinner et al., 1975). A scenario includes the rearing at J temperatures (T 1 , … T J ) of n j individuals. At each temperature T j , the mean value of development time (d j = 1/r j ) was computed from the theoretical curve given by equation 1, and the development times of individuals (D ij ) were drawn in a log-normal distribution of parameters μ j and σ² where μ j = log(d j ) -(σ²/2):

log ( D ij ) ∼ N ( log (d j )-(σ 2 /2) , σ 2 ) (3)
The variance of the development time is given by σ Colaphellus bowringi (Coleoptera: Chrysomelidae) (Shi et al. 2017c). We simulated in silico biological scenarios with three different values of the coefficients of variation (10%, 30%, 50%) and chose σ² accordingly. Because experimental measurements of development time are usually made on a daily basis (e.g. Tochen et al., 2014), the simulated development times (in days) were rounded up to the next whole number. To assess the impact of the experimental rearing effort on the accuracy of parameter estimations, we explored the effect of the number of individuals in a range of plausible values. The range was defined based on observed reference values in the literature, with 10, 50, 100 or 500 individuals per temperature [START_REF] Aghdam | Temperaturedependent development and temperature thresholds of codling moth (Lepidoptera: Tortricidae) in Iran[END_REF][START_REF] Briere | Comparison of temperaturedependent growth models with the development of Lobesia botrana (Lepidoptera: Tortricidae)[END_REF][START_REF] Howell | Physiological development time and zero development temperature of the codling moth (Lepidoptera: Tortricidae)[END_REF]Pollard et al., 2020;[START_REF] Roy | Relationship between temperature and developmental rate of Stethorus punctillum (Coleoptera: Coccinellidae) and its prey Tetranychus mcdanieli (Acarina: Tetranychidae)[END_REF]Zahiri et al., 2010). Based on a sample of the references within a review on the relationship between temperature and development rate (Rebaudo and Rabhi, 2018), we chose to explore the impact of the number of experimental temperatures with 5, 7 and 15 different temperatures in the range of the expected values (from CT min to CT max ), with different distributions of experimental temperatures. Temperatures were distributed differently depending on three scenarios. In the first one (referred to as "uniform"), temperatures were equally distributed within the range of the critical thermal thresholds, as exemplified in numerous studies (e.g. Kontodimas et al., 2004;Messenger and Flitters, 1958;Noor-ul-Ane et al., 2018;[START_REF] Roy | Relationship between temperature and developmental rate of Stethorus punctillum (Coleoptera: Coccinellidae) and its prey Tetranychus mcdanieli (Acarina: Tetranychidae)[END_REF]Zahiri et al., 2010). In the second scenario (referred to as "linear"), temperatures were equally distributed within a range of temperatures corresponding to the linear zone. This choice of experimental design is common in the literature because extreme temperatures are rarely observed in the habitat, and because there is a high mortality rate at extreme temperatures often leading to a lack of data (e.g. [START_REF] Campbell | Temperature requirements of some aphids and their parasites[END_REF][START_REF] Howell | Physiological development time and zero development temperature of the codling moth (Lepidoptera: Tortricidae)[END_REF]Zuim et al., 2017). To obtain temperatures in this range, the thermal safety margin (TSM), defined in Sinclair et al. (2016) as the difference between CT max and T opt (TSM = CT max -T opt ), was used. On the basis of the hypothesis that insects avoid temperatures close to critical thermal thresholds (Martin and Huey, 2008), temperatures were distributed between CT min + TSM and CT max -TSM in this scenario. The third scenario (referred to as "specific") corresponds to an alternative to classical experimental designs, with temperatures unequally distributed within the range of the critical thermal thresholds. To generate this scenario, we sampled additional temperatures in a uniform distribution and selected the requested number of temperatures closest to the thermal thresholds (see Supplementary Material 2 lines 95 to 113). As suggested by Kingsolver and Buckley (2020), the nonlinearity of development curves in extreme temperatures could be better described by experimental design including more low and high temperatures than mid-temperatures.

Assessing bias, precision and accuracy of parameter estimators

We combined the resulting temperature experimental designs with the four possible sample sizes at each temperature and the three values for the coefficient of variation in a complete factorial design leading to 108 scenarios. For each scenario, 500 repetitions were simulated, leading to 54000 in silico simulations. For each simulation, the Brière-2 model was fitted to the simulated data. With each model fit, we estimated the values of the four parameters of the model (a, b, CT min , CT max ) and T opt , and compared these estimates to the theoretical values with:

error ( pkz ) = pkz -p k (4)
where pkz is the estimation of the parameter p k (p denoting one of the parameters) with dataset k of scenario z.

For each parameter and each of the 108 scenarios, we quantified the bias of the estimator as the median error (ME), its precision as its interquartile range (IQR) and its accuracy by the mean of absolute errors (MAE). Bias evaluates the systematic error in parameter estimations, precision depends on the variance of the estimator, and accuracy is defined as the overall difference between the estimated parameter values and their true value (Walther and Moore, 2005). Bias and accuracy were measured using median error and mean of absolute errors, respectively, because they are both less sensitive to outlier values than more common measures like mean error and mean of squared errors. The complete code is available in Supplementary Material 2.

Impacts of misspecification

In order to estimate the impacts of misspecification (i.e. fitting another model than Brière-2), we fitted seven models for which two parameters corresponded to CT min and CT max , in addition to the Brière-2 model (Analytis model from [START_REF] Analytis | Über die relation zwischen biologischer entwicklung und temperatur bei phytopathogenen pilzen[END_REF]; Performance model from [START_REF] Shi | A simple model for describing the effect of temperature on insect developmental rate[END_REF]Shi model from Shi et al., 2011; Brière-1 model from [START_REF] Briere | A novel rate model of temperature-dependent development for arthropods[END_REF]Wang model from Wang et al., 1982;Ratkowsky model from Ratkowsky et al., 1983;Kontodimas model from Kontodimas et al., 2004). The models were fitted to 500 simulated datasets with CV = 10% and a uniform distribution of 15 temperatures and 500 individuals per temperature. For each fit, the parameters corresponding to CT min and CT max were extracted, and T opt was approximated from the fitted model. Then, the bias was measured with the median of errors (equation 8).

Model fit, assessment and comparison

Models characterizing the relationship between development rate and temperature in insects differ in complexity, as well as in their biological assumptions [START_REF] Davidson | On the relationship between temperature and rate of development of insects at constant temperatures[END_REF][START_REF] Sharpe | Reaction kinetics of poikilotherm development[END_REF], Scoolfield et al. 1981, Brière 1999[START_REF] Ikemoto | Intrinsic optimum temperature for development of insects and mites[END_REF][START_REF] Régnière | Effects of temperature on development, survival and reproduction of insects: experimental design, data analysis and modeling[END_REF]. As a result of the variety of available models, the process of model selection is a complex task. Here we propose guidelines for (1) model fitting, (2) statistical model selection, and (3) biological relevance of fits.

Methods for model fitting

From the vast literature on insect development rate (see the review by Rebaudo and Rabhi, 2018), the nonlinear least squares method (NLS) emerges as the main method to adjust insect development rate models to empirical data (but see [START_REF] Got | Temperature-dependent model for European corn borer (Lepidoptera: Pyralidae) development[END_REF][START_REF] Hansen | Prepupal diapause and instar IV developmental rates of the spruce beetle, Dendroctonus rufipennis (Coleoptera: Curculionidae, Scolytinae)[END_REF][START_REF] Régnière | Effects of temperature on development, survival and reproduction of insects: experimental design, data analysis and modeling[END_REF]Yurk and Powell, 2010 for examples of adjustments based on maximum likelihood estimations and Shi et al. 2017c for NLS limitations). NLS consists in iterative optimization using algorithms to identify equation parameter values minimizing the sum of squares of errors. These are traditional methods not specific to development rate curves for insects, and several algorithms exist for this purpose (Nocedal and Wright, 2006), the most widespread being the Gauss-Newton algorithm [START_REF] Björck | Numerical methods for least squares problems[END_REF] and the Levenberg-Marquardt algorithm (Levenberg, 1944;Marquardt, 1963). In R statistical software (R Core Team, 2021), the Gauss-Newton and Levenberg-Marquardt algorithms are respectively implemented through the packages stats (R Core Team, 2021) and minpack.lm [START_REF] Elzhov | minpack.lm: R interface to the Levenberg-Marquardt nonlinear least-squares algorithm found in minpack, plus support for bounds[END_REF].

Statistical model selection

The first step in fit assessment is to analyze the precision of parameter estimates. For that purpose, standard errors of estimates can be computed. In a second step, statistical selection criteria can be used to select one or several model candidates among all fitted models. These criteria are grounded on goodness-of-fit and the parsimony principle (i.e. models must fit the data accurately with as few parameters as possible). Goodness-of-fit quantifies how well a model fits to a set of experimental data using statistical criteria. Here, we present two measures of the discrepancy between the predictions of the model and the observed development rates, and two information criteria that can be used for model selection.

The Residual Sum of Squares (RSS) is the sum of the squares of the differences between values of the development rate predicted by the model ( r j ) and the observed values (r ij ). It has been vastly used for insect development rate curve characterization over recent decades (e.g. Kontodimas et al., 2004;Pollard et al., 2020;[START_REF] Shi | A comparison of different thermal performance functions describing temperaturedependent development rates[END_REF]. RSS value is scale-dependent, meaning that a threshold discriminating goodness-of-fit cannot be defined. However, RSS can be used as a comparison tool in order to identify models with a higher difference between predicted and empirical data and is calculated with:

RSS= ∑ j=1 J ∑ i=1 n j ( r ij -r j ) 2 (5)
The root-mean squared error (RMSE; equation 2) is defined as:

RMSE= √ RSS n (6)
Its usage is similar to RSS for characterizing insect development rate model accuracy (e.g. [START_REF] Damos | Temperature-dependent bionomics and modeling of Anarsia lineatella (Lepidoptera: Gelechiidae) in the laboratory[END_REF]Mirhosseini et al., 2017;[START_REF] Shi | Comparison of thermal performance equations in describing temperature-dependent developmental rates of insects: (I) empirical models[END_REF]. In contrast with RSS, RMSE has the advantage of being expressed in the same unit as the outcome, making it easier to interpret. It can also be normalized using the mean or median development rate. Both RMSE and RSS can be used to rank model fits and to select the model that minimizes RMSE and RSS values.

To compare models fitted on the same dataset integrating a trade-off between model complexity (i.e. the number of parameters K) and goodness-of-fit, the Akaike Information Criterion (AIC; [START_REF] Akaike | A new look at the statistical model identification[END_REF] and Bayesian Information Criterion (BIC; [START_REF] Schwarz | Estimating the dimension of a model[END_REF]) can be used. Because goodness-of-fit tends to be maximized by adding more parameters to a model, both of these criteria consider a penalty parameter dependent on the number of parameters. In the case of least squares estimations of model parameters [START_REF] Burnham | Model selection and multimodel inference: a practical informationtheoretic approach[END_REF], the AIC and BIC can be calculated with:

AIC=n × log ( RSS n ) + 2 K (7) BIC=n ×log ( RSS n ) + log (n) × K (8)
The BIC involves a larger penalty than the AIC and penalizes models with a high number of parameters more than the AIC. The AIC is a common tool for selecting a model for insect development rate curve characterization (e.g. Orang et al., 2014;Quinn, 2017;Tran et al., 2012). The BIC is less popular but has been used several times in the same way (e.g. [START_REF] Arbab | Determining suitability of thermal development models to estimate temperature parameters for embryonic development of Sitona lepidus Gyll. (Coleoptera: Curculionidae)[END_REF].

According to [START_REF] Burnham | Model selection and multimodel inference: a practical informationtheoretic approach[END_REF] the goodness-of-fit of models resulting in close AIC values can be considered equivalent. The study proposed a rule of thumb to discriminate between models using the differences between the lowest AIC value and each model's AIC value (noted Δ AIC ): the goodness-of-fit of models resulting in Δ AIC < 2 can be considered equivalent to the model with the lowest AIC, while models resulting in Δ AIC > 10 can be discarded as their goodness-of-fit is low in comparison with the model with the lowest AIC. In the case of Tecia solanivora, we chose to select models based on their Δ AIC values, and models resulting in Δ AIC > 10 were discarded. The advantage of this method is that several models can be selected and the resulting biological trait estimations can be compared.

Biological relevance of fits

The AIC and BIC allow us to choose the model that represents the best trade-off between model complexity and goodness-of-fit among a set of models, but not necessarily a model that would be a good representation of the biological process. A good description of data according to statistical measures can result in biologically unrealistic development rate curves with aberrant values of CT min and CT max , notably when data are sparse in low and high temperatures [START_REF] Ikemoto | Novel method of specifying low and high threshold temperatures using thermodynamic SSI model of insect development[END_REF]. Several examples of evaluation of biological traits values can be found in the literature, but the practice does not seem to be widespread (Orang et al. 2014;Kontodimas et al. 2004). Complementary knowledge about the biology of the studied species could help to clarify such criteria. For instance, in the case of Tecia solanivora, we chose to discard models that resulted in CT min inferior to 5°C, based on measured survival rates in low temperatures [START_REF] Crespo-Pérez | Modeling invasive species spread in complex landscapes: the case of potato moth in Ecuador[END_REF]. Because individuals cannot survive in temperatures lower than 5°C, we made the assumption that estimations of minimal critical threshold should be higher than 5°C. Thus, a limitation on the use of biological criteria for model selection is that it requires preliminary results on the temperaturedependent development rate and on the species' response to temperature in order to discard aberrant model adjustments. Here we used published results on survival rates in low temperatures, as seen in Nephus includens and Nephus bisignatus (Kontodimas et al. 2004). However, the most relevant biological criteria could depend on the species and on available knowledge. Moreover, filtering fits based on biological criteria requires strong assumptions and can influence the interpretation of results. In our case we cannot conclude that development is not possible in temperatures under 5°C based on development rate curve interpretation, since we made this assumption for model selection. Thus, biological criteria should always be justified with regard to the knowledge of the species, and vigilance is required in interpreting results. To add biological information in the modeling process, other methods could be used. For instance, it is possible to add a constraint on a parameter estimation in an NLS procedure. The algorithm look for the parameter estimations that minimizes the residual sum of squares in respect of the constraint (e.g. CT min ≤ 5°C). Alternatively, parameters could be estimated in a bayesian framework with informative priors for parameters such as CT min . Such approach was used to estimate upper and lower thresholds in the mountain pine beetle Dendroctonus ponderosae (McManis et al., 2018). Both methods however prevent the use of models that do not include parameters corresponding to biological traits, and thus the generalization of our approach. In the example presented here, we used 14 models that do not directly estimate CT min and/or CT max and/or T opt (e.g. the polynomial model of degree four has 5 parameters but none correspond to a biological trait). For these models, we computed CT min and CT max by taking respectively the minimal and maximal temperature where the predicted development rate was null, and we computed T opt using the temperature where the predicted development rate took its maximal value (local maximal value). Yet, we acknowledge that for some models, these biological traits could be written analytically as a function of the model parameters (see Brière-2 model and equation 2 which was given in Brière et al., 1999).

Another limitation for assessing the biological relevance of fits is that some models result in an asymptotic curve around zero in low and/or high temperatures (e.g. models in Logan et al., 1976;[START_REF] Sharpe | Reaction kinetics of poikilotherm development[END_REF]Taylor, 1981). In order to overcome this limitation, we propose to use alternative metrics corresponding to threshold temperatures where the development rate is below a given fraction of the maximal development rate. The biological meaning of these metrics is different from the temperature thresholds where development rate is null, but they can be used similarly to CT min and CT max for discarding biologically unrealistic curves, and can be computed with every model. In the case of T. solanivora, we measured threshold temperatures where the development rate is inferior to 10% of the maximal development rate, noted CT10 min and CT10 max . Biological trait values can be known either with the value estimated through a parameter of a model, or directly from the fitted curve. When available, standard errors of biological traits estimations were computed, but in some cases no standard errors could be computed because their value were computed using the model predictions, as it is the case for CT10 min and CT10 max .

To illustrate the whole process of model fitting, statistical selection and biological assessment, we estimated and selected models for the larval stage of Tecia solanivora (Lepidoptera: Gelechiidae) using data published in [START_REF] Crespo-Pérez | Modeling invasive species spread in complex landscapes: the case of potato moth in Ecuador[END_REF]. The dataset is composed of 14 observed mean development rates measured at 9 different temperatures (see [START_REF] Crespo-Pérez | Modeling invasive species spread in complex landscapes: the case of potato moth in Ecuador[END_REF]. Tecia solanivora is a potato tuber moth present in Central and South america. No diapause has been reported in the Andean region near the equator probably due to the small thermal differences observed between the seasons at these latitudes (Dangles et al. 2008).

For all the analyses, we used the R software version 4.0.3 and the devRate package version 0.2.1 [START_REF] Rebaudo | Modelling temperaturedependent development rate and phenology in arthropods: The devRate package for R[END_REF]. The complete code is available in Supplementary Material 1. The dataset for T. solanivora is included in the package devRate.

Results

Experimental design impacts on estimations of biological traits 3.1.1. Bias and precision of estimations

In silico insect rearing experiments simulated with a theoretical curve of the Brière-2 model and fitting of the same model on simulated datasets showed that accuracy of parameter estimations depended on the coefficient of variation of the distribution of development times (fig. 3b), and on the characteristics of experimental design (fig. 3cde). Overall, biological trait estimations were unbiased, as the medians of errors were close to 0 for the three biological traits and the three coefficients of variation of development times (fig. 3b). However, the precision of estimations depended on the biological traits and scenario. CT max and T opt had higher interquartile ranges of errors than CT min (from 0.1204°C to 1137.27°C for CT max ; from 0.0724°C to 4.44°C for CT min ; from 0.0421°C to 1132.67°C for T opt depending on the scenario). With the experimental design with linear distribution of 5 temperatures and 10 individuals per temperature, which led to the highest interquartile range of errors for the three biological traits, the maximum error found for CT min was equal to 5.14°C. With this scenario resulting in a lack of data in high temperatures, T opt and CT max estimations took extreme values resulting in high error values superior to 10 4 °C.

Impacts of experimental design

Accuracy of biological trait estimations depended on the rearing effort. For the three values of coefficient of variation, accuracy increased with increasing number of temperatures, and with increasing number of individuals per temperature. With CV = 10%, CT max mean absolute error (MAE) varied from 23.87°C with 5 temperatures and 10 individuals, to 0.15°C with 15 temperatures and 500 individuals. With CV = 30%, CT max MAE varied from 851.74°C with 5 temperatures and 10 individuals, to 0.26°C with 15 temperatures and 500 individuals. With CV = 50%, CT max MAE varied from 1206.09°C with 5 temperatures and 10 individuals, to 0.42°C with 15 temperatures and 500 individuals. Similar results were observed for T opt estimations, but the values of MAE were lower than for CT max estimations (fig. 3de). For CT min , similar results were observed, however the values of MAE were generally low (< 1°C), except with poor experimental designs (e.g. 5 temperatures and 10 individuals per temperature resulted in the highest value of MAE for CT min , equal to 2.9°C). The distribution of temperatures also affected the accuracy of biological trait estimations. Linear distribution of temperatures resulted in higher MAE values than uniform and specific distributions with every scenario (fig. 3cde). Specific distribution of temperatures resulted in lower or equal values of MAE compared to uniform distribution (fig. 3cde).

With CV = 10%, a total of 150 individuals (with 15 temperatures and 10 individuals per temperature and specific or uniform distributions of temperatures) resulted in CT max MAE < 1°C (fig. 3d). With a specific distribution of 7 temperatures, at least 100 individuals per temperature (for a total rearing of 700 individuals) were required to obtain a CT max MAE < 1°C (fig. 3d). With the same number of temperatures and a uniform distribution, at least 500 individuals (for a total of 3500 individuals) were required to obtain CT max MAE < 1°C (fig. 3d). The minimal value of CT max MAE for the linear distribution (MAE = 1.2°C) was found with 15 temperatures and 500 individuals per temperature (fig. 3d). With CV = 30%, a total of 750 individuals (with 15 temperatures and 50 individuals per temperature and specific distribution of temperatures) resulted in CT max MAE < 1°C. With a uniform distribution, at least the double of individuals were required to obtain CT max MAE < 1°C (with 15 temperatures and 100 individuals per temperature). With a specific distribution of 7 temperatures and 500 individuals per temperature (for a total of 3500 individuals), CT max MAE was equal to 1.15°C. With a linear distribution of temperatures, the minimum value of MAE (2.9°C) was found with 15 temperatures and 500 individuals per temperature.

With CV = 50%, the minimal required experimental design to result in MAE < 1°C was a specific distribution of 15 temperatures and 100 individuals per temperature (total of 1500 individuals). With a uniform distribution of temperatures, 15 temperatures and 500 individuals per temperature (total of 7500 individuals) resulted in MAE < 1°C. With a linear distribution, the minimum value of CT max MAE found (MAE = 5.79°C) was for 15 temperatures and 500 individuals per temperature.

Table 1: Models fitted to development rate data of larval stage of T. solanivora. The reference for each model is shown, as well as all the criteria for model selection. RSS and RMSE are measures of the discrepancy between observations and predictions. AIC and BIC are used to select the most parsimonious models. CTmin, CTmax, Topt, CT10min and CT10max are used to discard models with unrealistic biological interpretations. Models are ordered according to their AIC values. Column K corresponds to the number of parameters. Standard error of parameter estimations are given in parenthesis when biological traits corresponded to a model's parameter. Other biological trait values were computed from the curves. Models marked with *** are those selected on the basis of statistical and biological criteria.

Model name

Impacts of model misspecification

We fitted eight models to datasets simulated with CV = 10%, a uniform distribution of 15 temperatures and 500 individuals per temperature in order to measure the bias of estimated values of CT min , CT max and T opt when another model than Brière-2 is fitted. The simulations resulted in different values of median error for CT min , T opt and CT max estimations (fig. 4). For CT min , negative median errors were found for four models, and varied from -3.1°C with the Rakowsky model and -0.07 with the Brière-2 model. Positive median errors varied between 0.5°C with the Shi and Wang models, and 1.8°C with the Performance model. For CT max , the Analytis model resulted in a negative median error equal to -0.17°C and for the seven other models, positive median errors were found between 0.28°C with the Brière-2 and 10.4°C with the Kontodimas model. For T opt , four models resulted in negative median errors from -0.035°C with the Brière-2 model and -0.18°C with the Shi and Wang models. The four other models resulted in positive median errors, between 0.028°C with the Analytis model and 1.2°C with the Kontodimas model. 14 Fig. 2: Relationship between development rate and temperature of Tecia solanivora larval stage, with eight models. Model selection with statistical and biological criteria resulted in eight models (see table 1 for model names and references). For each model, the fitted curve is represented. The grey area corresponds to the area between CT10 min and CT10 max , corresponding to critical thresholds where development rate is inferior to 10% of the maximal development rate. The vertical black line corresponds to the temperature at which development rate takes its maximal value (T opt ).

Model selection for temperature-dependent development rate of Tecia solanivora larval stage

The dataset of Tecia solanivora larval stages was used to estimate parameters of 26 models using the Levenberg-Marquardt algorithm (table 1). Standard errors of parameter estimations were important for most models (e.g. standard errors of CT min estimations varied between 1.86°C and 1735.4°C; for CT max , they varied from 0.01°C to 28.6°C; see table 10 in Supplementary Material 1 for a complete overview of standard errors of model parameters) because of the small sample size (N = 14), suggesting that parameters estimations may be biased. The model fits were sorted based on statistical measures of goodness-of-fit. Among the 26 models, the Janisch model minimized both the AIC and BIC [START_REF] Janisch | Relationship between temperature and development rate of Copitarsia incommoda (Lepidoptera: Noctuidae) in the Bolivian Andes[END_REF]. Fifteen model fits resulted in AIC values with a difference from the lowest value higher than ten and were therefore excluded. The eleven resulting model fits revealing no anomalies in their RSS and RMSE values (table 1) were then filtered based on biological relevance of biological traits. Because CT min and CT max could not be computed for some models, we analyzed the values of CT10 min and CT10 max . Based on observed development rates and survival at low temperatures [START_REF] Crespo-Pérez | Modeling invasive species spread in complex landscapes: the case of potato moth in Ecuador[END_REF], the minimal temperature threshold where development rate is lower than 10% of the maximal development rate (CT10 min ) should be higher than 5°C and lower than 10°C, as it is the minimal temperature at which positive values of development rate were measured. Two models demonstrated a value lower than 5°C (the Janisch model and the Damos model, Damos and Savopoulou-Soultani, 2012; table 1) and were discarded. The polynomial model of degree four was also discarded as no CT min nor CT10 min values were predictable because the model fit predicted an increasing development rate in temperatures inferior to 10°C. Based on the statistical and biological filtering, seven models [START_REF] Analytis | Über die relation zwischen biologischer entwicklung und temperatur bei phytopathogenen pilzen[END_REF]Kontodimas et al., 2004;Ratkowsky et al., 1983;Regniere et al., 2012;[START_REF] Schoolfield | Non-linear regression of biological temperaturedependent rate models based on absolute reaction-rate theory[END_REF]Taylor, 1981;Wang et al., 1982) and the polynomial model of degree two were proposed as potential candidates to represent T. Solanivora TPC. Those model fits differed with regard of their biological traits (fig. 2). Depending on the model, the predicted optimal developmental temperature (T opt ) for the larval stages of T. solanivora varied between 21.4°C, for the polynomial model of degree two and the Taylor model, and 27°C with the Schoolfield model. Except for the Schoolfield model, all the models estimated T opt values between 21.4°C and 24.4°C. CT10 min varied between 6.1°C with the Kontodimas model and 9.9°C with the Régnière model. CT10 max varied between 28.1°C with the Wang model and 33.6°C with the Kontodimas model. The MAE measures the accuracy of estimations, the higher the value, the lower the accuracy. The dotted line corresponds to a total number of individuals equal to 500, so that points above the line correspond to a total number of individuals greater than 500, and points below the line correspond to a total number of individuals lower than 500. d) Same results for CT max estimations, e) same results for T opt estimations. Fig. 4: Bias for CT min , T opt and CT max estimated with eight models, including seven misspecified models. Bias was measured with median errors for three biological traits (CT min , T opt and CT max ) with eight models fitted to 500 datasets simulated with CV = 10%, 15 temperatures distributed uniformly inside the range of theoretical critical thresholds, and 500 individuals per temperature. See table 1 for model names and references.

Discussion

In this article, we first investigate how species characteristics and experimental design impact inference power and accuracy. We then illustrate the whole process of development rate model selection using a dataset on Tecia solanivora. Model selection using statistical and biological criteria on Tecia solanivora larval stages led to the selection of eight models, resulting in different curves and values of biological traits (fig. 2,table 1). This example highlights the importance of model comparison, as eight candidate models were selected among the 26 models that were originally fitted. In the absence of additional data, the selection could not be refined, yet the assessment of the biological relevance of fits allowed us to discard two models that were biologically unrealistic based on knowledge of the survival rate at low temperatures [START_REF] Crespo-Pérez | Modeling invasive species spread in complex landscapes: the case of potato moth in Ecuador[END_REF]. This result exemplifies the benefits of using biological criteria for model selection in insect temperature-dependent development rate characterization. Because of the biological criteria we chose (CT min > 5°C), the best fitting model according to statistical criteria (Janisch model) measures of goodness-of-fit may not be sufficient for model selection, as they could lead to the selection of model fits that are biologically unrealistic. Despite this filtering of model fits with statistical and biological criteria, several models were validated and discrepancies in the predicted values of biological traits were observed (table 1). Such discrepancies could have consequences on the predictions of phenology models, depending on the chosen temperature-dependent development rate model fit. On the basis of available experimental data, this highlights the uncertainty concerning T. solanivora responses in low and high temperatures, so that vigilance is required when interpreting the predictions of phenology models with low and high temperatures. Despite high variations of T opt estimated values (from 21.4°C to 27°C), development rate curves were similar in mid-temperatures (fig. 2), where development rate data are less sparse. This result raises the question of whether the choice of model is important when temperatures are mild. Yet, additional data would be required to simulate the development of individuals in low and high temperatures, and a complementary rearing experiment of T. solanivora at different constant temperatures could be required. This raises the question whether data obtained for simulating insects phenology, as is the case for T. solanirova [START_REF] Crespo-Pérez | Modeling invasive species spread in complex landscapes: the case of potato moth in Ecuador[END_REF], are adequate for predicting the impacts of global warming, as estimating the development rate curve in high temperatures is uncertain due to data sparsity. In addition, the standard errors of parameter estimations were generally high, suggesting that fits were not precise. Thus, the observed differences in the estimated values of biological traits may not be directly interpreted. This result highlights an uncertainty in the T. solanivora temperaturedependent development rate. The cause of the high standard errors may be the small sample size (N = 14), due to the use of mean development rates instead of raw individual dataset. A dataset including individual values of the development rate may be needed for further investigations on the T. solanivora development rate in function of temperature, in order to estimate model parameters more precisely. For the generalization of the procedure, here we used the AIC knowing that it can favor complex models with small dataset. For specific cases with small datasets, the AIC corrected (AICc) could benefit to model selection [START_REF] Burnham | Model selection and multimodel inference: a practical informationtheoretic approach[END_REF].

The model selection process presented here allows the characterization of development rate as a function of temperature. As an additional step toward developing phenological models, diapause induction, duration and termination should be taken into account when required that often depends on day length (Pollard et al., 2020) and temperature [START_REF] Gray | Diapause in the gypsy moth: a model of inhibition and development[END_REF]. As exemplified by our results on T. solanivora, it is not always possible to determine which model is the most suitable, and several model adjustments may be used for development rate characterization. The use of multiple models to forecast phenology, as seen for the prediction of species distributions [START_REF] Araujo | Ensemble forecasting of species distributions[END_REF], could be an interesting development for insect development rate curve characterization and phenology modeling. The guidelines proposed here summarize established methods for model selection in the case of temperature-dependent development rate. None of these methods can assess the predictive power of models, which could be useful for the evaluation of models of the temperature-dependent development rate, as they are mostly used for predicting development and phenology in function of environmental temperatures. For example, cross-validation methods [START_REF] Arlot | A survey of cross-validation procedures for model selection[END_REF] could be an approach for assessing the predictive power of models of the temperature-dependent development rate. Also, developing tools to obtain confidence intervals for biological traits values would improve their interpretation.

The impacts of experimental design on biological trait estimations were explored using simulations of different biological and experimental scenarios. The results showed that accuracy of biological trait estimations depended on the variance in individual development times, as well as the characteristics of the experimental design (number and range of temperatures, and number of individuals per temperature; fig. 3). Accuracy of biological traits estimations can be maximized by maximizing the rearing effort i.e. by increasing the number of temperatures, the number of individuals per temperature and the range of temperatures (fig. 3). Our results suggest that the minimal required rearing effort to obtain low MAE values depended on the variation in development times among individuals. In addition, our results suggest that for the same total number of individuals, a high number of temperatures instead of a high number of individuals per temperature result in lower values of MAE. With a low value of coefficient of variation of development times (CV = 10%), the required rearing effort to obtain a MAE lower than 1°C with uniform or specific distributions of temperatures was of only 150 individuals (with 15 temperatures and 10 individuals per temperature). With 7 temperatures, the required number of individuals per temperature to obtain equivalent MAE values was ten times higher for a specific distribution and 50 times higher for a uniform distribution. This result can be explained by the fact that a high number of temperatures might allow us to obtain more data in the nonlinear parts of the curve, resulting in a more accurate characterization of the temperature-dependent development rate. However, we did not simulate the mortality of individuals, which is generally higher in extreme temperatures than in mid-temperatures (e.g. Rebaudo et al., 2017;[START_REF] Régnière | Effects of temperature on development, survival and reproduction of insects: experimental design, data analysis and modeling[END_REF]Zahiri et al., 2010). The required number of individuals per temperature might be higher in experimental conditions, as mortality in extreme temperatures can vastly reduce the number of individuals per temperature (e.g. in Rebaudo et al., 2017, 49 individuals at pupal stage survived at 30°C against 98 at 18.1°C).

Mean absolute errors of CT min estimations were always far lower than for CT max and T opt estimations, which suggests that the required experimental design depends on the pursued objectives. If the estimation of CT min is required, the rearing effort could be reduced: for instance, with 5 temperatures distributed uniformly and 50 individuals per temperature, the mean absolute error found for CT min was equal to 0.31°C (with a coefficient of variation in development times equal to 10%). Our results suggest that temperatures in the nonlinear zone close to CT min are not required in order to estimate accurately its value with the nonlinear Brière-2 model [START_REF] Briere | A novel rate model of temperature-dependent development for arthropods[END_REF], as the distribution of temperatures in the linear zone resulted in low mean absolute errors of CT min (from 1.22°C to 0.1°C with CV = 10%). However, if CT max and T opt estimations are needed, the rearing effort can be high when the coefficient of variation in development times is superior to 30% (at least 750 individuals with 15 temperatures and 50 individuals per temperature). In addition, data are required in the zone where development rates decrease (between T opt and CT max ). In the absence of data in this zone, the errors for the estimations of CT max and T opt could take values higher than 10 4 °C, which suggests that the fitting of nonlinear models in this case is not recommended. The cost of rearing individuals is directly related to the time necessary for insects to develop. Because at low temperatures the development times is generally high (e.g. the mean development time at 10°C was of 272 days in our simulations), a high effort may be necessary to obtain data in extreme temperatures.

In the absence of mortality, a distribution of temperatures with an effort on low and high temperatures, as suggested in Kingsolver and Buckley (2020), resulted in lower values of mean absolute error for CT max estimation (fig. 3). It is a consistent result, given that this type of temperature distribution allows us to obtain more data in the nonlinear zones of the curve, which could improve the fitting of models in extreme temperatures. However, because mortality is usually higher in extreme temperatures than in midtemperatures, this type of temperature distribution might be associated with a high number of individuals per temperature, in order to compensate for the higher mortality rate. Another limitation of our simulation study is that only one development stage is simulated. In laboratory experiments, individuals are raised at constant temperatures throughout their development. Mortality during development often leads to a higher number of individuals per temperature in the first stages of development than in the last stages (e.g., in Zahiri et al., 2010, the experiment started with 100 eggs, and at pupal stage there were between 13 and 82 individuals that survived depending on the temperature). Our simulations suggest that the accuracy of the estimates of biological traits depended on the number of individuals per temperature, so it is questionable whether the fitted development rate curves of the late immature stages are generally less accurate than those of the early stages.

For characterizing the relationship between development rate and temperature, our results suggest that the required rearing effort could be estimated with preliminary experiments to identify the variance in development times. Yet, general knowledge of the inter-individual variations in development times among insects is still needed. For instance, we made the hypothesis that the coefficient of variation of development times did not vary with temperature, yet this hypothesis still needs to be verified. A constant coefficient of variation may not be a valid assumption, especially for species with a diapause induced by temperature, as a certain temperature threshold may induce diapause in all individuals (e.g. Wall, 1974;Kimberling and Miller, 1988) and lead to a limited variance in development times, which is incompatible with a constant coefficient of variation that implies that the variance of development times increases with mean development time. The NLS method assumes that the residuals of a fit are distributed normally and that their variance is homogeneous, which differ from the design we chose for the simulations where we assumed that the duration of each stage has a constant insects (Kingsolver and Buckley, 2020). Therefore, our results on the impact of the experimental design on the estimation of CT max may contribute to the study of insect response to climate change. In addition, our results suggest that inappropriate model choice can lead to biased values of biological traits, and thus highlight the need for shared and standardized model selection methods.

  of variation by cv j =σ d j /d j =√( e σ

Fig. 3 :

 3 Fig. 3: Impacts of species attributes (coefficient of variation of development times distribution) and experimental design (number and range of temperatures, number of individuals per temperature) on estimations of biological traits. a) Summary of the different aspects of a dataset explored in this simulation study. b) Distributions of errors for each biological trait, in function of the CV of development times. Outlier values are not represented. c) Mean absolute errors (MAE) for CT min estimations with a CV of development times of 10%, in function of experiment attributes (number of temperatures, number of individuals per temperature and distribution of temperatures).The MAE measures the accuracy of estimations, the higher the value, the lower the accuracy. The dotted line corresponds to a total number of individuals equal to 500, so that points above the line correspond to a total number of individuals greater than 500, and points below the line correspond to a total number of individuals lower than 500. d) Same results for CT max estimations, e) same results for T opt estimations.

  was discarded. This result suggests that statistical 17

		10	CT min error T opt error
			CT max error
	Median error	5 0
		-5
		analytis_77	briere2_99	perf2_11	shi_11	briere1_99	wang_82	ratkowsky_83	kontodimas_04

-1). For the sake of simplicity, we considered a constant coefficient of variation for the development time at all temperatures. This constant CV was notably found for the different development stages of the cabbage beetle

0.0067 0.0218 -97.1 -93.9 1.1 (25.4) 35.1 (0.01) 30.4 7.7

Acknowledgments

This work benefited from the support of the project PI2P ANR-19-CE32-0001 of the French National Research Agency (ANR).

coefficient of variation. The constant coefficient of variation for the duration of each stage produces heteroscedasticity as the variance of development rate depends on the temperature and the mean development rate (Shi et al. 2017c). If the NLS method may lead to biased estimations of model parameters, it has been chosen for its broad acceptation in the literature while acknowledging its limitations (Shi et al. 2017c). In addition, we chose a vast range of coefficients of variation (10%, 30%, 50%) in the absence of knowledge about inter-individual variations, and the real range of coefficient of variation remains uncertain. Our result suggests that with a low coefficient of variation (CV = 10%), the required rearing effort for an accurate characterization of biological traits might be lower than what is generally achieved (e.g. a total of 1300 individuals were reared in Zahiri et al., 2010; 840 individuals were reared in Rebaudo et al., 2017), provided that a high number of temperatures inside the range of critical thresholds are tested (fig. 3de). However, when the coefficient of variation is high (CV = 50%), the required rearing effort for an accurate characterization of CT max and T opt can be high and potentially technically impossible to achieve (e.g. with a uniform distribution of temperatures inside the range of critical thresholds, a total of 7500 individuals were required to obtain values of mean absolute errors lower than 1°C).

In our simulations we fitted seven alternative models in order to estimate the impacts of misspecification. Depending on the biological trait and on the fitted model, positive and negative biases were observed (fig. 4). This result suggests that model selection can be biased by experimental effort, resulting in biased values of biological traits. The results suggest that T opt estimation is more robust to the fitted model than CT min and CT max . Bias for CT min estimations was lower in comparison with CT max (fig. 4), for which median error took values close to zero for two models (Brière-2 and Analytis). Thus, our results suggest that the estimation of CT max requires a significant effort in experimental design to obtain unbiased estimations.

Our two analysis on simulated and experimental dataset show the importance of sample size in the estimations of biological traits. Because the experimental dataset used is composed of mean development rates, we could not evaluate the role of the variance of development rate on the estimation of parameters. The simulations show that the accuracy of estimations depends on the coefficient of variation. Both analysis suggested that the estimation of biological traits (attributes of development rate thermal performance curve) may include uncertainty. Also, the results of the simulations suggest that the lack of development rate data near the critical thresholds could lead to the observed discrepancies in estimated values of attributes of development rate thermal performance curve.

If insect response to temperature has been studied until now mostly by assessing lower critical thresholds (Kingsolver and Buckley, 2020;Stejskal et al., 2019)