# Modeling Temperature-Dependent Development Rate in Insects and Implications of Experimental Design

Baptiste Régnier, Judith Legrand, François Rebaudo

<figure><img src="image_1.png" /></figure>

► To cite this version:

Baptiste Régnier, Judith Legrand, François Rebaudo. Modeling Temperature-Dependent Development Rate in Insects and Implications of Experimental Design. Environmental Entomology, 2022, 51 (1), pp.132-144. 10.1093/ee/nvab115. hal-03982548

HAL Id: hal-03982548

https://agroparistech.hal.science/hal-03982548

Submitted on 30 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
**Title:** Modeling temperature-dependent development rate in insects and implications of experimental design

**Authors' names and addresses:** Baptiste Régnier<fnref n="1" />, Judith Legrand<fnref n="2" />, François Rebaudo<fnref n="1" />*

**Abstract**

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.

**Keywords:** development rate; laboratory; nonlinear models; temperature

**1. 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 (Ewing et al., 2016; Haridas et al., 2016; Nielsen et al., 2016). Insect phenological models have received much attention in entomology (Chuine and Régnière, 2017, Shi et al., 2017b). Used primarily for pest management (e.g. Crespo-Pérez et al., 2011; Rebaudo et

<footnote n="1">¹ Université Paris-Saclay, CNRS, IRD, UMR Évolution, Génomes, Comportement et Écologie, 91190, Gif-sur-Yvette, France.</footnote>

<footnote n="2">² Université Paris-Saclay, INRAE, CNRS, AgroParisTech, GQE - Le Moulon, 91190, Gif-sur-Yvette, France.</footnote>

<footnote>*Corresponding author: francois.rebaudo@ird.fr ; tel. +33 (0)1 69 82 37 48 ; IRD - UMR 247 EGCE, 1 Avenue de la Terrasse, bât. 13, BP1, 91190 Gif-sur-Yvette cedex, France.</footnote>
al., 2017; Roy et al., 2002), 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 (Shi et al. 2017b). 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; Huey and Stevenson, 1979; 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<sub>min</sub> and CT<sub>max</sub>) and the optimal temperature for development (T<sub>opt</sub>). 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).

<figure><img src="image_3.png" /><figcaption>Fig. 1: Theoretical performance curve for temperature-dependent development rate. Models generally assume that development rate is null under a critical temperature (CT<sub>min</sub>) from which it increases slowly as temperature increases, reaches a temperature zone where the relationship can be considered as linear (Campbell et al. 1974) until the development rate reaches a maximum at an optimal temperature (T<sub>opt</sub>). It then rapidly decreases and is considered null under a critical temperature (CT<sub>max</sub>).</figcaption></figure>

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
individuals), despite known limitations on the use of constant versus fluctuating
temperatures (Colinet et al., 2015). 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; Roy et al., 2002). 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 (Campbell et al., 1974). 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; Shi et al., 2012). 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. Briere et al., 1999; Lactin et
al., 1995; Sharpe and DeMichele, 1977), 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 Damos and Savopoulou-Soultani, 2012; Kontodimas et al., 2004; Quinn,
2017; Régnière et al., 2012; Shi and Ge, 2010; Shi et al., 2016, 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², 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.

## 2. 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*<sub>j</sub> the number of insects reared at temperature T<sub>j</sub> (j=1, ..., J) and D<sub>ij</sub> the length of the stage for the ith individual reared at temperature T<sub>j</sub>. For each individual, the development rate is computed as the inverse of the length of the stage r<sub>ij</sub>=1/D<sub>ij</sub>. Below, we denote as r<sub>ij</sub> the observed values of the development rate of individual i at temperature T<sub>j</sub> and as r̂<sub>j</sub> the prediction from the model r(T<sub>j</sub>).

### 2.1. 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.

### 2.1.1. *In silico* biological and experimental scenarios

All *in silico* biological scenarios were simulated with the Brière-2 model (Briere et al., 1999), chosen for its wide acceptance in the literature (e.g. Aghdam et al., 2009; Roy et al., 2002; Tochen et al., 2014).

$$r(T)=aT(T-CT_{min})(CT_{max}-T)^{(1/b)} \quad (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} = \frac{(2 b CT_{max} + (b+1) CT_{min}) + \sqrt{4 b^2 CT_{max}^2 + (b+1)^2 CT_{min}^2 - 4 b^2 CT_{min} CT_{max}}}{4 b + 2} \quad (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 Aghdam et al. (2009) ($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 (Régnière et al., 2012; Régnière, 1984; 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 $\mu_J$ and $\sigma^2$ where $\mu_J = \log(d_J) - (\sigma^2/2)$:

$$\log(D_{ij}) \sim N(\log(d_j) - (\sigma^2/2), \sigma^2) \quad (3)$$

The variance of the development time is given by $\sigma_{d_j}^2 = (e^{\sigma^2} - 1)d_j^2$ and the coefficient of variation by $cv_j = \sigma_{d_j}/d_j = \sqrt{(e^{\sigma^2} - 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
*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 $\sigma^2$ 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 (Aghdam et al., 2009; Briere and Pracros, 1998; Howell and Neven, 2000; Pollard et al., 2020; Roy et al., 2002; 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; Roy et al., 2002; 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. Campbell et al., 1974; Howell and Neven, 2000; 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.

### 2.1.2. 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(\hat{p}_{kz}) = \hat{p}_{kz} - p_k \quad (4) $$

where $\hat{p}_{kz}$ 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.

### 2.1.3. 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 Analytis, 1977; Performance model from Shi et al., 2011; Shi model from Shi et al., 2011; Brière-1 model from Brière et al., 1999; 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).

## 2.2. 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 (Davidson 1944, Sharpe and DeMichele 1977, Scoolfield et al. 1981, Brière 1999, Ikemoto 2005, Régnière et al. 2012). 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.

### 2.2.1. 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 Got and Rodolphe, 1989; Hansen et al., 2011; Régnière et al., 2012; 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 (Björck, 1996) 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 (Elzhov et al., 2016).

### 2.2.2. 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 ($\hat{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; Shi and Ge, 2010). 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 = \sum_{j=1}^{J} \sum_{i=1}^{n_j} (r_{ij} - \hat{r}_j)^2 \qquad (5)$$

The root-mean squared error (RMSE; equation 2) is defined as:

$$RMSE = \sqrt{\frac{RSS}{n}} \qquad (6)$$

Its usage is similar to RSS for characterizing insect development rate model accuracy (e.g. Damos and Savopoulou-Soultani, 2008; Mirhosseini et al., 2017; Shi et al., 2016). 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; Akaike, 1974) and Bayesian Information Criterion (BIC;
Schwarz, 1978) 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 (Burnham and Anderson, 2002), the AIC and BIC can be calculated with:

$$AIC = n \times \log\left(\frac{RSS}{n}\right) + 2K \quad (7)$$

$$BIC = n \times \log\left(\frac{RSS}{n}\right) + \log|n| \times K \quad (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. Arbab and Mcneill, 2011).

According to Burnham and Anderson (2002) 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.

### 2.2.3. 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 CTmin and CTmax, notably when data are sparse in low and high temperatures (Ikemoto and Kiritani, 2019). 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 CTmin inferior to 5°C, based on measured survival rates in low temperatures (Crespo-Pérez et al., 2011). 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 temperature-
dependent 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 looks for the parameter estimations that minimize the residual sum of squares
in respect of the constraint (e.g. $CT_{min} \le 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; Sharpe and DeMichele, 1977; 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*<sub>min</sub> and *CT10*<sub>max</sub>.

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 Crespo-Pérez et al. (2011). The dataset is composed of 14 observed mean development rates measured at 9 different temperatures (see Crespo-Pérez et al., 2011). *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 (Rebaudo et al., 2018). The complete code is available in Supplementary Material 1. The dataset for *T. solanivora* is included in the package *devRate*.

# 3. Results

## 3.1 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*<sub>max</sub> and *T*<sub>opt</sub> had higher interquartile ranges of errors than *CT*<sub>min</sub> (from 0.1204°C to 1137.27°C for *CT*<sub>max</sub>; from 0.0724°C to 4.44°C for *CT*<sub>min</sub>; from 0.0421°C to 1132.67°C for *T*<sub>opt</sub> 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*<sub>min</sub> was equal to 5.14°C. With this scenario resulting in a lack of data in high temperatures, *T*<sub>opt</sub> and *CT*<sub>max</sub> estimations took extreme values resulting in high error values superior to 10<sup>4</sup>°C.

### 3.1.2. 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*<sub>max</sub> 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*<sub>max</sub> 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*<sub>max</sub> 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*<sub>opt</sub> estimations, but the values of MAE were lower than for *CT*<sub>max</sub> estimations (fig. 3de). For *CT*<sub>min</sub>, 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*<sub>min</sub>, 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*<sub>max</sub> 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*<sub>max</sub> 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*<sub>max</sub> MAE < 1°C (fig. 3d). The minimal value of *CT*<sub>max</sub> 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*<sub>max</sub> MAE < 1°C. With a uniform distribution, at least the double of individuals were required to obtain *CT*<sub>max</sub> 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*<sub>max</sub> 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 CTmax MAE found (MAE = 5.79°C) was for 15 temperatures and 500 individuals per temperature.
<table><caption>**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. *CT*<sub>min</sub>, *CT*<sub>max</sub>, *T*<sub>opt</sub>, *CT*10<sub>min</sub> and *CT*10<sub>max</sub> 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.</caption><thead><tr><th>Model name</th><th>Reference</th><th>K</th><th>RSS</th><th>RMSE</th><th>AIC</th><th>BIC</th><th>CT<sub>min</sub></th><th>CT<sub>max</sub></th><th>T<sub>opt</sub></th><th>CT10<sub>min</sub></th><th>CT10<sub>max</sub></th></tr></thead><tbody><tr><td>janisch_32</td><td>Janisch (1932)</td><td>4</td><td>0.0012</td><td>0.0093</td><td>-121</td><td>-117.8</td><td></td><td></td><td>23.5</td><td>2.8</td><td>28.2</td></tr><tr><td>schoolfield_81***</td><td>Schoolfield et al. (1981)</td><td>6</td><td>0.0011</td><td>0.0087</td><td>-118.8</td><td>-114.4</td><td></td><td></td><td>27</td><td>6.7</td><td>29.3</td></tr><tr><td>wang_82***</td><td>Wang et al. (1982)</td><td>6</td><td>0.0011</td><td>0.0088</td><td>-118.6</td><td>-114.1</td><td>7 (13.3)</td><td></td><td>24.4</td><td>8.4</td><td>28.1</td></tr><tr><td>regniere_12***</td><td>Regniere et al. (2012)</td><td>6</td><td>0.0011</td><td>0.0088</td><td>-118.5</td><td>-114.1</td><td>9.9 (1735.4)</td><td>30.1 (0.54)</td><td>23.5</td><td>9.9</td><td>29.8</td></tr><tr><td>taylor_81***</td><td>Taylor (1981)</td><td>3</td><td>0.0020</td><td>0.0118</td><td>-116.3</td><td>-113.8</td><td></td><td></td><td>21.4 (0.73)</td><td>9.5</td><td>33.4</td></tr><tr><td>poly4</td><td>-</td><td>5</td><td>0.0015</td><td>0.0103</td><td>-116</td><td>-112.2</td><td></td><td>31.1</td><td>21.6</td><td></td><td>30.2</td></tr><tr><td>ratkowsky_83***</td><td>Ratkowsky et al. (1983)</td><td>4</td><td>0.0018</td><td>0.0112</td><td>-115.7</td><td>-112.5</td><td>5 (3.47)</td><td>33.5 (0.91)</td><td>22.2</td><td>8.6</td><td>31.9</td></tr><tr><td>damos_11</td><td>Damos and Savopoulou (2011)</td><td>3</td><td>0.0022</td><td>0.0126</td><td>-114.5</td><td>-112</td><td></td><td></td><td>21.8</td><td>4.1</td><td>39.4</td></tr><tr><td>analytis_77***</td><td>Analytis (1977)</td><td>5</td><td>0.0018</td><td>0.0115</td><td>-113.1</td><td>-109.3</td><td>4.5 (21.8)</td><td>35.1 (21.1)</td><td>21.7</td><td>9.7</td><td>32.1</td></tr><tr><td>poly2***</td><td>-</td><td>3</td><td>0.0027</td><td>0.0138</td><td>-112</td><td>-109.5</td><td>8.6</td><td>34.2</td><td>21.4</td><td>9.2</td><td>33.6</td></tr><tr><td>kontodimas_04***</td><td>Kontodimas et al. (2004)</td><td>3</td><td>0.0028</td><td>0.0142</td><td>-111.2</td><td>-108.7</td><td></td><td>34.1 (1.09)</td><td>23.3</td><td>6.1</td><td>33.6</td></tr><tr><td>damos_08</td><td>Damos and Savopoulou (2008)</td><td>3</td><td>0.0029</td><td>0.0143</td><td>-110.9</td><td>-108.4</td><td>0</td><td>34.1</td><td>23.3</td><td>5</td><td>33.6</td></tr><tr><td>harcourtYee_82</td><td>Harcourt and Yee (1982)</td><td>4</td><td>0.0026</td><td>0.0137</td><td>-110.1</td><td>-106.9</td><td>8.2</td><td>34.1</td><td>22</td><td>9</td><td>33.5</td></tr><tr><td>perf2_11</td><td>Shi et al. (2011)</td><td>4</td><td>0.0026</td><td>0.0137</td><td>-110.1</td><td>-106.9</td><td>8.3 (1.86)</td><td>34.1 (1.42)</td><td>21.8</td><td>9.1</td><td>33.5</td></tr><tr><td>lactin2_95</td><td>Lactin et al. (1995)</td><td>4</td><td>0.0027</td><td>0.0139</td><td>-109.8</td><td>-106.6</td><td>7.6</td><td>34.1</td><td>22.6</td><td>8.5</td><td>33.6</td></tr><tr><td>logan10_76</td><td>Logan et al. (1976)</td><td>5</td><td>0.0024</td><td>0.0132</td><td>-109.3</td><td>-105.4</td><td>6.5</td><td>34 (3.3)</td><td>21.8</td><td>8.4</td><td>33.4</td></tr><tr><td>beta_95</td><td>Yin et al. (1995)</td><td>5</td><td>0.0026</td><td>0.0136</td><td>-108.3</td><td>-104.5</td><td>8.5 (7.25)</td><td>35.1 (28.6)</td><td>21</td><td>9.2</td><td>33.9</td></tr><tr><td>hilbertLogan_83</td><td>Hilbert and Logan (1983)</td><td>5</td><td>0.0026</td><td>0.0136</td><td>-108.3</td><td>-104.5</td><td>8</td><td>34.1</td><td>21.9</td><td>8.9</td><td>33.5</td></tr><tr><td>shi_11</td><td>Shi et al. (2011)</td><td>5</td><td>0.0026</td><td>0.0137</td><td>-108.1</td><td>-104.3</td><td>8.3 (2.83)</td><td>34.1 (1.58)</td><td>21.8</td><td>9.1</td><td>33.5</td></tr><tr><td>briere1_99</td><td>Briere et al. (1999)</td><td>3</td><td>0.0041</td><td>0.0170</td><td>-106</td><td>-103.5</td><td></td><td>35.1 (0.01)</td><td>25.3</td><td>1.9</td><td>35</td></tr><tr><td>lactin1_95</td><td>Lactin et al. (1995)</td><td>3</td><td>0.0042</td><td>0.0174</td><td>-105.5</td><td>-103</td><td>9</td><td></td><td>18.4</td><td>9.4</td><td></td></tr><tr><td>logan6_76</td><td>Logan et al. (1976)</td><td>4</td><td>0.0038</td><td>0.0165</td><td>-104.9</td><td>-101.7</td><td></td><td>34.4 (1.17)</td><td>24.9</td><td></td><td>34.1</td></tr><tr><td>beta_16</td><td>Shi et al. (2016)</td><td>4</td><td>0.0038</td><td>0.0165</td><td>-104.9</td><td>-101.7</td><td></td><td>34.4 (1.16)</td><td>24.9</td><td></td><td>34.1</td></tr><tr><td>bayoh_03</td><td>Bayoh and Lindsay (2003)</td><td>4</td><td>0.0044</td><td>0.0177</td><td>-102.9</td><td>-99.8</td><td>9.5</td><td>35</td><td>16.3</td><td>9.6</td><td>34.9</td></tr><tr><td>davidson_44</td><td>Davidson (1944)</td><td>3</td><td>0.0072</td><td>0.0227</td><td>-98.1</td><td>-95.5</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>briere2_99</td><td>Briere et al. (1999)</td><td>4</td><td>0.0067</td><td>0.0218</td><td>-97.1</td><td>-93.9</td><td>1.1 (25.4)</td><td>35.1 (0.01)</td><td>30.4</td><td>7.7</td><td>35</td></tr></tbody></table>
### 3.1.3. 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*<sub>min</sub>, *CT*<sub>max</sub> and *T*<sub>opt</sub> when another model than Brière-2 is fitted. The simulations resulted in different values of median error for *CT*<sub>min</sub>, *T*<sub>opt</sub> and *CT*<sub>max</sub> estimations (fig. 4). For *CT*<sub>min</sub>, 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*<sub>max</sub>, 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*<sub>opt</sub>, 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.

<figure><img src="image_4.png" /><figcaption>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 *CT*10<sub>min</sub> and *CT*10<sub>max</sub>, 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*<sub>opt</sub>).</figcaption></figure>
## 3.2. 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*<sub>min</sub> estimations varied between 1.86°C and 1735.4°C; for *CT*<sub>max</sub>, 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 (Janisch, 1932). 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*<sub>min</sub> and *CT*<sub>max</sub> could not be computed for some models, we analyzed the values of *CT10*<sub>min</sub> and *CT10*<sub>max</sub>. Based on observed development rates and survival at low temperatures (Crespo-Pérez et al., 2011), the minimal temperature threshold where development rate is lower than 10% of the maximal development rate (*CT10*<sub>min</sub>) 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*<sub>min</sub> nor *CT10*<sub>min</sub> 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 (Analytis, 1977; Kontodimas et al., 2004; Ratkowsky et al., 1983; Regniere et al., 2012; Schoolfield et al., 1981; 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*<sub>opt</sub>) 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*<sub>opt</sub> values between 21.4°C and 24.4°C. *CT10*<sub>min</sub> varied between 6.1°C with the Kontodimas model and 9.9°C with the Régnière model. *CT10*<sub>max</sub> varied between 28.1°C with the Wang model and 33.6°C with the Kontodimas model.
<figure><img src="image_5.png" /><figcaption>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<sub>min</sub> 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<sub>max</sub> estimations, e) same results for T<sub>opt</sub> estimations.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 4: Bias for CT<sub>min</sub>, T<sub>opt</sub> and CT<sub>max</sub> estimated with eight models, including seven misspecified models. Bias was measured with median errors for three biological traits (CT<sub>min</sub>, T<sub>opt</sub> and CT<sub>max</sub>) 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.</figcaption></figure>

**4. 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 (Crespo-Pérez et
al., 2011). 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<sub>min</sub> > 5°C), the best fitting model according to
statistical criteria (Janisch model) was discarded. This result suggests that statistical
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*<sub>opt</sub> 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* (Crespo-Pérez et al., 2011), 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* temperature-dependent 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 (Burnham and Anderson, 2002).

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 (Gray et al., 2001). 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 (Araujo
and New, 2007), 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 (Arlot
and Celisse, 2010) 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; Régnière et al., 2012; 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<sub>min</sub> estimations were always far lower than for CT<sub>max</sub> and
T<sub>opt</sub> estimations, which suggests that the required experimental design depends on the
pursued objectives. If the estimation of CT<sub>min</sub> 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<sub>min</sub> 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<sub>min</sub> are not required in order to estimate
accurately its value with the nonlinear Brière-2 model (Briere et al., 1999), as the
distribution of temperatures in the linear zone resulted in low mean absolute errors of
CT<sub>min</sub> (from 1.22°C to 0.1°C with CV = 10%). However, if CT<sub>max</sub> and T<sub>opt</sub> 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<sub>opt</sub> and CT<sub>max</sub>). In the absence of data in this zone, the errors for the
estimations of CT<sub>max</sub> and T<sub>opt</sub> could take values higher than 10<sup>4</sup>°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<sub>max</sub> 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 mid-temperatures, 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
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 CTmax and Topt 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<sub>opt</sub> estimation is more robust to the fitted model than CT<sub>min</sub> and CT<sub>max</sub>. Bias for CT<sub>min</sub> estimations was lower in comparison with CT<sub>max</sub> (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<sub>max</sub> 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), insect response to high temperatures is an essential component to study development rate in the context of global warming (Ma et al., 2021). Thus, accurate characterizations of CTmax values for temperature-dependent development rate are necessary for a better understanding of the influence of high temperature events on the development of
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.

## Acknowledgments

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

## References

Aghdam, H.R., Y. Fathipour, G. Radjabi, and M. Rezapanah. 2009. Temperature-dependent development and temperature thresholds of codling moth (Lepidoptera: Tortricidae) in Iran. Environ. Entomol. 38:885-895.

Akaike, H. 1974. A new look at the statistical model identification. IEEE Trans. Automat. Contr. 19:716-723.

Analytis, S. 1977. Über die relation zwischen biologischer entwicklung und temperatur bei phytopathogenen pilzen. J. Phytopathol. 90:64-76.

Araujo, M., and M. New. 2007. Ensemble forecasting of species distributions. Trends Ecol. Evol. 22:42-47.

Arbab, A., and M.R. Mcneill. 2011. Determining suitability of thermal development models to estimate temperature parameters for embryonic development of *Sitona lepidus* Gyll. (Coleoptera: Curculionidae). J. Pest. Sci. 84:303-311.

Arlot, S., and A. Celisse. 2010. A survey of cross-validation procedures for model selection. Stat. Surv. 4:40-79.

Björck, Å., 1996. Numerical methods for least squares problems. Society for

Industrial and Applied Mathematics, Philadelphia, Pennsylvania, United States.

Briere, J.-F., and P. Pracros. 1998. Comparison of temperature-dependent growth models with the development of *Lobesia botrana* (Lepidoptera: Tortricidae). Environ. Entomol. 27:94-101.

Briere, J.-F., P. Pracros, A.-Y. Le Roux, and J.-S. Pierre. 1999. A novel rate model of temperature-dependent development for arthropods. Environ. Entomol. 28:22-29.

Burnham, K.P., and D.R. Anderson. 2002. Model selection and multimodel inference: a practical information-theoretic approach, 2nd ed. Springer, New York.

Campbell, A., B.D. Frazer, N. Gilbert, A.P. Gutierrez, and M. Mackauer. 1974. Temperature requirements of some aphids and their parasites. J. Appl. Ecol. 11:431-438.

Chuine, I., and J. Régnière. 2017. Process-based models of phenology for plants and animals. Annu. Rev. Ecol. Evol. Syst. 48:159-182.
Colinet, H., B.J. Sinclair, P. Vernon, and D. Renault. 2015. Insects in fluctuating thermal environments. Annu. Rev. Entomol. 60:123-140.

Crespo-Pérez, V., F. Rebaudo, J.-F. Silvain, and O. Dangles. 2011. Modeling invasive species spread in complex landscapes: the case of potato moth in Ecuador. Landsc. Ecol. 26:1447-1461.

Damos, P., and M. Savopoulou-Soultani. 2012. Temperature-driven models for insect development and vital thermal requirements. Psyche. 2012:123405.

Damos, P., and M. Savopoulou-Soultani. 2008. Temperature-dependent bionomics and modeling of *Anarsia lineatella* (Lepidoptera: Gelechiidae) in the laboratory. J. Econ. Entomol. 101:1557-1567.

Davidson, J., 1944. On the relationship between temperature and rate of development of insects at constant temperatures. J. Anim. Ecol. 13:26-38.

Elzhov, T.V., K.M. Mullen, A.-N. Spiess, and B. Bolker. 2016. minpack.lm: R interface to the Levenberg-Marquardt nonlinear least-squares algorithm found in minpack, plus support for bounds.

Ewing, D.A., C.A. Cobbold, B.V. Purse, M.A. Nunn, and S.M. White. 2016. Modelling the effect of temperature on the seasonal population dynamics of temperate mosquitoes. J. Theor. Biol. 400:65-79.

Got, B., and F. Rodolphe. 1989. Temperature-dependent model for European corn borer (Lepidoptera: Pyralidae) development. Environ. Entomol. 18:85-93.

Gray, D.R., F.W. Ravlin, and J.A. Braine. 2001. Diapause in the gypsy moth: a model of inhibition and development. J. Insect Physiol. 47:173-184.

Hansen, E.M., B.J. Bentz, J.A. Powell, D.R. Gray, and J.C. Vandygriff. 2011. Prepupal diapause and instar IV developmental rates of the spruce beetle, *Dendroctonus rufipennis* (Coleoptera: Curculionidae, Scolytinae). J. Insect Physiol. 57:1347-1357.

Haridas, C.V., L.J. Meinke, B.E. Hibbard, B.D. Siegfried, and B. Tenhumberg. 2016. Effects of temporal variation in temperature and density dependence on insect population dynamics. Ecosphere. 7:e01287.

Howell, J.F., and L.G. Neven. 2000. Physiological development time and zero development temperature of the codling moth (Lepidoptera: Tortricidae). Environ. Entomol. 29:766-772.

Huey, R.B., and R.D. Stevenson. 1979. Integrating thermal physiology and ecology of ectotherms: a discussion of approaches. Am. Zool. 19:357-366.

Ikemoto, T., 2005. Intrinsic optimum temperature for development of insects and mites. Environ. Entomol. 34:1377-1387.

Ikemoto, T., and K. Kiritani. 2019. Novel method of specifying low and high threshold temperatures using thermodynamic SSI model of insect development. Environ. Entomol. 48:479-488.

Janisch, E., 1932. The influence of temperature on the life-history of
insects. Trans. R. Entomol. Soc. Lond.
80:137-168.

Kingsolver, J.G., and L.B. Buckley. 2020.
Ontogenetic variation in thermal
sensitivity shapes insect ecological
responses to climate change. Curr.
Opin. Insect. Sci. 41:17-24.

Kontodimas, D.C., P.A. Eliopoulos, G.J.
Stathas, and L.P. Economou. 2004.
Comparative temperature-dependent
development of *Nephus includens*
(Kirsch) and *Nephus bisignatus*
(Boheman) (Coleoptera: Coccinellidae)
preying on *Planococcus citri* (Risso)
(Homoptera: Pseudococcidae):
evaluation of a linear and various
nonlinear models using specific
criteria. Environ. Entomol. 33:1-11.

Kreppel, K.S., S. Telfer, M. Rajerison, A.
Morse, and M. Baylis. 2016. Effect of
temperature and relative humidity on
the development times and survival
of *Synopsyllus fonquerniei* and
*Xenopsylla cheopis*, the flea vectors of
plague in Madagascar. Parasit.
Vectors. 9:1-10.

Lactin, D.J., N.J. Holliday, D.L. Johnson, and R. Craigen. 1995. Improved rate model of temperature-dependent development by arthropods. Environ. Entomol. 24:68-75.

Lamb, R., G. Gerber, and G. Atkinson. 1984.
Comparison of developmental rate
curves applied to egg hatching data of
*Entomoscelis americana* Brown
(Coleoptera: Chrysomelidae). Environ.
Entomol. 13:868-872.

Levenberg, K. 1944. A method for the
solution of certain non-linear
problems in least squares. Q. Appl.
Math. 2:164-168.

Logan, J.A., D.J. Wollkind, S.C. Hoyt, and L.K.
Tanigoshi. 1976. An analytic model for
description of temperature
dependent rate phenomena in
arthropods. Environ. Entomol.
5:1133-1140.

Ma, C.-S., G. Ma, and S. Pincebourde. 2021.
Survive a warming climate: insect
responses to extreme high
temperatures. Annu. Rev. Entomol.
66:163-184.

McManis, A., J. Powell, and B. Bentz. (2018).
Developmental parameters of a
southern mountain pine beetle
(Coleoptera: Curculionidae)
population reveal potential source of
latitudinal differences in generation
time. Can. Entomol., 151:1-15.

Marquardt, D.W. 1963. An algorithm for
least-squares estimation of nonlinear
parameters. J. Soc. Ind. Appl. Math.
11:431-441.

Martin, T.L., and R.B. Huey. 2008. Why "suboptimal" is optimal: Jensen's inequality and ectotherm thermal preferences. Am. Nat. 171:E102-E118.

Martín-Vega, D., L.M. Díaz-Aranda, A. Baz, and B. Cifrián. 2017. Effect of temperature on the survival and development of three forensically relevant *Dermestes* species (Coleoptera: Dermestidae). J. Med. Entomol. 54:1140-1150.

Messenger, P.S., and N.E. Flitters. 1958.
Effect of constant temperature
environments on the egg stage of
three species of Hawaiian fruit flies.
Ann. Entomol. Soc. Am. 51:109-119.

Mirhosseini, M.A., Y. Fathipour, and G.V.P.
Reddy. 2017. Arthropod
development's response to
temperature: a review and new software for modeling. Ann. Entomol. Soc. Am. 110:507-520.

Mwalusepo, S., H.E.Z. Tonnang, E.S. Massawe, G.O. Okuku, N. Khadioli, T. Johansson, T., P.-A. Calatayud, and B.P. Le Ru. 2015. Predicting the impact of temperature change on the future distribution of maize stem borers and their natural enemies along east african mountain gradients using phenology models. PLoS One. 10:e0130427.

Nielsen, A.L., S. Chen, and S.J. Fleischer. 2016. Coupling developmental physiology, photoperiod, and temperature to model phenology and dynamics of an invasive heteropteran, *Halyomorpha halys*. Front. Physiol. 7:165.

Nocedal, J., and S.J. Wright. 2006. Numerical optimization, 2nd ed. Springer, New York, NY

Noor-ul-Ane, M., M. Ali Mirhosseini, N. Crickmore, S. Saeed, I. Noor, and M.P. Zalucki. 2018. Temperature-dependent development of *Helicoverpa armigera* (Hübner) (Lepidoptera: Noctuidae) and its larval parasitoid, *Habrobracon hebetor* (Say) (Hymenoptera: Braconidae): implications for species interactions. Bull. Entomol. Res. 108:295-304.

Orang, F.S., H.R. Aghdam, H. Abbasipour, and A. Askarianzadeh. 2014. Effect of temperature on developmental rate of *Sesamia cretica* (Lepidoptera: Noctuidae) immature stages. J. Insect Sci. 14:197.

IPCC (Intergovernmental Panel on Climate Change). 2014. Climate change 2014:

synthesis report. Contribution of Working Groups I, II, and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Geneva, Switzerland: IPCC.

Pincebourde, S., and J. Casas. 2015. Warming tolerance across insect ontogeny: influence of joint shifts in microclimates and thermal limits. Ecology. 96:986-997.

Pollard, C.P., C.T. Griffin, R. Andrade Moral, C. Duffy, J. Chuche, M.T. Gaffney, R.M. Fealy, and R. Fealy. 2020. phenModel: A temperature-dependent phenology/voltinism model for a herbivorous insect incorporating facultative diapause and budburst. Ecol. Modell. 416:108910.

Quinn, B.K. 2017. A critical review of the use and performance of different function types for modeling temperature-dependent development of arthropod larvae. J. Therm. Biol. 63:65-77.

R Core Team. 2021. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.

Ratkowsky, D.A., R.K. Lowry, T.A McMeekin, A.N. Stokes, and R.E. Chandler. 1983. Model for bacterial culture growth rate throughout the entire biokinetic temperature range. J. Bacteriol. 154:1222-1226.

Rebaudo, F., and V.-B. Rabhi. 2018. Modeling temperature-dependent development rate and phenology in insects: review of major developments, challenges, and future directions. Entomol. Exp. Appl. 166:607-617.
Rebaudo, F., Q. Struelens, F. Callizaya Condori, and R. Quispe. 2017. Relationship between temperature and development rate of *Copitarsia incommoda* (Lepidoptera: Noctuidae) in the Bolivian Andes. Appl. Entomol. Zool. 52:313-320.

Rebaudo, F., Q. Struelens, and O. Dangles. 2018. Modelling temperature-dependent development rate and phenology in arthropods: The devRate package for R. Methods Ecol. Evol. 9:1144-1150.

Régnière, J., 1984. A method of describing and using variability in development rate for the simulation of insect phenology. Can. Entomol. 116:1367-1376.

Régnière, J., J. Powell, B. Bentz, and V. Nealis. 2012. Effects of temperature on development, survival and reproduction of insects: experimental design, data analysis and modeling. J. Insect Physiol. 58:634-647.

Roy, M., J. Brodeur, and C. Cloutier. 2002. Relationship between temperature and developmental rate of *Stethorus punctillum* (Coleoptera: Coccinellidae) and its prey *Tetranychus mcdanieli* (Acarina: Tetranychidae). Environ. Entomol. 31:177-187.

Schoolfield, R.M., P.J.H. Sharpe, and C.E. Magnuson. 1981. Non-linear regression of biological temperature-dependent rate models based on absolute reaction-rate theory. J. Theor. Biol. 88:719-731.

Schwarz, G. 1978. Estimating the dimension of a model. Ann. Stat. 6:461-464.

Sharpe, P.J.H., and D.W. DeMichele. 1977. Reaction kinetics of poikilotherm development. J. Theor. Biol. 64:649-670.

Shi, P., and F. Ge. 2010. A comparison of different thermal performance functions describing temperature-dependent development rates. J. Therm. Biol. 35:225-231.

Shi, P., F. Ge, Y. Sun, and C. Chen. 2011. A simple model for describing the effect of temperature on insect developmental rate. J. Asia Pac. Entomol. 14:15-20.

Shi, P., B. Wang, M.P. Ayres, F. Ge, L. Zhong, and B.-L. Li 2012. Influence of temperature on the northern distribution limits of *Scirpophaga incertulas* Walker (Lepidoptera: Pyralidae) in China. J. Therm. Biol. 37:130-137

Shi, P.-J., G.V.P. Reddy, L. Chen, and F. Ge. 2016. Comparison of thermal performance equations in describing temperature-dependent developmental rates of insects: (I) empirical models. Ann. Entomol. Soc. Am. 109:211-215.

Shi, P.-J., G.V.P. Reddy, L. Chen, and F. Ge. 2017a. Comparison of thermal performance equations in describing temperature-dependent developmental rates of insects: (II) two thermodynamic models. Ann. Entomol. Soc. Am. 110:113-120.

Shi, P.-J., M.-L. Fan, and G.V.P. Reddy. 2017b. Comparison of thermal performance equations in describing temperature-dependent developmental rates of insects: (III)
phenological applications. Ann. Entomol. Soc. Am. 110:558-564.

Shi, P., D.A. Ratkowsky, N. Wang, Y. Li, L. Zhao, G.V.P. Reddy, and B.-L. Li. 2017c. Comparison of five methods for parameter estimation under Taylor's power law. *Ecol. Complex.* 32:121-130.

Sinclair, B.J., K.E. Marshall, M.A. Sewell, D.L. Levesque, C.S. Willett, S. Slotsbo, Y. Dong, C.D.G. Harley, D.J. Marshall, B.S. Helmuth, and R.B. Huey. 2016. Can we predict ectotherm responses to climate change using thermal performance curves and body temperatures? *Ecol. Lett.* 19:1372-1385.

Stejskal, V., T. Vendl, Z. Li, and R. Aulicky. 2019. Minimal thermal requirements for development and activity of stored product and food industry pests (Acari, Coleoptera, Lepidoptera, Psocoptera, Diptera and Blattodea): A review. *Insects*. 10:149.

Stinner, R.E., G.D. Butler, J.S. Bachelier, and C. Tuttle. 1975. Simulation of temperature-dependent development in population dynamics model. *Can. Entomol.* 107:1167-1174.

Taylor, F. 1981. Ecology and evolution of physiological time in insects. *Am. Nat.* 117:1-23.

Tochen, S., D.T. Dalton, N. Wiman, C. Hamm, P.W. Shearer, and V.M. Walton. 2014. Temperature-related development and population parameters for *Drosophila suzukii* (Diptera: Drosophilidae) on cherry and blueberry. *Environ. Entomol.* 43:501-510.

Tran, L.T., S.P. Worner, R.J. Hale, and D.A.J. Teulon. 2012. Estimating development rate and thermal requirements of *Bactericera cockerelli* (Hemiptera: Triozidae) reared on potato and tomato by using linear and nonlinear models. *Environ. Entomol.* 41:1190-1198.

Ungerer, M.J., M.P. Ayres, and M.J. Lombardero. 1999. Climate and the northern distribution limits of *Dendroctonus frontalis* Zimmermann (Coleoptera: Scolytidae). *J. Biogeogr.* 26:1133-1145.

Uvarov, B.P. 1931. Insects and climate. *Trans. R. Entomol. Soc. Lond.* 79:1-232.

Walther, B.A., and J.L. Moore. 2005. The concepts of bias, precision and accuracy, and their use in testing the performance of species richness estimators, with a literature review of estimator performance. *Ecography.* 28:815-829.

Wang, R., Z. Lan, and Y. Ding. 1982. Studies on mathematical models of the relationship between insect development and temperature. *Acta Ecol. Sin.* 2:47-57.

Yurk, B.P., and J.A. Powell. 2010. Modeling the effects of developmental variation on insect phenology. *Bull. Math. Biol.* 72:1334-1360.

Zahiri, B., Y. Fathipour, M. Khanjani, S. Moharramipour, and M.P. Zalucki. 2010. Preimaginal development response to constant temperatures in *Hypera postica* (Coleoptera: Curculionidae): picking the best model. *Environ. Entomol.* 39:177-189.
Zuim, V., H.S. Rodrigues, D. Pratissoli, and J.B. Torres. 2017. Thermal requirements and performance of the parasitoid *Trichogramma pretiosum* (Hymenoptera: Trichogrammatidae)

on *Helicoverpa armigera* (Lepidoptera: Noctuidae) eggs under variable temperatures. Environ. Entomol. 46:1156-1164.

**Supplementary Material**

**Supplementary Material 1:** Characterizing temperature-dependent development rate in insects: model adjustment, assessment of goodness of fit, and selection using the R package devRate and immature development of *Tecia solanivora* in function of temperature.

**Supplementary Material 2:** Impact of experimental design on temperature-dependent development rate mathematical characterization, a simulation study: R script.