Jean-Sébastien Pierre

Maurice Hullé

Jean-Pierre Gauthier

Claude Rispe

Critical windows: A method for detecting lagged variables in ecological time series

Keywords: critical periods, correlation, time series, climate, lag analysis

Many developmental processes in the life sciences, ecology and even in economics depend strongly on the environmental conditions occurring in a bounded time interval, the results occurring often far later. Examples are as diverse as plant phenology, grapewine maturation, diapause induction and so on. The method proposed here, aims at detecting quickly such effects.

The basic idea is to regress the recorded results of a series of replications of the process against a function of an independent time series. This variable is defined on a set of periods of time, systematically scanned by varying their lower and upper bounds. In simple cases when this function is the integral and the effect strictly limited in a time window, the response model, under the form of correlation coefficients, is tractable and its shape is predictable. It is the same when the window is a bell-shaped function and can be fitted with other weighting functions as the Beta and the polynomials. The null hypothesis of absence of influence of any past interval is tested by Monte-Carlo simulation. The most likely window of influence is determined by the maximum correlation coefficient, and the bivariate confidence interval is estimated by bootstrap. The period found with a rectangular shaped window can be used as a starting point for more specific windows. This technique has the advantage of avoiding to split the climatic series into arbitrary slices, thus multiplying the predictors and complicating the models selection. It is closely linked to continuous lag distributed models with the simplification that the variable of interest is not explicitly time dependent. Examples are given for the prediction of aphids population dynamics, male morphs induction in aphids, and the phenology of the ash tree.

Introduction: What are critical periods ?

Some processes involving so called "critical periods" are fairly common in various fields of biology, ecology, behavioural sciences and psychology. They share the following features: a) the effect of a driving factor is lagged; b) the driving factor only has an effect during a limited period; c) the effect of the factor is cumulative during this period. The concept comes from Ethology and was coined by Lorenz [START_REF] Lorenz | Studies in animal and human behaviour[END_REF]: the strength of imprinting in birds and mammals, depends on the presence of the mother-like stimuli during a short period after hatching, and has no effect thereafter. This concept is now commonly used in developmental psychology.

Many examples can be drawn also from Ecology. For example, one can quote the proportion of insects entering diapause and the moment where they do so and, similarly the proportion of endodormancy in plants which depends upon environmental cues like the decrease of photoperiod and temperature in a given period of time. The research on plant phenology constitutes a large frame for such investigations. Similar relationships can be found between the amount of cold received during winter and the proportion of insects emerging from diapause. The number of sexual forms induced in various species of aphids depends on the decrease of temperature , as well as of photoperiod during a short period of the fall. The quality of grapes in vineyards is known to require rainfall at fruit set and a sunny weather at maturation. Many other examples can be found emphasizing the interest of modeling critical periods, which were studied since several decades under different names and with different approaches.

The present article aims at detecting statistically such phenomenon and at finding the range of time in which environmental conditions influence them. Practically, the procedure is not new and arises from at least two statistical approaches. One is the systematic scanning of time periods, and another is the use of lag distribution models.

Systematic scanning

The idea of scanning systematically the periods of influence of an external factor on a biological phenomenon was firstly applied to the phenology of vegetable crops, affected by weather, by Goldwin [9] . This author outlined the principles of the method, without studying in depth the underlying statistical model. The idea was then applied to the outbreaks of aphids by Thomas et al. [START_REF] Thomas | Associations between Weather Factors and the Spring Migration of the Damson-Hop Aphid, Phorodon humuli[END_REF] , with a first proposition of graphical representation, and then, in a quite different manner by Pierre et al [START_REF] Pierre | Une methode statistique et graphique de recherche des periodes de l'annee ou les populations animales sont particulierement sensibles a une composante donnee du climat (periodes critiques). Application au cas des pucerons des cereales[END_REF] . These latter authors made the first attempt to define the method as a bivariate correlogram, to propose an interpretable graphical representation of it, and to find its theoretical shape. Several cases of critical climatic periods were then found in various aphids: nine different species of winged aphids(A'brook, [START_REF] Rispe | Forecasting the incidence of aphids using weather data[END_REF]) , Phorodon humuli,(Thomas et al [START_REF] Thomas | Associations between Weather Factors and the Spring Migration of the Damson-Hop Aphid, Phorodon humuli[END_REF]), Sitobion avenae F.

(Pierre and Dedryver [START_REF] Pierre | Un modèle de prévision des pullulations du puceron Sitobion avenae F. sur blé d'hiver[END_REF], Pierre et al [START_REF] Pierre | Une methode statistique et graphique de recherche des periodes de l'annee ou les populations animales sont particulierement sensibles a une composante donnee du climat (periodes critiques). Application au cas des pucerons des cereales[END_REF] , Rhopalosiphum padi L., Rispe et al [START_REF] Rispe | Effect of climate on the proportion of males in the autumn flight of the aphid Rhopalosiphum padi L. (Hom., Aphididae)[END_REF]). The procedure suggested by Goldwin [START_REF] Goldwin | A technique for studying the association between components of the weather and horticultural parameters[END_REF] was also developped further by Coakley et al. [START_REF] Coakley | Model for predicting severity of Septoria tritici blotch on winter-wheat[END_REF] under the name "Window" and later by Pietravalle et al.. [START_REF] Pietravalle | Modeling of Relationships Between Weather and Septoria tritici Epidemics on Winter Wheat: A Critical Approach[END_REF] under the appellation "Window Pane". Both of them applied the method to the incidence of the wheat disease Septoria nodorum. Coakley suggested the boostrap method to find the variance and confidence interval of the best correlation coeffcient. Hullé and Gauthier rewrote the FORTRAN program written by Pierre et al. [START_REF] Pierre | Une methode statistique et graphique de recherche des periodes de l'annee ou les populations animales sont particulierement sensibles a une composante donnee du climat (periodes critiques). Application au cas des pucerons des cereales[END_REF] in language C, adding bootstrapping and Monte-Carlo simulation to test correctly the correlation peaks. The software application, named Criticor was then distributed freely by the French National

Institute of Agronomic Research (INRA, now INRAE). The program description and notice is

available on the site Yumpu [START_REF] Anonymous | CRITICOR, notice d'utilisation[END_REF]. This very method was then used by Ennaïfar et al [START_REF] Ennaïfar | Evaluation of models to predict take-all incidence in winter wheat as a function of cropping practices, soil, and climate[END_REF] for the take-all desease (Gaeumannomyces graminis var. tritici) incidence of wheat, and Atlan et al [START_REF] Atlan | Phenotypic Plasticity in Reproductive Traits of the Perennial Shrub Ulex europaeus in Response to Shading: A Multi-Year Monitoring of Cultivated Clones[END_REF] for the flowering of the european gorse (Ulex europaeus L.). Recently, Pierre (2017) deposited in CRAN a version written in R [START_REF] Core | R: A language and environment for statistical computing[END_REF] language named Rcriticor in reference to the INRA (now INRAE) version.

Statistically, this method has some advantages, especially that of reducing the number of possible predictors only to one for each type of climatic series. A widespread alternative is to split the past climate into arbitrary slices (weeks, months, decades,etc.) and to achieve a multiple regression of any kind (linear model, neuron network, learning,..) between the variable of interest and the set of predictors so created (e.g. A'Brook 1981 [START_REF] Rispe | Forecasting the incidence of aphids using weather data[END_REF], Tsai et al.2016 [START_REF] Tsai | Phenological responses of ash (Fraxinus excelsior) and sycamore (Acer pseudoplatanus) to riparian thermal conditions[END_REF], Ji and Peters [START_REF] Ji | Lag and Seasonality Considerations in Evaluating AVHRR NDVI Response to Precipitation[END_REF] among many others). Then comes the tedious work of selecting models, to disentangle the correlations between predictors. The research of critical periods does directly most of the work ahead, as we will show on a phenologic example. The problem of model selection is reduced therefore to a small number of predictors, one or at most two by type of climatic series. The correlations between them are furthermore reduced providing that the critical periods for diverse climatic series do not overlap.

As described above, the method that we propose is close to lag analysis, a statistical field which has recently attracted the interest of various scientists, especially in the field of plant phenology.

There is however a fundamental difference: the variable of interest is not time-dependent and is only supposed to be observed later that its predictor. At the contrary, in the lag distributed models, the predictor does not cease its action and influences the whole time series of the variable of interest in its future (Sims[22]). This difference has some statistical consequences that will be described in the Material and Methods section of this article.

Lag distributed models

The other approach, largely developped in the field of econometrics (e.g. Seddighi [START_REF] Seddighi | Introductory Econometrics: A Practical Approach[END_REF]) consists in the discrete or continuous [START_REF] Sims | Discrete approximations to continuous time distributed lags in econometrics[END_REF] lag-distributed models, sometime called briely lag analysis. In these models, the dependent variable is itself a time series or at least a variable indexed by time. The phenomenon observed at time t depends on several values of another time series, lagged in the past with each one its coefficient. This approach was recently used with success in Ecology.

Among the most elaborate approach using lag analysis, we can cite the work of Ogle et al..

(2015) [START_REF] Ogle | Quantifying ecological memory in plant and ecosystem processes[END_REF], who propose a comprehensive modelling of lag analysis in ecology, with a bayesian estimation of all relevant parameters. Their work is more ambitious than our one which aims only at a fast research to detect "hot spots" in a climatic series during the year. Although quite close to the preceeding approach, there is a subtle but noticeable difference as outlined above. So, [START_REF] Ji | Lag and Seasonality Considerations in Evaluating AVHRR NDVI Response to Precipitation[END_REF] [START_REF] Ji | Lag and Seasonality Considerations in Evaluating AVHRR NDVI Response to Precipitation[END_REF], working on biweekly series, present in their equation [START_REF] Rispe | Forecasting the incidence of aphids using weather data[END_REF] an expression which would be the strict discrete equivalent of the equation (1) of the present article, excepted that the dependent variable is estimated in t while in our equation (1) Y is not indexed by the time but supposed to be observed only once during a period of time (typically the year). They refer also to the lag distribution model of Seddighi (2000) [START_REF] Seddighi | Introductory Econometrics: A Practical Approach[END_REF] in a a book on econometrics comprising a comprehensive chapter on lag analysis. There again, the model of equation ( 1) is the continuous version of the one defined by this author p 118, equation 4.28., with the same noticeable difference. In the materials and methods section of this article, we will assess more precisely the correspondence of this model with the "continuous time lag distributed models" approach defined by Sims [START_REF] Sims | Discrete approximations to continuous time distributed lags in econometrics[END_REF] and present some results about this continuous version. A domain where lag distributed models have shown their usefulness is that of plant phenology,especially when recorded by remote sensing. [START_REF] Wang | Spatial patterns of NDVI in response to precipitation and temperature in the central Great Plains[END_REF] [START_REF] Wang | Spatial patterns of NDVI in response to precipitation and temperature in the central Great Plains[END_REF] studied so the evolution of the Normalized Difference Vegetation Index (NDVI),in relation to climatic factors, namely precipitations and temperature with long time lags (up to 2 months). [START_REF] Ji | Lag and Seasonality Considerations in Evaluating AVHRR NDVI Response to Precipitation[END_REF] [START_REF] Ji | Lag and Seasonality Considerations in Evaluating AVHRR NDVI Response to Precipitation[END_REF] built a similar model for precipitations, referring to a lag model based upon a biweekly cutting of the rainfall data. Their Figure 1 is conceptually the same as our Figure 1 shown further.

Recently, [START_REF] Hufkens | An integrated phenology modelling framework in R[END_REF] [START_REF] Hufkens | An integrated phenology modelling framework in R[END_REF] incorporated such models in a comprehensive phenology modelling framework, available in the repository CRAN (R development team 2019) [START_REF] Core | R: A language and environment for statistical computing[END_REF], emphasising in their article on the importance of plant phenology for ecological processes. Keenan et al [START_REF] Keenan | On quantifying the apparent temperature sensitivity of plant phenology[END_REF], study the temperature sensitivity of plant phenology. They define it quite simply, following various authors, as the ratio of the variation of the date of occurence of a phenological stage, by the variation of temperature in a given year. They draw from their work interesting conclusions on the effect of the observation time (denominator of the ratio) on the measure of this sensitivity. Incidentally, studying the correlation of the phenological date with the duration of the observation time, they draw exactly the same correlogram as we do in their Figure 1c. In 2013, Chuine [START_REF] Chuine | Plant Development Models. Phenology: An Integrative Environmental Science[END_REF] provided a comprehensive overview of phenology models that enlights the importance of lag distributed models.

Aims of the article and examples provided

In this paper, we analyse some properties of the temporal bivariate correlogram resulting from the systematic scanning of many periods of time, and propose an underlying model in two simple cases (rectangular and bell-shaped window), as well as a technique of data analysis for detecting critical periods. We sketch also the properties of some other filtering windows, the Beta shape and the polynomial one. This technique is then applied, to various ecological problems, namely the outbreaks and biology of aphids and, to join the important field of plant phenology, the occurence of a given development stage of of the ash tree (Fraxinus excelsior L.) found in a public database.

Material and Methods

Modelling critical periods 1.The regressive process.

The idea of critical period can be expressed mathematically by a regressive process between a time series U and a random variable Y : The transfer function between  

Ut and Y

is the a scalar product (product of functions in a Sobolev space) between

 

Xt and a window

function   t  .
Y , the exogeneous variable, is therefore the sum of this transfer function and an error  .

  = ( ) T Y k t dU t     ( 1
)
where k is a regression coefficient,  an random variable associated to each realisation of the process,   .  a window function weighting the time series () Ut. When needed, Y , () dU t and  will be indexed by i to infer the model from replicates. With the rectangular window (a and b) the limits of the critical period is sharp, with a bell-shaped window, (c and d) they are fuzzy.

The model is illustrated by Figure 1, with two kinds of weighting functionss (.)

 . For statistical analysis, it is convenient to consider   Ut as a pathwise solution of some stochastic differential equation such as

    = ( ) dU t X t dt dW t   (2)
where

 

Xt is the derivative of a trend or drift, and

 

Wt the standard brownian motion amplified by a standard deviation or volatility  .

For example, the daily record of mean temperatures in a year can be considered as the sum of a seasonal sinusoidal function and of a random noise. This noise is more realistically modeled as a "pink" noise than by a white noise, taking in consideration some autocorrelation, but at the moment we do not enter into such accurate considerations.

By substituting (1) in (2), the model splits easily in the sum of a deterministic and of two stochastic parts: [START_REF] Cochrane | Continuous-Time Linear Models[END_REF], except that we focuse on the "far" past and not at the close past of the process.

  = ( ) ( ) ( ) TT Y k t X t
Furthermore, the Y response is completely dissociated temporally from the time series from which it originates. Actually, it is enough to consider that Y is observed at any time .

obs tT 

Y is unique, and only linked to a realisation of the series, what we shall explain later. We focus our interest on two simple cases of windows defined on R : i) a rectangular window

  ( ) = ( ) ( ) t H t H t        (4)
where

 

H t u  is the Heaviside unit step function switching from 0 to 1 on u. For proper scaling, we will set

1 =   
to ensure that the weights sum to one over R . We assume also that >> T  . This last assumption models the fact that Y is observed after that the effect of   We note that these two shapes of windows are not at all the only possible, an infinitely many others can be imagined. With those two shapes of kernel, some useful results can be mathematically obtained, which are fully developed in appendix A where two other types of windows are briefly discussed. In the first case, and when replicates are available the process reduces obviously to the simple form:

    = i i i i Y k X t dt k dW t      (7)
Where   i dW t denotes the pathwise i th realisation of the brownian motion. To distinguish clearly the two sorts of error, we can set:

  = ii A k dW t     (8)
as representing the process error due to the random noise in the case of a simple rectangular window. The properties of the brownian motion ensures that   =0 EA and

    22 = var Ak     (9)
The pathwise independence of the brownian motion ensures also the pathwise independence of the errors i A . Equation ( 7) is therefore an ordinary regression model, although implying both a stochastic integral, and an observation error i  on which we can make any kind of hypothesis, conditioning only the type of regression we can carry out (normal, Poisson, binomial, beta, etc.)

If the replicates consist in a set of years (or of any period of interest), that later process may be illustrated as in Figure 1b, in the case of a rectangular window: Each year a given variable is sampled one time, and is regressed on the past. The past is taken into account through the window

  t 
, and is therefore proportionnal to the area between  

Xt and the time axis in an interval   , . In the second case, the integral is weighted by all the past, but its contributions are concentrated around the middle time of influence * t . The boundaries of the influential period are fuzzy.

The major point of interest in such processes is to determine wether the shape and boundaries of the window could be estimated. For that purpose, we shall devise an empirical correlogram, and search its properties in the two cases described above, (rectangular and bell-shaped function). We hypothesize that the Bravais-Pearson correlation coefficient can be used as a criterion function to evaluate the window of influence of   Ut. This point is proven in appendix A1.

Exploring other shapes of the window

An infinity of functions (.)

 can be used in equation( 1) to filter the time series () Ut. Basically we use the simpler, the rectangular window, to find critical periods. After that a period is found, its bounds can be used as initial boundary guess to find a different window. Seddighi [START_REF] Seddighi | Introductory Econometrics: A Practical Approach[END_REF] and Sims [START_REF] Sims | Discrete approximations to continuous time distributed lags in econometrics[END_REF] insist on the general difficulty to find a good shape for the equivalent function in lag distribution models, without setting some strong constrains on them. In the discrete case, trying to find an unconstrained set of coefficients leads quickly to an overparametrisation of the model assorted with a weak gradient, making the convergence of most algorithms unsure and unstable. Some types of functions are classical: the exponential decreasing from present to past, and polynomial structure. We did not keep the decreasing exponential because, as we set it in § 1.1, our model is not fully of type "lag distribution", the dependent variable having no temporal link with the predictor. We considered three shapes (plus the rectangular one) either for their ecological relevance or for their versatility:

1. Rectangular (the default window).

2. Gaussian (bell-shaped, no precise limits).

3. Generalized Beta (strictly bounded, very versatile).

4. Polynomial (different shapes, bounded and normalised by its integral for summing to 1 in the boundaries of the window). This window can thus be strictly bounded and is also very versatile in shape. We chose to limit the degree of the polynomial to 5.

In appendix A, the rectangular and gaussian cases are discussed in depth. The beta and polynomial cases are just outlined as providing few analytic results but entering in the general frame of the equation (A23) in Appendix A. In the package Rcriticor, three novel functions were added to optimize the choices 2,3 and 4. We'll give the results of their use in the "Results" part.

Precision of estimations, tests

Testing the significance of the best correlation coefficient as found in § 2 is not straightforward as all the calculated coefficients are themselves correlated. The task seems mathematically untractable and that led us to adopt some simulation approaches. We chose two of them: Monte Carlo permutations for testing the significance of the higher peak and a bootstrap estimation for assessing significance limits. Those two strategies are used in both implementations described lower in § 1.3.

Monte Carlo permutation test for peaks significance

The choice of a maximum correlation coefficient in absolute value leads to a high risk of spurious correlation. The classical test of the correlation coefficient, using the t approximation

2 = 2 t r nt  (10)
is flawed in this case. Its significance value is however given in both softwares as an indication.

Alternatively, a p-value is calculated by random permutations of the dependent values with replacement. Given a threshold risk  , the observed peak (or sink, respectively) is considered as significant if it occurs in les then 100 percent of the permutations.

This procedure allows also to draw contour maps of the areas including significant correlations in the plane (a,d) defined above ( = d b a  ).

An example is given in Figure ??c showing the area where the correlation between the sums of temperatures above 3 o C and the peak of outbreak for the grain aphid Sitobion avenae F. are significant. The observed peak, figured as a black dot, is inside a significant area at risk = 0.05  . Criticor and Rcriticor give also an histogram of the expected distribution of the correlation coefficients under the null hypothesis that there is no correlation between the independent series and the dependent set. Examples are given in Figure 2.

Bootstrap tests for the confidence limits of critical periods and peaks

Another important question is to assess confidence areas for the best correlation coefficient and for the asociated bounds of the corresponding critical period. In both softwares, a bootstrap procedure is proposed to achieve this task. At each bootstrap run, a sample of the available cases of the dependent variable is drawn with replacement and each case is conveniently associated to its corresponding time series (typically a year). Each run produces a pseudovalue and a given number of runs (1000 or more) permits an estimation of the variance-covariance matrix of the maximum (respectively minimum) correlation coefficient, and of a and d , respectively the beginning and duration of the associated critical period. This is done by computing the Efron's bootstrap estimator [START_REF] Efron | Bootstrap Methods: Another Look at the Jackknife[END_REF]. Let us note that another estimator of the variance-covariance matrix of those coefficients may be obtained via the Fisher's information matrix (see appendix B).

Computer implementation

Two implementations of the method were achieved: one, maintained by the french National Research Institute for Agriculture food and Environment (previously INRA, now INRAE) under the acronym Criticor, and the other maintained by J.S. Pierre is available in CRAN as a package and named Rcriticor. Criticor is written in C and Rcriticor in the R programming language [START_REF] Core | R: A language and environment for statistical computing[END_REF]. As all time series available are discrete, the softwares implement discrete versions of the process described above, where integrals are replaced by discrete sums of trapezes. Some example of results are given here below.

Set of data analysed

Grain aphids peaks of population

Sitobion avenae F., the grain aphid,is a direct pest for cereal crops, especially for winter wheat.

Forcasting models were achieved by different method by our laboratory since 1975, and population in the field were recorded during 15 years among which 7 were intensive. The outbreaks of Sitobion avenae F. were firstly correlated to the sum of day degrees above 3 o C [START_REF] Pierre | Un modèle de prévision des pullulations du puceron Sitobion avenae F. sur blé d'hiver[END_REF] in the month of february. This was obtained by linear multiple regression in which the independent variables were the sum of temperatures in each month preceding the outbreak. This statistical procedure was then followed by the research of critical periods [START_REF] Pierre | Une methode statistique et graphique de recherche des periodes de l'annee ou les populations animales sont particulierement sensibles a une composante donnee du climat (periodes critiques). Application au cas des pucerons des cereales[END_REF], the example given here.

Aphids were sampled weekly, and the variable studied here is the recorded peak of population. This example is provided, along with the data, in the CRAN package Rcriticor. Let us note that we switched after that toward mechanistic models of population dynamics [START_REF] Plantegenest | Assessment of the relative impact of different natural enemies on population dynamics of the grain aphid Sitobion avenae in the field[END_REF] themselves incorporated in the commercial package Colibri® belonging now to the firm Bayer Cropscience.

Induction of male morphs of the cherry-oat tree aphid

This example comes from our common work with Rispe et al [START_REF] Rispe | Effect of climate on the proportion of males in the autumn flight of the aphid Rhopalosiphum padi L. (Hom., Aphididae)[END_REF]) trying to find the factors influencing the production of sexuate morphs in other cereal aphid, Rhopalosiphum Padi L, the bird cherry-oat aphid. This species spends most of its time as parthenogenetic morphs, but autumnal conditions induce the production of sexuate morphs. Among them the males are the easiest to record, as they are currently caught in the English and European suction traps network.

The photoperiod is well known as being the main factor influencing this cycle switch

Occurence of the stage 11 (first leaf unroll) in the ash tree

This series was taken from the public database PEP725, devoted to phenological data. In this database, we chose the ash tree in the village of Cardedeu (Catalunia, Spain) because it provided a long series of observations spanning from 1953 to 2000, that is 47 years. Unfortunately somme yearly observations lack and, ass a result, only 23 years are available: 1953,1955,1969,1973,[1980][1981][1982][1983][1984][1987][1988][1989][1990][1991][1992][1993][1994][1995][1996][1997][1998][1999][2000]. We chose as phenological data the date of leafs unfolding, a good marker of spring revival of the trees. The climatic data were obtained from the meteorologic service of Catalunya which give access freely to the data. Unfortunately, for the site Cardedeu, only monthly means were available and only for temperature and rainfall. To recover one part of the climatic variablity at a scale smaller then the month, we achieved a daily interpolation of these data by cubic splines before applying the research of critical periods.

Obviously, this procedure leads to work on a smoothed version of the climatic series and the stochasticity of climate is known less accurately than for the other examples. The R scripts of reconstitution of the data are provided as supplementary material.

Results

The sensitivity to late winter temperatures in the grain aphid populations

Aphids outbreaks occur generally in june, and thus february appeared as a critical period for temperature, preceding the event of interest (outbreak or population peak) of five to six months.

The present method was then used to define more accurately this period. Figure 3a shows the correlogram obtained in this case. The figure is much more noisy than the theoretical figures 6 and Figure 7 (in appendix A), but looks more like Figure 7, evoking a bell-shaped influence. The best correlation was found as = 0.981  for a period beginning at day 37 after january the first, i.e. february 6, and lasting 44 days, that is until march 22. This new period, longer than the month, was used further for forcasting [START_REF] Pierre | Une methode statistique et graphique de recherche des periodes de l'annee ou les populations animales sont particulierement sensibles a une composante donnee du climat (periodes critiques). Application au cas des pucerons des cereales[END_REF]. Figure 3b shows the confidence area of the maximum correlation position by bootstrap, while Figure 3c (significance map) and Figure 3d (histogram of permutations) ensure that the maximum correlation observed is significant at the risk = 0.05  Other shapes of windows.

On this example, with a gaussian window we get, at most, = 0.788  , a lower value than with the rectangular window. Trying the Beta function, we find = 0.981  , exactly like with the rectangular window, resulting from the fact that the optimization algorithm converges toward =1 

and =1

 , the values where the Beta distribution is identical to the uniform one. The best correlations found with a polynomial window were obtained at .9709 for the second degree and .9755 for the fifth degree. To summarize, there is no reason here to challenge the rectangular window.

The sensitivity to induction factors of sexuality in the bird cherry-oat tree aphid

The Criticor method shows that, surprisingly, high temperatures in late summer and autumn favor the proportion of males in the suction traps. The critical period for temperature, the span of which is estimated from August 24th and October 17th. Figure 4a shows the correlogram and Figure 4b the scatter plot corresponding to the estimated critical period. This shows that, if photoperiod drives the induction of sexuals, the importance of their production is strongly determined by autumn temperature. The observation of a positive correlation is the inverse of our initial hypothesis. This fact has a strong importance in aphids population dynamics: sexuals go back to the winter host, the cherry-oat tree where females lay diapausing winter eggs which are cold resistant.

The part of the population which remains parthenogenetic, at the opposite, stays on winter and volunteer cereals. Those latter are susceptible to be killed by harsh frosts. The knowledge of this fact has applications on risk forecasting: this species transmits viruses on winter cereals in autumn, and that as more as the population is constituted of more parthenogenetic morphs and less sexual morphs.

Other shapes of window.

All our attempts to fit other shapes of window than rectangular resulted in failures. In all cases, the absolute values of the correlation coefficients were lower than in the case of rectangular window.

In some cases (Beta, polynomials of degrees 2,3,4,5) the correlation coefficient turned to be negative. We do not present these results in details.

Ash tree phenology in Spain

The date of occurence of the stage 11 (leaf unfolding) of the ash tree was correlated with past temperature and rainfall.

Temperature

We found a negative relation (Figure 5a) between temperature (threshold 0 0 C) with a correlation coefficient of -0.667 (maximum in absolute value) for a period beginning at Julian day 43 (February 10) and during 74 days (till Julian day 117 i.e. April 26). The bootstrap estimation (1000 resamplings, Figure 5b) corrects somewhat this crude estimation, obtaining a slightly better correlation of = 0704   for a slightly shorter period beginning at day 55 (55.196 estimated) and lasting 59 days (exact estimate : 59.147).That is a period of two months between February 22 and April 23. Figure 5d shows the scatter plot obtained for this last period, a fairly linear relationship, rather surprising as, theoretically, following the theory of day-degrees, an inverse hyperbolic relation is expected. Figure 5c shows the significance map obtained, showing a very large area in which the influence of temperature is highly significant (black area, p<0.01). This result is obviously banal but shows that day-degrees have little influence on the leaf unfolding of the ash tree before late february. This is probably because before that, winter coldness is necessary for the breaking of endodormancy. The effect of day-degrees accumulation acts therefore only after that dormancy breaking.

Rainfall

Rainfall has a longer influence than temperature (Figure 5e), the period detected covering three months, from Julian day 61 (March 2) to Julian day 184 (July 4). The absolute value of the correlation coefficient is lower than for temperature, ( = .465  ). On the correlogram (Figure 5a) it appears as the merging of a short period beginning at mid March and lasting a few days and another much later not visible on the figure. The bootstrap then corrects strongly this estimation by setting a new one for a period ongoing from day 65 (bootstrap: 65.456 -March 4) to day 139 (bootstrap: 139.095 -May 18, Figure 5g). This estimation has the interest of being more predictive although with a lower correlation coefficient ( = 0.334  ). In effect, this period ends just at the earliest beginning of leaf unfolding. We therefore kept this estimation for further multivariate regression. The scatterplot show a positive relation but with a noticeable deviation from linearity (smoothing spline, red curve,Figure 5h). This encouraged us to use a General Additive Model (GAM) to combine temperature and rainfall in an explanative model. The positive effect of rainfall suggests that it delays the leaf unfolding of the ash tree. This is a little surprising as rain is supposed to be beneficial for the plants development. The inverse correlation of rainfall and temperature may be suspected here: the rainfall results often in a loss of temperature.

Bivariate model

To combine the effects of temperature and rainfall we used a General Additive Model. The choice of GAM was driven by a detectable nonlinearity of the residuals (Figure 5h). We compared two families of errors, the classic normal one and the Cox Proportional Hazard model as the variable of interest is of the survival time type.

Comparison with a monthly multivariate approach

An alternative frequently used is to introduce in a multiple regression a series of predictors obtained by cutting the time series into equal segments, on which the climatic variable is summed or averaged. Those segments can be months, half months, weeks, etc. The inconvenience is to introduce an artificial splitting of the climatic series of interest, and to manage a large set of predictors leading to the difficulty of finding a "good" subset of predictors. To illustrate this point, we present the results obtained in the ash tree case. The stepwise procedure was based on the Bayesian Information Criterion (BIC), more penalizing than AIC and leading to a slightly more parsimonious model (only three predictors, instead of four with AIC). We can then compare the characteristics of the two approaches. The comparison (Table 3) shows firstly that the correlogram based approach is basically more parsimonious than that by multiple regression, although, in the present case, the stepwise procedure does the job well. It shows also that there is a gain of accuracy in the definition of the periods, because the correlogram based approach allows to find bounds inside one month, at any place. In our example, we can see that the influence of rainfall in March,

Searching other shapes of the window.

We give here shortly the results of our trials to fit other shapes then rectangular for the filtering window. We tried, as announced in Materials and Methods, three cases besides the rectangular one: gaussian (without precise limits), Generalized Beta (strictly bounded) and polynomial (bounded and normalised). Table 4 summarizes the results obtained. As in the preceding cases, the results are somewhat deceiving. The gaussian windows does not improve the fit, the Beta case challenges very few the rectangular one (only giving less importance to the upper bound of the window) regarding the temperature, and converges to the rectangular one (uniform) regarding the rainfall. The polynomial shape gives the best results for degrees 3 and 5 for temperature and for degrees 3 and 4 for rainfall. The correlation coefficients get higher results than that of the rectangular shape only in the two polynomial degrees for rainfall. In the case of temperature, degrees 5 and 3 give contradictory results, the degree 5 leading to a positive correlation, and an artefactual scatterplot: one of the value gets the huge value of 1 billion, while all other points are flattened around zero. It is clear that the simulated annealing procedure found a spurious maximum which must be discarded. Concerning rainfall, the polynomial approaches are more interesting as giving more importance to the end of the period than to its beginning. At degree

Discussion

We provide here a method for exploring the specific question of cumulative and delayed effects of a time series onto a measurable event. This case is relatively frequent in biology, ecology and psychology. We show here that our method is efficient when such effects are suspected. This article has no other purpose than to contribute to statistical methods for exploring the past effects of some factors onto measurable events in biology. We bring some mathematical results about its statistical background. Under two sorts of transfer functions the maximum exists, and is theoretically obtainable in a quasi closed form (see appendix B for details). The maximum of correlation coincides with the maximum likelihood estimate (proven in appendix B). The visual inspection of the correlogram is useful to detect eventually a local maximum and/or a combination of separate critical periods. A R package is associated to achieve this exploration. We restricted our investigations to some simple cases. Some refinements are provided later: how to manage non-stationarity in the independent series such as regular trend, seasonality, truncation.

Multiplicative trends or other functional forms would require the Ito or Stratonovitch calculus which were not required in our simple models where the white noise can be considered separately and additively.We provide here a method for exploring the specific question of cumulative and delayed effects of a time series onto a measurable event. The modelling of critical periods that we propose here is purely phenomenological. It says actually nothing about the underlying mechanisms which determine these influences of the past. It has however the advantage to focus the study of these eventual mechanisms in precise periods of time. This can be of great interest for building relevant mechanistic models in the first step of an ecological research. We show the kinship with the lag distribution models, emphasising however on a difference which results in a simplification: the lag itself is not modelled, neither estimated. The interest is focused on the size of the response. Even when the response is the date of realisation of an event such as the occurence of a phenological stage, this date is considered as a quantitative response, not as a time lag between a cause and an effect. Although it can be useful for building models in which the effects are delayed, it is not conceived as an "all purpose" toolbox susceptible to solve all type of related problems.

All the examples given here are purely local. As more and more data become available along with the corresponding meteorological records, it would be interesting to consider multisite strategies to study critical periods. A first idea, the simplest, is to use our method (or another) separately in each place, obtain local predictors, and build a multisite variance-covariance model where the sites play the role of a random factor. This could be valuable for distant sites few correlated together. Other more sophisticated ideas can also be considered but risk to ruin the simplicity and exploratory character of the method. For multisite phenological data, the approach of Hufkens et al. seems more recommandable for such a purpose.

The method leads to parsimonious models of statistical explanation and prediction. This advantage is shared by all the published methods aiming at finding delayed effects in ecology and other sciences. It adds a tool for such investigations. It can take place in a toolbox aside with the proposals of Ogle et al.et al. [START_REF] Ogle | Quantifying ecological memory in plant and ecosystem processes[END_REF] and of Hufkens et al. [START_REF] Hufkens | An integrated phenology modelling framework in R[END_REF]. The system proposed by Ogleet al. [START_REF] Ogle | Quantifying ecological memory in plant and ecosystem processes[END_REF] is a comprehensive modelling framework for the explanation of an ecological phenomenon subject to several influences of the climate components. It requires a strong corpus of hypothesis about the event studied. When the ecological subject is relatively well known, the work of these authors can lead to a robust modelling of the system. This is generally not the case in animal population dynamics where statistical relationships between populations and climate are rarely explained by mechanical models, for diverse reasons. In the case of cereal aphids, for instance, a statistical relation emerges from the late winter temperature and the outbreaks in June [15][16]. A mechanistic model of population dynamics proved to be useful when an sampling of population is available in March-April (Plantegenest et al. 2001[18]). Both models have could never could be connected simply because cereal aphids in winter are so rare that no accurate sampling of their density is available. Such difficulties are frequent in the case of insects, which cross periods of very low density in the year. Ogle et al. [START_REF] Ogle | Quantifying ecological memory in plant and ecosystem processes[END_REF] insist also on the concept of ecological memory, a concept highly relevant to the case of critical periods as all goes as if the accumulation of the climatic series is retained during the critical period and forgotten after that.

The question arises also to extend the method to a direct multivariate approach by combining all meteorological data together. This can be considered by two ways: Firstly by eigenvector techniques such as PCA and try to maximize the effect of a linear combination of several weather records on the biological observations, secondly by searching in a 2 p space ( p being the number of meteorological series) the most influential set of periods caracterised by their boundaries. The inconvenience would be the difficulty of visualising such a space of responses.

Nevertheless, we are currently trying to achieve both of these approaches. The tool presented here can find its place in the armoury of statistical methods which allow to study the relationships between climate and populations. With few modifications it could also serve in other scientific fields where delayed effects are suspected whatever could be the time and space scale where they occur.

Deposit for the data

The data corresponding to Sitobion avenae L. are deposited in CRAN, along with the package Rcriticor. Those corresponding to Ulex europaeus L. and to Rhopalosiphum padi L. are deposited in Mendeley. The ash tree data are available in the PEP725 phenologic database (http://www.PEP725.eu), and the corresponding climatic data on the site of the .eu) Meteorological Service of Catalunya (https://es.meteocat.gencat.cat/?lang=es). the R scripts necessary to reproduce our calculations of these data are provided as supplementary material.

5 Aknowledgements Y . If such a period exists, then a maximum or minimum should be observed on a map plotting (.,.)

 against a and = d b a  . Such empirical map can be seen in Figure 2a.

We eliminate the trend or drift considered as fixed and keep only the stochastic parts. Then,  

  = ii Y k dW t      (A4)
The covariance implies to get the expectation of their scalar product, the value of which is proposed below:

Proposition:   2 , .= ab E S Y k (A5)
Where  is the length of the intersection between intervals [ , ] ab and [ , ] .

proof:

Let us consider a single replication:

        2 ,, .= bb a b i i i i i i aa E S Y k dW s dW u dW t          (A6)
and for any replication

          22 , .= bb a b i i i aa E S Y k E dW s dW u E dW t E                    2 = k (A7)
As the increments of

  i Wt are independent,     0 ii E dW t dW u    (A8)
if and only if = tu . in this case:

      2 == i i i E dW t dW u E dW t dt         (A9)
This condition is verified only in the common part of the intervals   . From the habitual properties of correlation coefficients it results also that a and b are least square estimates of  and  . Figure 5 shows , in perspective and contour plot the The shape of the surface (Figure 6) is straightforward from formula (A11). All is governed by the relative sizes of d and 

. The interval [ , ] ab , the length of which is d , the abscissa of figure 5, grows monotonously from left to right. Two lines, the horizontal = a  and the negative bissector = ad   share the plane in four quadrants. In the lower left quadrant, a is lower then  and  is null until da  . Then,  grows faster then d . In the lower right quadrant, > ad   and thus, =     , a fixed value, and the ratio decreases from left to right with d .

In the upper left quadrant, > A.2 The case of a bell-shaped window.

In case, the window has no definite bounds. The interesting things are to estimate * t and c the pseudo standard deviation of the gaussian window. Let us return to expression (A6). It takes now the form: And, for the same reasons of independence as in the previous case, and because

             
  =0 i E  (A15)
this reduces in:

        2 22 , .= b ab a E S Y k E t dW s      (A16) As   t 
is fixed for all t , we can factor the expectation in:

        2
      , , , 2 , 1 == 22 ab
F h c F h c F h c hh     (A26)
the derivative of which is:

    3 2 2 , 1 2 ' = , 2 4 h F h c f h c h h    (A27)
Equating to zero we find:

    2 , 1 ,= 4
F h c f h c h  (

A.3 Generalized Beta case

The Generalized Beta density distribution can be used as weighting function, just as the gaussian distribution. It has one advantage: that of being precisely restricted between two boundaries. Its other feature is an extreme versatility, with a shape varying from exactly rectangular window to humped one, symmetrically or not, L or J shape, and even in U shape depending on the values of its two parameters,  and  . The scheme of reasoning is essentially the same as for the gaussian case, as we can start at equation (A23), just replacing the function F which is the integral of the gaussian function by the integral of the generalized Beta distribution. This leads to :

  22 ( , , , ) ( , , ) ( , , ) = () ab Y k B b B a d          (A30)
where (., , )

B  designs the integral of the Generalised Beta density of parameters  and  and of support [ , ] ab but, by definition, this difference is equal to 1 and thus:

22 ( , , , ) = ab Y k d     (A31)
from this formula, the correlation coefficient seems to depend only on the length of the intersection of the trial window and of the "true" range of action of the Beta function but, this optimal range depends itself on the coefficients of the Beta distribution. There is thus no useful analytical solution but the range and the coefficients  and  of the Beta distribution can be found by an optimization procedure. We used for that the function optim of R , with as method either Nelder-Mead (the default), or SANN (simulated annealing) or L-BFGS-B. This last method allows to specify restrictions (box constraints) on the parameters and is necessary to maintain the two shape parameters of the Beta distribution greater than zero or sometimes greater or equal to 1 to obtain a humped shape.

A.4 The polynomial case

Let us choose for (.)  a polynomial of degree p scaled to unity on

( ) = ( ) ( ) ( ) 23 b a P s ds a b a b a b a       (A35) 2 2 3 3 0 () ( ) = 11 ( ) ( ) ( ) 23 Pt t a b a b a b a       (A36)
or, putting

2 2 3 3 0 11 = ( ) ( ) ( ) 23
I a b a b a b a      and = i i a A I : 2 0 1 2 ( ) = t A At A t    (A37)
and so on for higher degrees. Reasonably, to avoid overparametrisation and keep a simple shape for the window, we limited our investigations to the 5th degree. Under this form (A36), the function () t  sums to 1 in the interval [ , ] ab and is the normalised form of the polynomial. So as that for the Beta shape, there is no interesting analytical solution and the same procedure of optimisation was used.

Appendix B: Maximum correlation and maximum likelihood B.1: Identity of both (normal errors case)

Let us go back to formula (3) adapting it to a rectangular window and supposing an intercept  :

    = i i i i Y k X t dt k dW t          (B1)
 being a random variable that we'll suppose distributed distributed normally as   2 0, r N  with the standard hypothesis of independence of variance and no autocorrelation. This can be rewrote as: Which is realized in expectation. This has as consequences that the estimates of  and  are slightly biased and asymptotically unbiassed for large n s. The important point is that the determination of  and  by maximum likelihood coincides with the maximisation of the correlation coefficient. After that, it is straightforward to find  and k from equation (B5).

= i i i Yk     (

B.2: Consequences

From k we get  the maximum correlation coefficient from the common formula:

  var = () var k Y   (B16)
And we have estimates of the variances of ˆ,,    , and k by calculating the Fisher's information matrix from the log-likelihood value in the neighbourhood of these parameters. The difficulty comes from that the partial differentials of the log-likelihood with respect to  and  include the realization of a stochastic process ( ( , )    ) and thus are not the derivatives of continuous functions. We chose to treat the differentials in mean, one way to stabilise the Fisher information matrix. These estimates can be compared with those obtained by bootstrapping. Those are in development for the next version of Rcriticor.

Highlights

 The study of relationships between climate and animal or plants populations is of prominent importance especially in the context of climate change.

 We focuse here explicitely on the hypothesis that one series may have an influence on a further event by a bounded window of time, long time before the effect. This window is called "critical period".

 We propose a bivariate correlation function to detect such periods. Some close form estimators are proposed for it, bootstrapping and Monte-Carlo methods are proposed to estimate its variability, a CRAN package is available for that.

 We give three examples of ecological applications of the method.

Figure 1 :

 1 Figure 1: Illustration of the transfer function and shape of the window function. a: rectangular window function b: resulting transfer function c: bell-shaped window function d: resulting transfer function The red bar and red dot indicate the result of the process, and the red arrows the application from the integral to the final observed value of Y .With the rectangular window (a and b) the limits of the critical period is sharp, with a bell-shaped window, (c and d) they are fuzzy.

Fig 1

 1 Fig 1 illustrates its shapes in both cases.

Figure 2 :

 2 Figure 2: Histogram of the absolute values of the maximum correlation coefficients obtained through 1000 random permutations of the dependent variable. a: in the software CRITICOR (INRAE).The observed coefficient (-0.85) is thus declared significant at risk = .05 

Figure 3 :

 3 Figure 3: Intensity of outbreaks of the grain aphid Sitobion avenae F. a: Correlogram crossing the intensity of outbreaks of the grain aphid Sitobion avenae F. and the sum of temperatures above 3 o C during 7 years in Britanny (France). The observed maximum corresponds to a period beginning at day 37 (february 6th) and spending 44 days (until march 22th). b: The same correlogram with bootstraping (500 subsamplings). The red diamond indicates the bootstrap estimator of the maximum of correlation. Red dots figure individual pseudovalues. The

Figure 4 :

 4 Figure 4: a -Correlogram obtained between the proportions of males of R. padi in the autumn flight and mean summer temperatures at Rothamsted. Periods are from July 19th to October 17th, with duration from 1 to 60 days. The arrow indicates the highest correlation peak. b-Linear relationship between the corresponding sum of temperature and the rate of males production.

3 :

 3 April and May are discarded by the process of predictors selection, at the advantage of June alone, while the scanning process detects a critical period covering march to may. This advantage of parsimony is shared by all the methods which research precisely the lags of influence of the Comparison of the structures of the models obtained by Rcriticor (correlogram based research) and by multiple regression.

  defined as in equation (2), § 1.1.1 . This interval will be called further the trial period. ( , ) ab  is the correlation coefficient between Y and , ab S . This correlogram can be used empirically to search for the period which influences most the observed variable

  where  is the scale of   Xt, and d the length of the interval   , ab . To go further, we have to choose a shape for   t  . Let us treat the case of the rectangular window defined in equation (4) §1.1.1. We have to correlate a vector whose components are

  correlation coefficient is proportional to the length of the intersection of the window interval and of the trial interval and inversely proportional to the square root of the length of the  is maximum is an unbiassed estimator of the influence interval[ , ]

  shape of the function   , ad  , with its definite peak on the point   ,  . We chose the representation in the plan   , ad rather than   , ab not to leave the graph empty under the bissector.

  in the upper right quadrant. The whole set results in a pyramidal shape with an acute peak. Of course, this is only a theoretical result, altered by uncontrolled sources of noise in the real world.

Figure 6 :

 6 Figure 6: The shape of the expected correlogram under the hypothesis of a rectangular window: a: perspective view; b: contour plot.

Figure 7 :

 7 Figure 7: The shape of the expected correlogram under the hypothesis of a bell shaped window: a: perspective view; b: contour plot.

  Fc being the normal cumulative distribution function of standard deviation c . The expression is the covariance as both As the Gauss function is symmetric, let us conjecture that the interval ** , ab   , if any, that maximizes , ab  is symmetric, left and right of * t . The mean of   ., Fc is supposed to be * t meaning that   ., Fc must be understood as   is only possible if ( , ) = ( , ) f a c f b c , implying that b is symmetric to a with regards to * t QED Determination of h Equations (A25) however a and b indeterminate. So let us put =2 b a h  and, knowing that a and b are symmetric, let us try to optimize half span h . Without loss of generality, let us put * =0 t in equation (5,?2.2) (a simple change of origin), and considering that , k  and Y  are positive constants, we reduce the problem to maximize:



  involves the result of a stochastic integral, but, as the different replicates of the time series as predictor, must be considered as fixed. Supposing that  , i   is known ( and  are known), the maximum likelihood estimators of  and k are well known: Let us A be the two columns and n and  are unknown, and will be replaced, in what follows, by the trial parameters a and b . The preceding equation, however, means that once a maximum likelihood is found for those two parameters, the solution for  and k is straightforward, being given by equation (B5).

From 2 )

 2 for any other period   , ab such as > ba . We can wright the likelihood of the observed vector Y as follows:We have then to equate to zero the partial derivative of this expression with respect to the bounds a and b . Let us begin by a : firs term equals to 0 by construction of residuals, and the second is null in average, by independence of the increments   i dW a and of the residuals of the regression i  . So, = a  ensures the nullity of the derivative, in expectancy. Let us now derive with respect to b series of calculations as for a , we get:

Table 1

 1 , a and b summarises the results. Although theoretically more correct, the Cox Proportional Hazard model does not do better than the standard normal one.

	The standard regression model is more powerful, although its Akaike Information Criterion is
	higher than that of the Cox model. Furthermore, the residuals pass the Shapiro-Wilks test of
	normality. The lack of significance of the rainfall effect would incite to take only into account the

effect of temperature. We kept it, however, because of the improvement of the AIC it gives. At end, we get a parsimonious model, with only two predictors, directly issued from the analysis of

Table 1 :

 1

	model	components	BIC
	complete	12 predictors	184.76
		mean temperatures: January to June
		rainfall: January to June
	best predictor alone	mean temperature in March	177.26
	stepwise regression	temperatures: February,March; rainfall: June	173.60
	GAM	temperatures:February, March; rainfall: June (smoothed)	171.74

Results for the univariate (temperature alone) and bivariate (temperature and rainfall) models, as predictors for the stage 11 of the ash tree in Cardedeu (Spain). a: GAM model, gaussian family, b: GAM model, Cox Proportional Hazard model.( 2 R : pseudo 2 R of Cox and Snell) Figure 5: Research of predictors for the stage 11 of the ash tree in Cardedeu (Spain) -a to d: temperature -e to h: rainfall.From left to right: correlogram, bootstrap estimation with pseudovalues and confidence ellipse, significance map, scatterplot corresponding to the bootstrap estimator. In the frame h the red curve is a smoothing spline indicating the non linearity of the relation which led to use a GAM model.

Table 2 :

 2 Summary of model selection with monthly predictors

Table 4 :

 4 Fitting different critical windows on the ash tree phenology (stage 11) in Cardedeu

	4, especially, we

  The shape of the surface (Figure7a) looks like a smoothed version of the one obtained at §A.1 and shown Figure??a. It is organised by the same crossed lines as in the rectangular window case. This case is rather more realistic than that described in 2.1 as nobody expects that the influence of a time series on a phenomenon could begin precisely at a given date and cease at another precise date. The interest of the bell shaped window is to show that, with a somwhat simple research procedure, a peak is detectable, bordering a window of influence whose boundaries are functions of the deviation parameter c .

	h	60.479	. In the case of a gaussian bell-shaped window with the deviation parameter c	, the
	procedure detects the span of influence of the time series as an interval ** [ , t h t h ]  , where h	is
	solution of equation (A28).
						A28)
	h or with the change of variable =	s
			c
				  * zs	=	  * * 21 cs 4 Fs 	(A29)
	This equation has a solution in s which can be found numerically. For = 20 c	,	s	3.024	and

We thank warmly Jean-Pierre Masson, emeritus professor of statistics and probabilities at Agrocampus Rennes, for his carefull reading of the mathematical developments of this article. We thank Anne Atlan for allowing us to use her european gorse data, and Nicolas Parisey for his reading of the manuscript. We thank James Bell and Lynda Alderson of the Rothamsted Insect Survey (RIS, UK) for sending us the original data of aphids males flights. We thank also Dr Sarah Perryman e-RA curator at Rothamsted who made the English weather data available for us.

Appendix A: Shape and extrema of the empirical correlogram A.1 The case of a rectangular window.

Let us consider the following correlation function:


Chaima Boureghda

Raouf Boulcina

Vincent Dorcet

Fabienne Berrée

Bertrand Carboni

Abdelmadjid Debache

Facile synthesis of 5-arylidene rhodanine derivatives using Na 2 SO 3 as an eco-friendly catalyst. Access to 2-mercapto-3-aryl-acrylic acids and a benzoxaborole derivative

Knoevenagel reaction 5-Arylidene rhodanine 2-Mercaptoacrylic acids Sodium sulfite Benzoxaborole

A simple, efficient and environment-friendly procedure for the synthesis of 5-arylidene rhodanines derivatives via a Knoevenagel type reaction was developed using rhodamine, a variety of differently substituted aldehydes and Na 2 SO 3 as benign catalyst in ethanol. Selected 5-arylidene rhodanines were subjected to basic hydrolysis to afford 2-mercapto-3-substituted-acrylic acids. The presence of a boronic acid group is well tolerated in such transformations. In the case of 2-formylphenylboronic acid, this sequence opens access to a new benzoxaborole derivative.

Introduction

Heterocycles display an impressively high biological potential and were also used as important intermediates in the area of organic synthesis. Derivatives containing the 2-thioxothiazolidin-4-one (rhodanine) substructure are an interesting and important part of this wide library of bioactive compounds and exhibit a wide variety of biological activities, 1,2 as anti-fungal, 3 antibacterial, [4][5][6][7] anti-inflammatory, [8][9][10] anti-cancer, [11][12][13] anti-Alzheimer, 14,15 anti-viral [16][17][18] and anti-diabetic (Figure 1). [19][20][21][22] More precisely, they act, for example, as inhibitors of aldose reductase, protein tyrosine phosphatase PTP, aggregation of Tau proteins, HIV-1 integrase, hepatitis C virus, protease and - In addition, they are versatile synthetic intermediates since they can undergo many chemical transformations, such as hydrolysis in basic medium to lead to mercaptoacrylic acids, 23 Michael-type addition, 24 Diels Alder reactions to access dihydropyrans 25 and 2-thioxopyrano [2,3-d] [1,3]thiazoles. 26 For these reasons, 5-arylidene rhodamines have been the subject of numerous investigations and a number of procedures have been reported for their preparation. The main method of synthesis is the Knoevenagel reaction between rhodanine and aldehydes in the presence of various catalysts, such as, for example, ammonium acetate, 27 2-hydroxyethylammonium acetate, 28 NaOAc/HOAc, 29 2,2,6,6-tetramethyl piperidine, 30 glycine, 31 1-butyl-3-methyl imidazolium hydroxide, 32 NH 4 Cl/NH 4 OH, 33 piperidine, 34 supported K 2 CO 3 on Al 2 O 3 under microwave irradiation 35 or in the presence of an ionic liquid. 36 As a continuation of our research for an new environment-friendly approaches, [37][38][39] we herein propose a new simple and ecocompatible procedure for the preparation of 5-arylidene rhodanines using Na 2 SO 3 as catalyst in EtOH at reflux under conventional heating (Scheme 1). The hydrolysis of some of them was carried out in a second step to provide the corresponding 2-mercapto-3-substituted-acrylic acids 4.

Results and discussion

The optimal conditions to access 5-arylidene rhodanines 3 were first determined with rhodanine and 4-nitrobenzaldehyde as model compounds using different experimental conditions of temperature, solvent and amount of Na 2 SO 3 as catalyst (Table 1). The modification of the initial procedure (r.t., 20% Na 2 SO 3 , EtOH, Entry 1) by heating the reaction mixture at reflux 2h significantly improved the formation of the 5-arylidene rhodanines (Entry 2). This beneficial effect slightly diminished with a smaller quantity of catalyst (Entries 3 and 4), while an increase is also unfavorable (Entries 5-7). No improvement was found by using H 2 O, EtOH/H 2 O (1/1) or THF/H 2 O (1/1) as solvents (Entries 8-10). Having optimized conditions in hands, we then changed the nature of the aldehyde to evaluate the scope and eventual limitations of this method (Table 2). This procedure is tolerant with a large range of functional groups (halogen, alcohol, nitro, ether, acetal, amino, imidazole). Good yields were also observed in the case of heterocycles as imidazole or 1H-pyrazole-4-carbaldehydes, an important class of pharmacophores. [40][41][42][43][44] The excellent chemioselectivity of this method was illustrated by its compatibility, whatever the position on the aromatic nucleus, with a boronic acid moiety, a substituent which is becoming common into drug discovery processes 45 and is also a source of diversity given the versatility of this function.

3a 80% 3b 82% 3c 66% 3d 84% 3e 82% 3f 84% 3n 70% 3t 78% 3u 70% 3v 80% 3o 85% 3p 80% 3q 70% 3r 80% 3s 55% 3g 70% 3h 68% 3j 65% 3k 60% 3l 75% 3i 56%
a Aryl (Hetaryl)aldehyde (1 mmol), rhodanine (1 mmol) and 20 mol% Na 2 SO 3 in 5 mL of EtOH at reflux, 2h for 3a-3n and 3t-3v and 7h for 3o-3s. Amongst the most important transformations achieved on the thiazolidinone cycle are the alkaline hydrolysis to 2-mercapto-3aryl-acrylic acids, 46 which have interesting biological activities as, for example, in vitro chelation capacity of cadmium, 47 inhibition of metallo--lactamase 48 or calpains 1 and 2. 49 They have been also employed as precursors of amino acids, 50 benzo[b]thiophenes 51 and naptho[1,2-b]thiophenes. 52 Selected previously synthesized 5-arylidene rhodanine derivatives are subjected to basic hydrolysis to afford new 2-mercapto-3-arylacrylic acids in good to excellent yields, the boronic acid function being compatible with these experimental conditions (Table 3).

The (Z) stereochemistry of the double bond was confirmed by X-ray crystallographic analysis of 4h. Furthermore, it is worthy to note that the asymmetric unit effectively contains ten crystallographically independent molecules (Figure 2). Dimers are observed with a pseudo-inversion center due to the presence of intermolecular hydrogen bonds between two COOH groups, the dihedral angle between the planes of the aryl ring and the alkenyl double bond varying from 20 to 57°. Moreover, these dimers are stacked along the [10-1] direction of the unit cell with an average pseudo-translation of 5.39 Å. The particular case of 5-arylidene rhodanine 3t which has a boronic substituent at the ortho position deserves special comment. Contrary to what is expected, the 2-mercapto-3-arylacrylic acid 4t was not obtained. The new product was characterized by 1 H, 13 C, and 11 B NMR and mass spectroscopy. Based on these data, a benzoxaborole structure, whose interest in medicinal chemistry is now widely recognized 54,55 has been proposed. The formation of 5 can be rationalized by an oxa-Michael addition to the double bond of 3t, followed by an intramolecular cyclization to create the boroheterocycle and finally cleavage of the 2-thioxothiazolidin-4-one ring by acid treatment (Scheme 2). Scheme 2: Synthesis of benzoxaborole 5 and proposed mechanism for its formation.

Conclusion

In this work we propose a general procedure for the synthesis of 5-arylidene rhodanines with very good yields from a wide variety of functionalized aldehydes. Na 2 SO 3 was used as benign catalyst due to its use as a food additive, its low cost and its availability. Basic hydrolysis easily afforded the corresponding 2-mercaptoacrylic acids. Of special interest is the obtention of benzoxaborole when the boronic acid group is located at the ortho position of the aromatic substituent.

Fig. 1 :

 1 Fig. 1: Structures of some bioactive arylidene rhodanines.

RCHOScheme 1 :Facile synthesis of 5 -

 15 Scheme 1: Synthesis of 5-arylidene rhodanines 3 under Na 2 SO 3 catalysis.Conversion to 2-mercapto-3-substituted-acrylic acids 4.

Table 2 .

 2 Synthesis of 5-arylidene rhodanines 3. Scope and limitations a

a 1 )

 1 Aryl rhodanine (0.5 mmol), 20 mol% NaOH (4 mL), reflux. 2) HCl 2N (13 mL), r.t.

Fig. 2 :

 2 Fig. 2: Structures of (a) a dimer of 4h showing the intramolecular hydrogen bond (b) the asymmetric unit

Table 1 .

 1 Synthesis of 5-arylidene rhodanines 3. Optimization of the reaction conditions

			O	O				O
		HN S	S	+ O 2 N	H	different conditions	HN S	S
								NO 2
	Entry	Solvent	Temperature Cat (mol%)	Time (h)	Yield %
	1	EtOH		r.t.		20		16	44
	2	EtOH		reflux		20		2	84
	3	EtOH		reflux		5		2	82
	4	EtOH		reflux		10		2	75
	5	EtOH		reflux		30		2	67
	6	EtOH		reflux		40		2	68
	7	EtOH		reflux		50		2	68
	8	H 2 O		reflux		20		2	70
	9	EtOH/H 2 O (1/1)	reflux		20		2	73
	10	THF/H 2 O (1/1)	reflux		20		2	67

Table 3 .

 3 Synthesis of 2-mercapto-3-substituted-acrylic acids 4 a

	4x

Acknowledgments

C. B. thanks the Ministère de l'Enseignement Supérieur et de la Recherche (Algeria) for financial support (PNE program).

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org ……..