Patrick Coquillard
email: patrick.coquillard@unice.fr

Alexandre Muzy
email: a.muzy@univ-corse.fr.

Francine Diener
email: francine.diener@unice.fr

Optimal phenotypic plasticity in a stochastic environment

Keywords: Adaptive phenotypic plasticity, plasticity cost, population's genetic diversity, MaxEnt algorithm, exergy

This paper addresses the question of optimal phenotypic plasticity as a response to environmental fluctuations while optimizing the cost/benefit ratio, where the cost is energetic expense of plasticity, and benefit is fitness. The dispersion matrix Σ of the genes' response (H = ln|Σ|) is used: (i) in a numerical model as a metric of the phenotypic variance reduction in the course of fitness optimization, then (ii) in an analytical model, in order to optimize parameters under the constraint of limited energy availability. Results lead to speculate that such optimized organisms should maximize their exergy and thus the direct/indirect work they exert on the habitat. It is shown that the optimal cost/benefit ratio belongs to an interval in which differences between individuals should not substantially modify their fitness. Consequently, even in the case of an ideal population, close to the optimal plasticity, a certain level of genetic diversity should be long conserved, and a part, still to be determined, of intrapopulations genetic diversity probably stem from environment fluctuations. Species confronted to monotonous factors should be less plastic than vicariant species experiencing heterogeneous environments. Analogies with the MaxEnt algorithm of E.T. Jaynes (1957) are discussed, leading to the conjecture that this method may be applied even in case of multivariate but non multinormal distributions of the responses.

Introduction

Phenotypic plasticity is a general and common feature, probably shared by most organisms. This concept refers to the ability of genetically identical organisms to change their phenotype in response to environmental changes in space and time. These changes are favoured because they reduce the fitness variance from generation to generation. As a matter of fact, the phenotypic plasticity enables individuals of a population to colonize various ecological systems, to extend the geographical area of the species and thus to reduce its probability of extinction. For instance, a plant grown in a sunny habitat will exhibit broader leaves and a shorter stem than a genetically identical plant grown in a shadowed place. This physiological effect, called etiolation, is an adaptive response to a particular state of the environment, i.e., the quality and the intensity of light the plant receives. Etiolation allows the plant shadowed by neighbours competing with them in the struggle for light and to survive in spite of unfavourable conditions. Environmental changes do not only affect permanently the phenotype produced by the genotype such as morphological and/or life-history traits. They can also induce several reversible variations of the phenotype (non permanent effects) throughout the life of individuals (Lynch and Jones, 1998). Non permanent effects were particularly studied by behavioural ecologists and many studies were focused upon the fact that organisms can adaptively adjust their behaviour to the environmental heterogeneity during their lifetime.

The empirical concept of reaction norm (RN) is usually defined as "the set of phenotypes [including behavioural phenotypes (A/N)] that a single genotype produces in a given set of environments. Genotypes or individuals show phenotypic plasticity if their RN is nonhorizontal" [START_REF] Dingemanse | Behavioural reaction norm: animal personality meets individual plasticity[END_REF]. As shown in fig. 1, RN elevation and slope (plasticity) can vary among genotypes, and it turns out that some species are more plastic than others when confronted with environmental changes [START_REF] Via | Adaptive phenotypic plasticity: consensus and controversy[END_REF][START_REF] Dewitt | Costs and limits of phenotypic plasticity[END_REF]. In addition, some genes can considerably vary in their expression level, while others are more or less constant. These differences probably depend on the functions of the genes. . Phenotypic plasticity can be described as the amount of change across environmental states. The set of parameters (a vector) of a polynomial regression (dotted line) is used to describe the reaction norm (RN). The horizontal RN depicts a null plasticity (from Via et al. 1995 (modified)). (b). The relationship between phenotypic values and fitness shows the limits of plasticity in regard to one environment. A given genotype makes a phenotypic error if its mean phenotype differs from the optimum. The cost of such an error is the difference between the expected fitness and the optimal fitness (from DeWitt et al. 1998 (modified)).

The benefit of adaptive plasticity is that organisms can match across more states of the environment in comparison to organisms able to produce a single phenotype only. Currently, it is advocated by several authors that the most plastic phenotypes are harboured by organisms competing in heterogeneous environments relatively to those living in homogenous environments [START_REF] Ellers | Life history evolution in Asobara tabida: plasticity in allocation of fat reserves to survival and reproduction[END_REF][START_REF] Van Kleunen | Constraints on the evolution of adaptive phenotypic plasticity in plants[END_REF][START_REF] Weinig | Plasticity versus canalization: population differences in the timing of shadeavoidance responses[END_REF][START_REF] Snell-Rood | Toward a population genetic framework of developmental evolution: the costs, limits and consequences of phenotypic plasticity[END_REF]. It is now clear that adaptive phenotypic alters a variety of interactions between individuals and their environment through life-traits and/or behaviour : indirect interactions in multi-species assemblages, direct interactions between species such as inducible defences/offenses [START_REF] Miner | Ecological consequences of phenotypic plasticity[END_REF][START_REF] Ramos-Jiliberto | Population dynamics of prey exhibiting inducible defences: the role of associated costs and density-dependence[END_REF][START_REF] Ramos-Jiliberto | Role of inducible defences in the stability of a tritrophic system[END_REF] and niche construction [START_REF] Donohue | Niche construction through phenological plasticity: life history dynamics and ecological consequences[END_REF]. Because phenotypic plasticity is able to modify such interactions, it ultimately affects ecological processes, such as population stability, trophic relationships, population dynamics, species coexistence within communities and biodiversity [START_REF] Mouritsen | Parasites boosts biodiversity and changes animal community structure by trait-mediated indirect effects[END_REF]. Lastly, studies demonstrated that plasticity can promote stability and health of ecological systems experiencing environment fluctuations [START_REF] Miner | Ecological consequences of phenotypic plasticity[END_REF] and references therein). However, organisms able to exhibit a perfect plasticity are unlikely to exist. Such perfect organisms should produce the best phenotype in every condition of the environment. Indeed, most of organisms fail to produce the exact optimal phenotype for various reasons: inability to produce it (due to evolution constraints), inability to get reliable cues on the state of the environment and inability to pay the plasticity cost [START_REF] Auld | Re-evaluating the costs and limits of adaptive phenotypic plasticity[END_REF]. Energetic cost is a constraint which limits the evolution and the development of plasticity. In this respect, [START_REF] Dewitt | Costs and limits of phenotypic plasticity[END_REF] have listed a set of five main potential costs of phenotypic plasticity. Among them, these authors emphasize the "production cost", which should be considered only if the cost of production by plastic genomes exceeds those for fixed genotypes producing the same phenotype. This means that the variance associated to the average genome expression has an energetic price. The wider the variance, the higher the cost. However, the true nature of the undeniable cost of plasticity is still somewhat obscure and remains an open discussion. In this paper, we focus on the estimation of the best cost/benefit ratio of plasticity in the context of fluctuant (stochastic) environments. Since the energetic cost paid is positively correlated to the variance of the response to environmental fluctuations, the arising question is: what is the optimal variance of a phenotype undergoing a fluctuating environment which minimizes the cost/benefit ratio? First, we study the results obtained through the simulations of a model. Three genes are involved in a single phenotypic response driving the behaviour of a virtual wasp, which aims at maximising -by mean of a genetic algorithm -its progeny when confronted to a fluctuating environment. We find that in the course of the optimisation process, according to theoretical works in the domain, the logarithm of the determinant of the dispersion matrix of the three genes is an indicator of the convergence. Second, we build an analytical model of three genes having analogies with the previous one. This model is an attempt of generalizing the previous empirical approach. The optimal response (i.e. the minimisation of the cost/score ratio) is then found using the Lagrange's constrained optimization method, on the basis of a conjecture we made about the energy available and the covariance matrix of the 3 genes response. The analysis of the simulation results inductively corroborates the modelling assumptions of the analytical model and constitutes a qualitative cross validation of both models. On one hand, results concern some important conjectures already formulated by several authors about the origin of genetic diversity and phenotypic plasticity in various contexts. On the other hand, we argue that the method we used has some analogies with the MaxEnt algorithm (Jaynes, 1957a(Jaynes, , 1957b) ) applied to a multidimensional system and likely provide an estimator of the maximum exergy the system can own in a specific context. However, it does not constitute a rigorous application of E.T Jaynes' method, but rather a transposition of his ideas and some mathematical developments are still in course in order to identify theoretical underpinnings.

Numerical model

Model description

Animal behaviour

The model refers to the foraging behaviour of a parasitoid 2 insect (a small wasp). Some details of the algorithm are indicated in Fig. SD1. The model has been conceived as simple and neutral as possible. It is noteworthy that we did not aim at reproducing the behaviour and particular traits of a particular species but rather some general traits of thelytokous (i.e., without sexual reproduction), solitary (i.e., laying only one egg per host) and synovigenic females of parasitoid insects whose plasticity is likely a key factor of survival and fitness in fluctuant environments. It is based on the existence of a well documented trade-off between survival and reproduction resulting from competitive allocation of resources to either somatic maintenance or egg production. In synovigenic species females have the ability to mature eggs throughout their life. A dynamic control of egg load enables animals to retain some flexibility during the adult life and to minimize the risk of experiencing time or egg limitation. For simplification, we considered the case of female wasps that have an instantaneous vitellogenic activity. In addition, animals do not feed during the oviposition stage. Consequently, young females start their life with a limited amount of reserves. At initialization, a new wasp is instantiated. Immediately after the wasp has emerged, the animal (a female) has the task to lay its eggs. To achieve this task, the female has to forage within an environment in which the resource (hosts) it is looking for (e.g. eggs from another species) is aggregated into many patches dispersed within an area. The behaviour of the wasp consists in both finding patches and attacking the greatest possible number of hosts. When a wasp encounters a patch of hosts, the cumulative number of hosts it attacks (N t ) obeys to a saturation function of time t, admitting as parameter the initial number of hosts within the patch (N 0 ):

N t = N 0 (1-exp(-αt)),
(1) where α is a parameter. Thus, the velocity for attacking hosts decreases in time, following a negative exponential function. According to the Marginal Value Theorem [START_REF] Charnov | Optimal foraging: the marginal value theorem[END_REF], when its velocity has reached the average velocity calculated on the basis of the average environment richness, the wasp leaves the patch and tries to find a new one. Such an optimal patch-leaving policy was adopted since the goal of the model was to look for optimal reproductive strategies, and because most foraging animals, especially insect parasitoids, are usually behaving in good qualitative agreement with the theorem. The cycle "foraging for hosts on patches and travelling between patches" is repeated until the wasp has reached the end of its life or has exhausted its potential fecundity.

Plasticity.

Three processes driving the wasp behaviour are here represented by three real values coded on 32 bit-structures called genes. The wasp starts its life at an initial position along a linear trade-off (coded by gene G1) between its lifespan and its fecundity (fig. 2). Max lifetime and max fecundity are arbitrary fixed to 1000. Using different (realistic) values would lead to a change in scale without affecting qualitatively the results obtained. The animal is able to change its reproductive strategy throughout its life through gene G2, which defines the limits of its plasticity. The cost linearly reduces both fecundity and lifespan of the animal by means of a linear relationship with the value of G2. The wider is G2, the heavier will be the cost (in energy). The appropriateness of a non-symmetric plasticity cost was examined, but we decided to implement a symmetric cost so that the effects on the two characteristics are of the same order of magnitude and can be directly compared. In addition, the wasp can move within the range thanks to gene G3, which is a parameter of a linear estimator [START_REF] Mcnamara | Memory and the efficient use of information[END_REF], evaluated as follows:

2 An organism that lives at the expense of another (its host), impedes its growth and eventually kills it. Insect parasitoids, which are often very tiny, attack a single organism (plant or animal), from which they derive everything they need for their own growth and reproduction. One way a parasitoid does this is by laying its eggs in the body or eggs of the host insect (after Canadian Forest Service glossary).

1 3 3 ) 1 ( - - + = i i i G G µ λ µ
The animal has an initial estimate of the encounter rate of between the lowest (i.e., 0.0) host on a new discovered patch) possible instantaneous host The higher G 3 , the higher the effect of the past experience on the foraging strategy of the wasp. Every time step of its life, the animal uses range of possible values. Correlatively, the higher rate λ i , the more the animal will invest in fecundity and the less in longevity. 2)). The energetic cost is the cost of the plasticity which is proportional to G2.

Optimization process

Each wasp holds a single chromosome simulation consists in finding the vector ( wasp, i.e., the number of eggs laid throughout its life. means of a genetic algorithm and mutation rates were fixed to the number of generations 500. The con carefully verified by runs starting inter-patch travel time, probability to be born on a patch and

Environmental fluctuations

The environment in which each simulated animal has to forage is defined by travel time T (dilution of patches) and by the average number of hosts on patches ( patch quality). Each animal of a popula However, stochasticity is introduced experience environments in which each patch differs from each other by the number of hosts they own. The number of hosts per ) 50 , 100 (

2 0 = = → σ N N
. Second, the animals can randomly born inside or outside a patch with probabilities p and (1-p) travel to accomplish. These features aim at simulating the fluctuation levels of the environment. The simulation of each wasp is repeated chromosome is obtained by aver ty in a stochastic environment

The animal has an initial estimate of the encounter rate of µ 0 that was fixed to the midpoint ) and the highest (i.e., 1/t 1 , where t 1 is the time to find the first host on a new discovered patch) possible instantaneous host-encounter rate.

, the higher the effect of the past experience on the foraging strategy of the Every time step of its life, the animal uses µ i to decide which strategy to use within the range of possible values. Correlatively, the higher the instantaneous estimated host encounter , the more the animal will invest in fecundity and the less in longevity.

fecundity trade-off of a parasitoïd wasp [START_REF] Coquillard | Complex Systems: Activity-Based Modelling and Simulation[END_REF] lifespan trade-off. G2: Range of the plasticity. The wasp can move within G2 throughout its life to optimize its fitness. G3: a parameter equivalent to sensitivity to the environment richness (Eq.( 2)). The energetic cost is the cost of the plasticity which is proportional to G2.

Each wasp holds a single chromosome, which encapsulates the three genes. The goal of the simulation consists in finding the vector (G1, G2, G3), which maximises the score of the the number of eggs laid throughout its life. The score maximization is obtained by means of a genetic algorithm based on the GENITOR model [START_REF] Whitley | The GENITOR algorithm and selective pressure: why rank-based allocation of reproductive trials is best[END_REF] were fixed to 60% and 2.5% respectively. The population size is 300 and the number of generations 500. The convergence (stability of scores) of the algorithm was carefully verified by runs starting in various initial conditions (average richness of patches, obability to be born on a patch and energetic cost)

Environmental fluctuations

The environment in which each simulated animal has to forage is defined by dilution of patches) and by the average number of hosts on patches ( patch quality). Each animal of a population should experience the same environment.

is introduced in the model at two levels. First, the wasps can experience environments in which each patch differs from each other by the number of hosts they own. The number of hosts per patch is then drawn from a Gaussian distribution . Second, the animals can randomly born inside or outside a patch ). Thus, animals start their life with or without an inter travel to accomplish. These features aim at simulating the fluctuation levels of the environment. The simulation of each wasp is repeated 20 times, and the score related to each chromosome is obtained by averaging 20 scores. Lastly, 20 replicates of each simulation are 5

(2) was fixed to the midpoint is the time to find the first ounter rate. , the higher the effect of the past experience on the foraging strategy of the to decide which strategy to use within the estimated host encounter , the more the animal will invest in fecundity and the less in longevity.

al., 2010). G1: initial position of off. G2: Range of the plasticity. The wasp can move within G2 throughout its life to optimize its fitness. G3: a parameter equivalent to sensitivity to the environment richness which encapsulates the three genes. The goal of the which maximises the score of the The score maximization is obtained by [START_REF] Whitley | The GENITOR algorithm and selective pressure: why rank-based allocation of reproductive trials is best[END_REF]. Crossing-over . The population size is 300 and vergence (stability of scores) of the algorithm was average richness of patches, ).

The environment in which each simulated animal has to forage is defined by both inter-patch dilution of patches) and by the average number of hosts on patches (N 0 , i.e., the tion should experience the same environment. in the model at two levels. First, the wasps can experience environments in which each patch differs from each other by the number of hosts patch is then drawn from a Gaussian distribution . Second, the animals can randomly born inside or outside a patch start their life with or without an inter-patch travel to accomplish. These features aim at simulating the fluctuation levels of the times, and the score related to each replicates of each simulation are performed, each of them differing by the initial values of the populations' chromosomes. The best values obtained in the population are recorded at used to compute the dispersion simulation.

Results

A polynomial approach is usually used to describe the RN ( genotypes confronted to a defined set of environments this context, the genetic dispersion coefficients are related by:

T XGX Σ = ,
where the matrix X contains a polynomial series of each variable of the environment. The relation between the vector of the states of the environments and the coefficients of the polynomial function for a given g

Xg z =
Thus, we retained ln[det(Σ)] as a descriptor of the to the environment stochasticity in the simulator. i.e. H = ln[det(Σ)]. The results we present hereafter were obtained using a single of parameters (maximal plasticity cost, inter representative of the model behaviour metric of the convergence. A complete analysis of the [START_REF] Wajnberg | Optimal resource allocation to survival and reproduction in parasitic wasps foraging in fragmented habitats[END_REF]. Results show that variance for each gene drastically diminishes but the simulation there still persist residual variance is interpreted which prevents the algorithm converging

Table I. Mean and standard convergence of the genetic algorithm

Mean SD

ty in a stochastic environment performed, each of them differing by the initial values of the populations' chromosomes. The the population are recorded at each time step for each generation and dispersion matrix of the genes (see hereafter) over the 20 replicates of the A polynomial approach is usually used to describe the RN (i.e., the set of phenotypes) of genotypes confronted to a defined set of environments [START_REF] Via | Adaptive phenotypic plasticity: consensus and controversy[END_REF]De Jong, 1995) dispersion matrix Σ of character states and the matrix of polynomial a polynomial series of each variable x which measures the state of the environment. The relation between the vector of the states of the environments and the coefficients of the polynomial function for a given genotype is then: as a descriptor of the response variance of the to the environment stochasticity in the simulator. For simplicity we will call it the "H value",

The results we present hereafter were obtained using a single of parameters (maximal plasticity cost, inter-patch travel time = 100, behaviour and well suited to test the validity of the H value as a A complete analysis of the results is the subject of Results show that when the algorithm converges to the best scores each gene drastically diminishes but never vanish (fig 3). Indeed, at the end of the simulation there still persists a variance of both responses and scores (table II). This interpreted as the effect of the stochastic fluctuations o converging toward a single and optimal value of the scores.

Dispersion Matrix and Cost/score ratio. Left: convergence of the genetic algorithm to best scores ( ) of the three genes system. Right: In the same time, the Cost/Score ratio diminishes toward an optimum value. Fitted curves are only indicative. and standard deviation of genes and scores obtained after the convergence of the genetic algorithm (500 th generation of the numerical model). performed, each of them differing by the initial values of the populations' chromosomes. The each generation and over the 20 replicates of the the set of phenotypes) of [START_REF] Via | Adaptive phenotypic plasticity: consensus and controversy[END_REF]De Jong, 1995). In of character states and the matrix of polynomial

G1

(3) which measures the state of the environment. The relation between the vector of the states of the environments and the (4) variance of the genes confronted For simplicity we will call it the "H value", The results we present hereafter were obtained using a single combination patch travel time = 100, p = 0.5) are and well suited to test the validity of the H value as a the subject of another paper algorithm converges to the best scores, the . Indeed, at the end of a variance of both responses and scores (table II). This as the effect of the stochastic fluctuations of the environment, optimal value of the scores.

onvergence of the genetic algorithm to best scores (•) : In the same time, the Cost/Score ratio of genes and scores obtained after the numerical model).

Analytical model

Model description

By means of an analytical model we now attempt to generalize the results obtained in the previous section. Let {G 1 , G 2 , G 3 } be 3 genes that contribute to a phenotype in response to an environmental factor. Individual responses of the genes are modelled by means of a saturation function:

i t g y t g y i Gi Gi β α ))] ( exp( 1 [ )) ( ( max - - = (5)
where, y Gi is the response of the gene G i , max Gi y is the maximal response of G i , α i and β i are, respectively, the velocity parameter and shape parameter of y Gi and g(t) is the state of an environmental factor at the t instant. As a consequence, even if several genes share a unique α, but different shapes β i , the responses will quantitatively differ. This equation is of concave type and is one of the most commonly used to model genes' responses. It results from Eq. ( 5) that there is no time-lag between experiencing a new signal g and the corresponding phenotype. Similarly, we are not considering the case of convex functions since they usually either describe the toxicity effect of some compounds (degradation of various metabolic pathways) or are the result of an inhibitory effect of an epistatic gene.

Figure 2a exemplifies the responses of the genes to g, each response depending on α i and β i . In the following, it is considered that the velocity parameter α i is homologous to the variance (plasticity) of the gene i response: α i ∝ σ i ². The equivalence between the parameter αi and the variance σi² of the response is justified by considering that a gene which holds a high variance/plasticity of its response has, on average, a higher response to the environment signal (plasticity in the breadth of adaptation [START_REF] Gabriel | Environmental tolerance, heterogeneity, and the evolution of reversible plastic responses[END_REF]fig. 4b), whatever the value of the environment, than a gene that can only react to a small range of the environment. The energetic cost (E i ), associated to the i th gene response, according to DeWitt (1998) and [START_REF] Svanbäck | Fluctuating population dynamics promotes the evolution of phenotypic plasticity[END_REF] is conceived as the sum of several components. Although these authors distinguish five plasticity costs, we aggregated these five items into two categories: the maintenance cost min i E (i.e., sensory and regulatory mechanisms) and the cost associated to the variance of the phenotype (production costs, information acquisition costs, developmental instability costs and genetic costs). Svanback et al. use a linear function to quantify the energetic cost associated to the variance. Here, an exponential or a power function is preferred, since the variance is a second order moment. Lastly, an additional cost (ξ) paid for the response itself is considered. This cost should be equivalent to the cost paid by a fixed genotype to produce an equivalent phenotype. The resulting equation for cost takes the general form:

)) ( ( )) exp( min t g y E E Gi i i i i i i ξ α τ γ + + = , (6)
where γ,τ and ξ are positive constants.

The overall energetic cost E is then accounted as the sum of the 3 individual costs E i , over the interval of time ∆t:

)) ( ( , t g E E t i i t ∑ = ∆ (7)
We now introduce a constraint on variances, and consequently on costs, that can be understood either as a limitation of the available energy or as an elementary representation of an epistatic link between the genes: It turns out now that the problem consists in finding the best combination maximizes the response t S ∆ , minimizes the there are an infinite number of goal, we formulate the following conjecture:

Z 3 2 1 = + + α α α (Z ≥ 0) (8)

The phenotypic plasticity, measured by jointly contributing to the response to environment fluctuations particular value Z of the available energy

We recall that the determinant of a covariance matrix has the property: (Hadamard's inequality) and that Thomas, 2006); see also (maximisation) of H = ln|Σ Σ Σ Σ| subject to by means of the standard Lagrange's

Results

We partially explored the model and the constrained optimization of varying from 0.2 to 4.

Table II Parameter a ω φ {β 1 , β 2 , β 3 } {2.0, 1.8, 1} Z [0.2,…, min i E i ξ i i i ∂ , ,τ γ t ∆ ty in a stochastic environment
. Examples of the response of three genes over an interval of time ∆t = 25, submitted to environment fluctuations. Left: periodic environment (Eq. ( 9)). Right: stochastic environment (Eq. ( 10)

It turns out now that the problem consists in finding the best combination minimizes the cost t E ∆ and satisfies the constraint of possible combinations of α i for each value of formulate the following conjecture:

The phenotypic plasticity, measured by the determinant of the dispersion response to environment fluctuations, has its optimum for of the available energy which minimizes the E/S ratio.

We recall that the determinant of a covariance matrix has the property:

) and that ln|Σ Σ Σ Σ| is concave (a proof can be found in ee also fig. SD4). Consequently, we performed the optimisation subject to the constraint (8), {α 1 , α 2 , α 3 } taken as free variables, by means of the standard Lagrange's constrained optimization method. It turns out now that the problem consists in finding the best combination {α 1 , α 2 , α 3 } which the constraint (8). Notice that for each value of Z. To reach this dispersion matrix of genes has its optimum for a

We recall that the determinant of a covariance matrix has the property: (Cover and we performed the optimisation taken as free variables, to neutral values (table II) relaxing the constraint ( 8 .77, stochastic (uniform) environment; H meets the minimal cost/score ratio (E/S ≈ calculations gave identical results (Fig. SD4). The corresponding values compromise between the extreme values (Fig. SD2). The best combinations of {0.151, 0.374, 0.724}, {0.150, 0.373, 0.726} stochastic environments. These combinations e the best response in fluctuating environments at the lowest ). Figure 6d shows clearly that the relative contributions of the do not evolve linearly with Z (see also Fig. SD3). The variation depends on the shape . The maximum (H ≈ 2.28, H ≈ 4.18, stochastic 0.224). Numerical e corresponding values of E and S form a combinations of {α 1 , α 2 , α 3 } and {0.143, 0.364, should warrant the the lowest possible cost rly that the relative contributions of the depends on the shape parameter β i of the y i responses. When the system is highly constrained, the gene G2 (β = 1.8) is the most solicited. Next, its contribution diminishes in parallel with G1 (β = 2) while the contribution of G3 (β = 1.0) increases when Z is relaxed. The variations of weights both in costs and scores (γ i , τ i , δ i , ξ i ) should likely modify this pattern.

∏ = ≤ < n i ii 1 2 0 σ Σ proof can be found in

Discussion & Conclusion

Optimal plasticity, optimal energy

Results show that in the same time the numerical model converges to the best scores, the determinant of the dispersion matrix diminishes toward a residual value, which is of comparable magnitude (H = 4.61) to the one obtained with the analytical model using a Gaussian-stochastic environment (H= 4.18). However, these similar outcomes cannot lead to clear and definite conclusions. Simply, we notice that in both cases, it persists (even after the convergence of the systems) a non negligible variance. It results from the latter consideration that the system never converges to a unique solution but to a set of solutions, each of them forming the best combination of the three variances associated to genes to give the best adaptive response to a particular environmental fluctuation in time. Residual value of H can thus be understood as a description of the optimal variance for best scores at lowest cost in the context of stochastic, but limited, environmental fluctuations.

Since the "niche construction" process (Odling-Smee et al., 2003) refers to the ability of an organism to alter its environment, the plasticity of characters seems to be particularly involved in the modifications of habitats: 1.

Optimal phenotypic plasticity allows providing the best response and consequently performing a maximal direct work on the environment (e.g., seed dispersal or germination timing for plants, foraging activity or habitat choice for animals...) 2.

The minimization of the cost/benefit ratio allows saving a certain amount of energy (compared to a non-optimized organism) which can be reallocated to other physiological compartments resulting in an indirect work on the environment. Hence, according to the definition of exergy, one can assume that the optimum of phenotypic plasticity coincides with the optimum of exergy the organism can own [START_REF] Jørgensen | Towards a thermodynamics theory for ecological systems[END_REF] with respect to (i) the environmental signal the organism experiences and (ii) the genetic pathway involved in the response.

Obviously, the E/S value obtained is a theoretical optimum of an ideal organism. However, only few individuals of populations are likely to respond optimally. For several reasons: intrapopulation genetic variation, but also individuals are not in the same physiological state at a given time, to respond to environmental stimuli and therefore to devote the energy required for an optimal response. Hence, when an organism (or a population) has a very different plasticity than predicted by the H value, an accurate review of organisms living conditions must be made to determine which processes are acting to maintain them far from this value. Lastly, plasticity limits can probably lay in functional and internal constraints that reduce the benefit of plasticity in comparison with a perfect and optimal plasticity (see [START_REF] Gabriel | Environmental tolerance, heterogeneity, and the evolution of reversible plastic responses[END_REF] for a discussion about the plasticity limits, and van [START_REF] Van Kleunen | Constraints on the evolution of adaptive phenotypic plasticity in plants[END_REF] for an analysis in the context of plants plasticity).

The H value calculated is used to describe the volume in which the organism is likely to deploy its activity to ensure survival, at the lowest cost, in the limits of environmental fluctuations it encounters. In other words, the H value describes the activity (a state space) of the organism, under survival constraints, and provides a metric for the dispersion of the genes' responses, taken as random variables. While many aspects of gene expression require an accurate knowledge of molecular mechanisms, the distribution of responses of genes is largely explained by a simple model based on statistics rather than a theory of molecular biology.

Plasticity and environmental variations

The accurate knowledge of functional response equations of genes involved in a phenotype and the associated costs should allow, applying this method, to predict -ad minima-the environmental fluctuations to which may successfully face the organisms. This feature could help in predicting the potentiality of species to become invasive when artificially transplanted into exogenous ecosystems or, on the contrary, to fail to survive [START_REF] Donohue | Adaptive divergence in plasticity in natural populations of Impatiens capensis and its consequences for performance in novel habitats[END_REF]. In addition, according to [START_REF] Kaneko | Relationship among phenotypic plasticity, phenotypic fluctuations, robustness and evolvability; Waddington's legacy revisited under the spirit of Einstein[END_REF] who found a theoretic positive correlation between genetic variance and phenotypic variance over many genes, one can find through the results presented here a basis to better understand the sustainability of such a level of variation, in spite of the natural selection that should act in direction of a reduced plasticity in the context of fixed environments. Thus, according to van [START_REF] Van Kleunen | Constraints on the evolution of adaptive phenotypic plasticity in plants[END_REF] one can make the conjecture that species confronting harsh and monotonous factors close to their extremum (arid areas, cold climates, salted soils…) or mesic but stable factors (tropical forests: [START_REF] Figueiredo-Goulart | Variability in Fruit and Seed Morphology Among and Within Populations of Plathymenia (Leguminosae -Mimosoideae) in Areas of the Cerrado, the Atlantic Forest, and Transitional Sites[END_REF], should be less plastic than vicariant species accustomed to environments that are heterogeneous in space and time [START_REF] Ellers | Life history evolution in Asobara tabida: plasticity in allocation of fat reserves to survival and reproduction[END_REF][START_REF] Snell-Rood | Toward a population genetic framework of developmental evolution: the costs, limits and consequences of phenotypic plasticity[END_REF][START_REF] Tieleman | Phenotypic variation of larks along an aridity gradient: are desert birds more flexible?[END_REF].

Environmental variations as a source of genetic diversity

The results obtained here provide a restricted set of solutions close to the absolute optimum (in terms of cost and scores). This is clearly visible on fig. 6e, where acceptable solutions range approximately from Z = 1 to Z = 1.4. The numerical model does not contradict this observation since scores (= fitness) after convergence of the algorithm exhibit a residual variance (SD = 3.88). Consequently, individuals can adopt any one of nearly equiprobable solutions. Besides, it can be inferred that some small differences between the genes' responses, or transitory variations of energy availability (Z), will be of feeble consequence on their fitness (Fig. SD2). Thus, even among individuals exhibiting a quasi-optimal phenotypic plasticity with respect to their environment, there might still persist a small, but hardly reducible, genetic variance, keeping an open door for further modifications of the genes' apparatus. This is an important point, since, "…higher [genetic] biodiversity means that a wider spectrum of properties is available for survival under changing conditions" [START_REF] Jørgensen | Towards a thermodynamics theory for ecological systems[END_REF]. Hence, we agree with [START_REF] Kaneko | Relationship among phenotypic plasticity, phenotypic fluctuations, robustness and evolvability; Waddington's legacy revisited under the spirit of Einstein[END_REF] and [START_REF] Svanbäck | Fluctuating population dynamics promotes the evolution of phenotypic plasticity[END_REF] who advocate for environmental variations (among other processes) as a fundamental root of evolvability. In a same manner, Lande (2009) using a quantitative genetic model of plasticity found that major changes in the environment should first induce a substantial increase in genetic variance followed by a second phase where "the expected mean phenotype undergoes the last small fraction of the adaptation to the new expected optimum". However, the author also indicates that on a long time scale, "stabilizing selection around the expected optimum in the new average environment would re-canalize the genetic variance". But, this effect should remain a slow process. It results from these few considerations that a substantial part of genetic diversity could probably be hidden by an apparent homogeneity of the phenotypic plasticity in an ever changing world.

Lastly, we have to underline that further theoretical or empirical investigations on phenotypic plasticity should pay greater attention to the effects of both grain of environmental variation (coarse to fine-grained) and intensity of selection. Particularly, the role of these two interacting constraints on plasticity versus specialization is still in debate [START_REF] Snell-Rood | Toward a population genetic framework of developmental evolution: the costs, limits and consequences of phenotypic plasticity[END_REF][START_REF] Merilä | Variation in the degree and costs of adaptive phenotypic plasticity among Rana temporaria populations[END_REF]) and its elucidation should be of greatest importance to clarify cost, limits, dynamics and properties of plasticity.

Alternative method

An alternative to the maximization of the determinant of the covariance matrix (section 3) could have been the use of a Structurally Dynamic Model (SDM). The latter is a heuristic method recently applied with success to more than 20 ecological problems (particularly to the process of lake eutophication [START_REF] Zhang | Application of eco-exergy for assessment of ecosystem health and development of structurally dynamic models[END_REF]). In SDM, parameters of functions vary continuously in order to reproduce the adaptation of the system to forcing variables that are steadily moving through time. These changes are obtained by means of the optimisation of a goal function. In most cases, the goal function used is, among others, the eco-exergy. The main idea of SDM is to find, at each step of the simulation, the set of parameters "that is better fitted for the prevailing conditions of the ecosystem" (Jørgensen, 2009). Obviously, there are many similarities between the present work and the adaptive programming way of SDM. However, we notice that: 1. The present work mainly lies on statistics (we fitted some functional responses) rather than on detailed and realistic mechanistic relationships between the components of the system. 2. The optimization procedure we used is deterministic (Lagrange's constrained optimization method) and does not require prohibitive computing time. 3. The method presented here is unable to reproduce any dynamics and to produce prognoses such as SDMs can do. The method can only predict what should be the ideal parameters of an ideal system experiencing a fluctuant environment, and the distribution of the responses of its components.

Finally, it is clear that the choice of one of the two methods will depend on several considerations. A high complexity and an accurate knowledge of the main relationships between elements composing the system, objectives of the model such as shift of species composition and/or adaptation through time, the prognoses of the impact of environment modifications and the ability to elaborate a goal function should lead the modeller to the SDM choice.

Phenotypic plasticity and the MaxEnt algorithm

The most simple and convenient method for finding distribution functions, about which few is known, in statistical physics and information theory was proposed by E. T. Jaynes (1957aJaynes ( , 1957b) (see Fig. SD5 for a succinct description of the method and [START_REF] Martyushev | Maximum entropy production principle in physics, chemistry and biology[END_REF] for a complete review about Maximum Entropy Production Principle (MEPP) and its substantiation by Jaynes' Maximum Entropy Principle (MaxEnt)). Recently, the MaxEnt has been used to explain ecological patterns [START_REF] Dewar | Statistical mechanics unifies different ecological patterns[END_REF][START_REF] Harte | Maximum entropy and the state-variable approach to macroecology[END_REF]Azale et al., 2010;[START_REF] Banavar | Applications of the principle of maximum entropy: from physics to ecology[END_REF], biodiversity [START_REF] Shipley | From plant traits to plant communities: a statistical mechanistic approach to biodiversity[END_REF] and foodweb distribution [START_REF] Williams | Simple MaxEnt models explain food web degree distributions[END_REF]. However, some limitations in the application of the method to ecology have been pointed out [START_REF] Banavar | Applications of the principle of maximum entropy: from physics to ecology[END_REF][START_REF] Haegeman | Trivial and non-trivial applications of entropy maximization in ecology: a reply to Shipley[END_REF].

In the present case, the problem consists in finding the unknown distribution of {α 1 , α 2 , α 3 } about which we only know the energetic cost this distribution induces. Let us suppose that the genes' response is multinormal. In such conditions, H should be maximized, with respect to the distributions of the 3 genes' responses given by (4), satisfying the constraint (7), where H (expressed in nats) is the differential entropy [START_REF] Cover | Elements of information theory[END_REF] of the three genes (see below, Eq. ( 12)).

A brief description of the MaxEnt algorithm is given in supplementary information. Entropy of multivariate distributions has not been the subject of many works to date. A multivariate normal distribution admits two parameters µ (a vector), which defines its centre and a positive-definite symmetric matrix Σ Σ Σ Σ, which is the dispersion matrix of the random vector X (

3 ℜ ∈ X
in this case). The entropy of a multivariate distribution has the property: with equality iff X ~N (0, K), where n = 3 is the dimension of the distribution and K=EXX t .

Obviously, there are some strong similarities between the calculation of H = ln|Σ| we made and the differential entropy (Eq. ( 12)). In fact, both calculations lead to identical results under the assumption that responses are multinormal distributions p(S i ), which is, unfortunately, not the case since environment fluctuations g(t) are not Gaussian (excepted in one case) and genes' responses induced by the signal are calculated through saturation functions. This application can be considered as a transposition of the MaxEnt idea to a specific case, rather than a rigorous and standard application of the method. For two reasons at least: 1. We did not compute a true entropy since the dispersion matrix Σ is obtained from the gene responses and not from probability distributions. 2. Doing so, one does not make any assumption about the distributions p(S i ). Indeed, in such case, it is only supposed that the desired maximum entropy should be inferior to the maximum entropy of a multinormal-zeromean distribution, whatever the distributions under study.

The results presented here, somehow bode a new way of applying the MaxEnt, with some restrictions (concave functions) but theoretical underpinnings of such a scheme require now further developments.

Figure 1 .

 1 Figure 1. Phenotypic plasticity and fitness. (a). Phenotypic plasticity can be described as the amount of change across environmental states. The set of parameters (a vector) of a polynomial regression (dotted line) is used to describe the reaction norm (RN). The horizontal RN depicts a null plasticity (from Via et al. 1995 (modified)). (b). The relationship between phenotypic values and fitness shows the limits of plasticity in regard to one environment. A given genotype makes a phenotypic error if its mean phenotype differs from the optimum. The cost of such an error is the difference between the expected fitness and the optimal fitness (from DeWitt et al. 1998 (modified)).

Figure 2 .

 2 Figure2. The life span-fecundity trade the wasp along the fecundity-lifespan trade throughout its life to optimize its fitness. G3: a parameter equivalent to sensitivity to the environment richness (Eq.(2)). The energetic cost is the cost of the plasticity which is proportional to G2.

Figure 3 .

 3 Figure 3. Dispersion Matrix and Cost/score ratio. and corresponding H value (♦) of the three genes diminishes toward an optimum value.

Figure 4 .Figure 5 .

 45 Figure 4. Gene response and plasticity. = 4.0). The arrow indicates an example of variation due to the velocity parameter plasticity will have, on average, a smaller response than a gene owning a wider one, when confronted to a fluctuant environment over a large enough interval of time.

  the model: most of the parameters are fixed to neutral values optimization of H = ln|Σ Σ Σ Σ| is performed relaxing the constraint I. Parameters' value used for computations

Figure 6 .

 6 Figure 6. Optimal plasticity for the three genes model (a, b, c) Optimal cost and score are neither minimal nor maximal, but are intermediate values representing the best compromise. The minimal ratios obtained were: E/S = 0.224 = 43 0.224 = 43.45/193.87 (uniform-stochastic environment); E/S = 0.224 = 42.97/191.38 (Gaussian environment). (d) The contribution of each variance to the overall sum varies as Z is relaxed (vertical dashed l corresponds to the optimal response Z = 1.25). the constraint Z (scales of vertical axis are arbitrary). Maximum 1.25 (vertical line).

Figure

  Figure 6e depicts the effect of relaxing the constraint periodic environment; H ≈ (Gauss) environment) meets the calculations gave identical results ( compromise between the extreme values parameters are respectively, {0.151, 0.374, 0.724 0.743} for the periodic and stochastic environment organism to produce the best response in fluctuating of energy (fig 6a, 6b and 6c). (α i )s' do not evolve linearly with Z

(

  example of Gaussian fluctuations of the environment). Optimal cost and score are neither minimal nor maximal, but are intermediate values representing the .67/194.97 (periodic environment); E/S = stochastic environment); E/S = 0.224 = 42.97/191.38 (Gaussian-stochastic The contribution of each variance to the overall sum varies as Z is relaxed (vertical dashed line and Cost/Score ratio in function of the minimal E/S ratio (0.224) at Z =

Acknowledgments. We would like to warmly thank Drs. René Feyereisen and Eric Wajnberg (INRA, Sophia-Antipolis), Dr James Nutaro (Oak Ridge National Laboratory, TN), and Prof. R. Dewar (Australian National University of Canberra) for fructuous discussions and suggestions about this paper. Simulations were performed thanks to the free softwares: IDE cross-platform Code::Blocks, the C++ Borland™ 5.5 compiler, the MIT's Galib library (Wall, 2007) and the environment for statistical computing R.