The environmental zero‐point problem in evolutionary reaction norm modeling

Abstract There is a potential problem in present quantitative genetics evolutionary modeling based on reaction norms. Such models are state‐space models, where the multivariate breeder's equation in some form is used as the state equation that propagates the population state forward in time. These models use the implicit assumption of a constant reference environment, in many cases set to zero. This zero‐point is often the environment a population is adapted to, that is, where the expected geometric mean fitness is maximized. Such environmental reference values follow from the state of the population system, and they are thus population properties. The environment the population is adapted to, is, in other words, an internal population property, independent of the external environment. It is only when the external environment coincides with the internal reference environment, or vice versa, that the population is adapted to the current environment. This is formally a result of state‐space modeling theory, which is an important theoretical basis for evolutionary modeling. The potential zero‐point problem is present in all types of reaction norm models, parametrized as well as function‐valued, and the problem does not disappear when the reference environment is set to zero. As the environmental reference values are population characteristics, they ought to be modeled as such. Whether such characteristics are evolvable is an open question, but considering the complexity of evolutionary processes, such evolvability cannot be excluded without good arguments. As a straightforward solution, I propose to model the reference values as evolvable mean traits in their own right, in addition to other reaction norm traits. However, solutions based on an evolvable G matrix are also possible.

model describes how a multivariate individual phenotype y i,t is expressed as a linear or nonlinear function of quantitative traits z 0,i,t and a continuously varying multivariate developmental environment (environmental cue) u t , Here, z 0,i,t may be the individual parameter vector as function of time t (generations) in a parametrized model of the reaction norm, or alternatively the individual phenotypic values at discrete index environments. Interpolation between index environments results in a function-valued or infinite-dimensional individual reaction norm model y i,t = γ(u t − u ref ) (Kingsolver, Gomulkiewicz, & Carter, 2001;Kirkpatrick & Heckman, 1989;Kirkpatrick, Lofsvold, & Bulmer, 1990). The reference environment is often set to u ref = 0 (Gavrilets & Scheiner, 1993a,b;Lande, 2009), but that disguises the problem at hand. Second, the individual fitness function is where θ t is the vector of phenotypic expression that maximizes fitness in the given generation. Note that in the univariate and linear case, the covariance between u t and θ t determines the mean reaction norm slope in a stationary stochastic environment (Ergon & Ergon, 2017;McNamara, Barta, Klaassen, & Bauer, 2011). Third, the state equation that propagates the mean trait values may under given assumptions be the multivariate breeder's equation (Lande, 1979) Equation (3) is based on the assumption that the phenotypic traits can be split into two mutually independent and multinormally distributed parts, z 0,i,t = x 0,i,t + e 0,i,t , with the covariance matrices G = E x 0,i,t −x 0,t x 0,i,t −x 0,t T and E = E e 0,i,t −ē 0,t e 0,i,t −ē 0,t T , respectively. As a consequence also z 0,i,t is multinormally distributed, with the covariance matrix P = E z 0,i,t −z 0,t z 0,i,t −z 0,t T . I will here assume P and G to be constant, which is common in theoretical work (e.g., Lande, 2009), although it is unrealistic over longer time periods (Steppan, Phillips, & Houle, 2002). I will assume populations with non-overlapping generations, where all individuals live in the same time-varying environment, and make standard assumptions necessary for the multivariate breeder's equation (3) to be valid (Lande, 1979). For analytical purposes, expressions for mean values ȳ t and W t can be found from equations (1) and (2), but they are not needed for simulations.
The fundamental insight is that the reference environment in reaction norm models is an inherent part of the population state, independent of the actual environment where the individuals develop. The state of the population thus determines which environment it is adapted to, that is, where the expected geometric mean fitness is maximized (the location of the adaptive peak or the growth rate peak along the environmental axes). The environment the population is adapted to is, in other words, an internal population property, independent of the external environment. It is, however, only when the external environment coincides with the internal reference environment, or vice versa, that the population is adapted to the current environment. This is formally a result of state-space modeling theory, which is an important theoretical basis for quantitative genetics evolutionary modeling. As a consequence, the reference environments should be modeled as part of the evolutionary models (1) to (3). How this should be done is an open question, where the best answer may depend on the problem under study. One alternative is to let u ref be a function of an evolvable G matrix (e.g., Arnold, Bürger, Hohenlohe, Ajie, & Jones, 2008). That would give a complex solution, especially in the multivariate and nonlinear case, and this alternative is not further discussed (except in a simple numerical example in Section 4). As a straightforward solution, I propose that the reference environment vector may be modeled as a vector z c,t of mean traits in their own right, just as other reaction norm traits. Equation (3) must accordingly be augmented with the z c,t state variables. The details of this for parametrized models are developed in Section 2, while augmented function-valued models are discussed in Appendix 1.
Whether the reference traits in z c,t are evolvable is also an open question, but considering the complexity of evolutionary processes, such evolvability cannot be excluded without good arguments. Also note that evolvable reference traits may be combined with an evolvable G matrix.
The idea of an evolvable reference trait was introduced in Ergon and Ergon (2017), but then based on biological arguments, as a result of the novel idea of a perception trait as a means of relaxing constraints on the evolution of reaction norms. A main purpose of the present article is to show that the plasticity reference environment not only may be modeled but that it in principle must be modeled, in one way or another, as part of the quantitative genetics state-space model (although this is not necessary if the reference environment is not evolvable).
The reference environment vector z c,t is closely related to the environment the population is adapted to, which we may denote u 0 . As discussed in detail for the special case in Ergon and Ergon (2017), an unsymmetrical distribution of the phenotype y results in a difference between z c,t and u 0 , but at equilibrium in a stationary stochastic environment the expected deviation is independent of the mean values μ U and μ Θ of u t and θ t , respectively.
Under the assumption that some elements in the environmental reference trait vector are genetically variable, these elements must be included in the state equation (3), or its function-valued counterpart. In Section 2, I show how this can be done for multivariate and nonlinear parametrized reaction norms. If all elements in the reference environment have zero genetic variance in the population, they can without consequence be set to zero, and this is thus an implicit assumption in traditional reaction norm models.
As discussed in Ergon and Ergon (2017), an important result of a fully evolvable plasticity reference environment is the property of complete genetic assimilation, by which "selection can act in such a manner as to turn an environmentally stimulated phenotype (i.e., plasticity) into a fixed response to prevalent environmental conditions (assimilation)" (Pigliucci & Murren, 2003). I here use the term

ERGON
"complete genetic assimilation" as in Ergon and Ergon (2017), to describe the evolutionary scenarios where, after an abrupt environmental change, there is an initial increase in phenotypic plasticity, after which the mean plasticity is reduced and the environment range, or value, to which the population is adapted moves toward the current mean environment. This entails that all elements in the reference environment vector have genetic variability, such that they are evolvable.
A major difficulty of the approach with evolvable reference traits is to find empirical measures of these latent parameters. In the linear and univariate example in Ergon and Ergon (2017), for example, individual reference traits z c (horizontal reaction norm variation) cannot be distinguished from z a (vertical reaction norm variation) by means of a static breeding design (Hill, 2010;Thompson, Brotherstone, & White, 2005). Different values of the variance G cc of z c will, however, give different dynamical responses to environmental variations, and assuming that the variance G aa of z a is known this can be used to identify G cc . It is also possible to identify several parameters (Appendix 2).
In Section 3, I simulate a set of linear reaction norms, to clarify why the environment u 0 where the phenotypic variance has a minimum must be seen as a population characteristic. I also include a simulation example with multivariate and nonlinear reaction norms, where the environment changes in a sudden step, and where the property of complete genetic assimilation is demonstrated. In a third simulation example, I show the effect of evolvable environmental reference values on the results in Chevin and Lande (2015), regarding the plastic response in a population system with a single phenotype, and two correlated environmental cues. Finally, I include a discussion in Section 4.
In Appendix 1, I show that the plasticity reference environment needs to be modeled also in function-valued models, and how that can be done for univariate and nonlinear reaction norms based on environmental index values. I also describe two additional problems in such cases.
In Appendix 2, I finally present a preliminary example showing how the variances of and covariances between quantitative reference traits may be identified from dynamical experiments.
An example Matlab code for the extended Chevin-Lande simulation is provided in Data S1.

| Background state-space theory
As a background and reference for the theoretical development, I include a summary of the underlying state-space theory for discretetime systems. Caswell (2001, Ch. 3) refers to Zadeh's formal theory of state (Zadeh & Desoer, 1963), but state-space modeling of dynamical systems is older than that. Of special historical importance is the seminal paper of Kalman (1960), concerning the discrete-data linear filtering problem (Kalman filtering), although linear statespace models are special cases.
The starting point for a general discrete-time state-space model is the idea of an abstract discrete-time system that interacts with its environment through a vector φ t of input variables and a vector y t of response variables. A vector x t of variables that takes its values in some set X (a state-space) is a state vector if it satisfies the following two requirements: 1. There exists a function g (·) that uniquely determines the response at any time t as a function of the input and the state at t,

2.
There exists a function f (·) that uniquely determines the state at any time t as a function of the state at any earlier time t 0 and the input sequence from t 0 to t − 1, for any t 0 and sequence φ 0 , φ 1 ,…, φ t−1 , that . From this follows that x 1 = f(x 0 , φ 0 ), and generally that x t at any time t can be propagated one step forward in time according to (Åström & Murray, 2008) The function g (·) is known as the output or observation function, and the function f (·) as the state function, while x t is the state. At t = t 0 the state variables will have or be given some initial values, but from then on all information from the past is carried by the state variables. It should be noted that any current state may be the result of a large number of different initial states and input sequences, especially if t 0 is far back in time, and the initial state cannot therefore be reconstructed from the current state without detailed knowledge of the entire input sequence. Also note that the excitation φ t may be a combination of deterministic and stochastic signals, and that the functions g (·) and f (·) may include different parts of a common input vector φ t .
When g (·) and f (·) are linear, and when the stochastic part of φ t is white (no autocorrelation) and normally distributed, the optimal mean value x t can be found from y t using the Kalman filter, such that (Lewis, Xie, & Popa, 2008;Newman et al., 2014). In the general case, estimates of the distribution of x t can be found from y t using the Chapman-Kolmogorov

| Evolutionary state-space models
Assuming sufficient genetic variation, the mean phenotypic values in a population will evolve when the environment varies from generation to generation. As summarized in the Introduction, mathematical modeling of this evolution for plastic organisms involves a statespace model, which assuming non-overlapping generations require three equations. First, equation (1) describes how a multivariate individual phenotype y i,t is expressed as a linear or nonlinear function of quantitative traits z 0,i,t and a continuously varying developmental environment (cue vector) u t . Second, equation (2) describes how the individual fitness depends on the difference between the phenotype y i,t and the vector θ t of phenotypic expressions that maximizes fitness in the given generation. Third, the state equation may under given assumptions be the multivariate breeder's equation (3) (Lande, 1979).
When equation (1) is compared with the general state-space output function (4), it is apparent that the environmental reference vector u ref must be part of either the current state or the current input. As equation (4) describes how the abstract discrete-time system interacts with the current environment through the vector φ t of input variables, and as a reference environment possibly far away from the current environment cannot be part of the current input, it must necessarily be an inherent part of the current state of the population (as illustrated in Figure 1 in Section 3). The current , which leaves u t as the current input in equation (4). Note, however, that also θ t in the fitness function (2) is an input variable, such that the total current input In traditional reaction norm models, the reference environment is assumed to be the same for all individuals in the population, and the current mean state is then , that is, the reference environment is in principle a population state variable, although it is implicitly assumed be constant. The environment the population is adapted to, is, in other words, an internal population property, independent of the external environment. It is, however, only when the external environment coincides with the internal reference environment, or vice versa, that the population is adapted to the current environment.
Again, note that u ref cannot be a part of the current input, which according to equation (4) interacts with the system. The state variable u ref thus determines which environment the population is adapted to, whether it coincides with the current environment or not.
Any population state variable must be modeled as a population mean value, a variance or a higher order statistical moment, or functions of the statistical moments. As we must assume that the population may be adapted to different stationary stochastic environments, independent of constant G and P matrices, and as the elements in u ref must have the same dimensions as the elements in u t (e.g., temperature and salinity), the remaining choice is a mean trait vector, which we may denote z c,t . Note that u ref should be modeled in this way also when it is set to zero, and that we in general must assume that z c,t may be evolvable. As described in Section 4, this way of modeling a possibly varying input is natural also from an engineering control point of view.
As mentioned in Section 1, we could alternatively let u ref be a function of an evolvable G matrix (e.g., Arnold et al., 2008), but that possibility is not discussed further in this article (except in a simple numerical example in the Section 4). As also mentioned in the Introduction, an  (1). Assuming that z c,i,t , just as z 0,i,t , can be split into two independent and multinormally distributed parts, z c,i,t = x c,i,t + e c,i,t , and that the additive genetic covariance matrix in z c,t will be evolvable. This results in a dynamical reference environment, which in a stationary stochastic environment will evolve into an equilibrium.
With u ref = z c,i,t , the model (1, 2, 3) will according to the multivariate breeder's equation result in the augmented state-space model where , and thus a constant mean state variable z c,t+1 =z c,t . In that special case, we may without further consequences set In case only some of the traits in z c,i,t have genetic variability, only such traits should be included in equation (7), while the others may be set to zero. In equation (7), W i,t and W t are still computed from equation (2). Evolution in a stationary stochastic environment will lead to an equilibrium, where E cov(W i,t ,z 0,i,t = 0 and E cov(W i,t ,z c,i,t = 0, that is, where the expected selection gradient is E t = 0.

| Parametric reaction norm modeling
With z 0,i,t split into elevation traits z a,i,t and slope and shape traits z b,i,t , the reaction norm function in equation (6) becomes Following Gavrilets and Scheiner (1993a), this function can be approximated by a power series in terms of the components of q environmental cues, with p different products of u 1,t − z c,1,i,t , length m × p. Note that all of z a,i,t , z b,i,t and z c,i,t may have independent additive genetic and non-additive parts. When equation (6) is replaced by equation (9), equation (7) must be replaced by The total number of state variables is thus m + m × p + q, where q is the number of environmental cues. Note that the system (9, 10) has the external references μ Θ , μ U , and cov(U,Θ) through the fitness function (2). It follows from Ergon and Ergon (2017) that a symmetric phenotypic distribution p (y) at equilibrium in a stationary stochastic environment results in E z a,t = and E z c,t = U , while an unsymmetrical p(y) leads to deviations from μ Θ and μ U . These deviations will, however, be independent of the actual values of μ U and μ Θ , such that a positive definite matrix G cc gives complete genetic assimilation in any stationary stochastic environment. It also follows from Ergon and Ergon (2017) and McNamara et al. (2011), that the mean slope values around the origin in a stationary stochastic environment is a function of cov(U,Θ).

| Adaptive peak as population characteristic
Theoretical and simulation results for a simple linear example system with an evolvable plasticity reference environment are discussed in detail in Ergon and Ergon (2017). Here, I take a closer look at the linear reaction norms in that example, to show why the environmental cue u 0 where the phenotypic variance is minimized, that is, the location of the adaptive peak, is a population characteristic. In the example system simulated in Ergon and Ergon (2017) which with c = 0 and G cc = 0 gives the two-trait reaction norm in Lande (2009). As shown in Ergon and Ergon (2017), the environment where the phenotypic variance has a minimum is that is, u 0 =c for G bc = G ab = 0.This implies that u 0 is a population property, which as shown in Figure 1 may be located far away from the current external environment. (10)

| A multivariate and nonlinear case
As also discussed in Ergon and Ergon (2017), as well as in the Introduction, an important consequence of an evolvable reference environment is complete genetic assimilation in any stationary environment. Here, I simulate a multivariate and nonlinear system, where complete genetic assimilation as defined in the Introduction takes place. Figure 2 shows step response phase portraits, that is, ā 1 = f c 1 and ā 2 = f c 2 , for a system with the individual reaction norm model with correlated cues u 1 and u 2 , and with independent and zero mean white noise components e 1 and e 2 . The fitness function was with correlated values of θ 1 and θ 2 . The state equation (10) Note that the plots show that complete genetic assimilation takes place. Figure 3 shows the corresponding mean plasticity slope plots.
Note that only b 11 is different from zero in a stationary stochastic environment, which may have implications for the possibilities to find parameter values from collected data (see Section 4 and Appendix 2).

| Extended Chevin-Lande example
More generally, evolvable reference environments will have profound effects on all types of evolutionary modeling involving reaction norms. Here, I show how it will affect the analysis of Chevin and Lande (2015), regarding how reaction norm slope values respond to correlated multiple environmental variables. Chevin and Lande studied the plastic response in a population system with a single phenotype and two environmental cues (environments of development) u 1,t = ε 1,d,t and u 2,t = ɛ 2,d,t , and the phenotypic expression that maximizes fitness θ t = ε 1,s,t +ε 2,s,t , where ɛ 1,s,t and ɛ 2,s,t are the environments of selection, and where μ θ = B μ U 1 +μ U 2 . They used the traditional approach with reference environments equal to zero, that is, with an individual reac-

F I G U R E 2
Step response phase portraits, that is, ā 1 = f(c 1 ) and ā 2 = f(c 2 ), for a system with the individual reaction norm model (12) and fitness functon (13), with steps in μ U 1 and μ U 2 from 0 to 6, and in μ Θ 1 and μ Θ 2 from 0 to 12, applied at t = 5,000 generations. The simulation ended at t = 10,000 generations. The G matrix was diagonal with G a 1 a 1 = G a 2 a 2 = G c 1 c 1 = G c 2 c 2 = 0.5 and G b 11 b 11 = G b 12 b 12 = G b  = 1.6, cov θ 1 ,θ 2 = 0.05, cov u 1 ,θ 1 = cov u 2 ,θ 1 = cov u 1 ,θ 2 = cov u 2 ,θ 2 = 0.2, and ω 2 = 10 from equal initial values b 1 ∕B =b 2 ∕B = 1. The reason for the different final slope values is that the two cues are correlated, and also correlated with the phenotypic expression θ that maximizes fitness, and the main point in the paper was thus that interpretation of the reaction norm slopes must take these correlations into account. traits. The G cc matrix was diagonal. I used σ 2 e = 0.5 in all simulations (as in Lande, 2009). See Data S1 for Matlab code.
The simulation results in Figure 4 show that interpretation of the reaction norm slopes also must take the variances (and covariance) of the traits c 1 and c 2 into account. For G c 1 c 1 = G c 2 c 2 = 0, the results for b 1,t ∕B and b 2 ∕B are the same as in a simulation using the method in Chevin and Lande (2015) (Figure 1a). With increased values of G c 1 c 1 (blue) and b 12 (magenta), and lower panel shows b 22 (blue) and b 23 (magenta). All initial parameter values were set to zero F I G U R E 4 Evolution of scaled mean reaction norm slopes for the system in Equations (15), (16) and (10) in constant mean environments, from initial values equal to one to final stationary values b 1,∞ ∕B and b 2,∞ ∕B (blue). The G matrix was block-diagonal with G aa = 0.5, G b 1 b 1 = G b 2 b 2 = 0.04, G b 1 b 2 = 0.01 and G c 1 c 2 = 0, while G c 1 c 1 and G c 2 c 2 varied. The residual variance was σ 2 e = 0.5 in all simulations. Left panels show results for G c 1 c 1 = G c 2 c 2 = 0, that is, for the case in Chevin and Lande (2015) (Figure 1a). The central and right panels show results for G c 1 c 1 = G c 2 c 2 = 0.5 and G c 1 c 1 = G c 2 c 2 = 5, respectively. For comparison, results using the Chevin-Lande algorithm with G c 1 c 1 = G c 2 c 2 = 0 are included in the left panels (magenta, not easy to see as it is overlapped by blue curve) and G c 2 c 2 , the final absolute values of b 1,t ∕B and b 2,t ∕B were reduced.
Very large values of G c 1 c 1 and G c 2 c 2 gave b 1,t ∕B → 0 and b 2,t ∕B → 0 for t → ∞.

| D ISCUSS I ON
The main point in this article is that the plasticity reference environment u ref is a population characteristic, that ought to be modeled as such, and this is the case also if it is set to zero. Under the assumption of constant additive genetic and phenotypic covariance matrices, the remaining choice is to model u ref as a vector z c of mean traits. The corresponding additive genetic covariance matrix G cc may be zero, and we may then set u ref = 0. However, if G cc ≠ 0, at least some of the "reference traits" will evolve in a changing environment, and they must then be included in the augmented state equation (7).
One may ask why not the covariance matrices G and P also should be modeled and included as state variables in the augmented state equation (7), and the answer is yes, in principle they should. In such cases, evolvability of these matrices cannot be based on individual selection, but on, for example, mutations. Here, however, I assume that G and P are constant and not evolvable, such that augmentation with these matrices is not necessary. See Arnold et al. (2008) for a review of empirical, analytical, and simulation studies of the G matrix, with a focus on its stability and evolution.
The biological mechanism behind evolvable "reference traits" may be that individuals perceive the environment differently, as discussed in Ergon and Ergon (2017), and we could accordingly introduce individual "perception traits" z c . As shown, such perception traits may be used also in multivariate and nonlinear cases, leading to parametrized models according equations (2), (9), and (10). As shown in Appendix 1, perception traits may be used also in models based on index environment phenotypes, which through interpolation leads to function-valued models. In such models, however, G cc > 0 leads to non-normal distributions, which is in conflict with the assumptions behind the multivariate breeder's equations (7) and (10). Another added difficulty is that the individual state variable z c,i,t does not fit into a function found through interpolation between phenotypic index traits z 1,i,t to z r,i,t . A similar problem in a life-history trait setting is discussed in Irwin and Carter (2013).
The state-space model (9,10) could have been formulated just as a generalization of the model in Ergon and Ergon (2017), based on biological arguments for perception traits. In addition to that, however, my intention has here been to show that, independent of these arguments, modeling of the reference environment is in principle necessary from a basic state-space modeling point of view.
The most important result from a practical point of view is that population systems with a positive definite covariance matrix G cc obtain complete genetic assimilation in any stationary stochastic environment, as discussed in the Introduction. This means that the reaction norms at equilibrium after a change from one mean environment to another will be shifted to the new environment without any change in slope and shape. The adaptive peak, as determined by the state of the population, thus moves such that the population becomes adapted to the new environment. This movement is illustrated in a phase plane plot in Ergon and Ergon (2017), as well as in Figure 2. Long after the change in mean environment, the complete genetic assimilation will return the mean fitness to its original value, which is an essential difference from the partial genetic assimilation obtained in Lande (2009). More generally, the mean phase space position values z a and z c in equation (10) will evolve to new equilibrium values, while the mean slope and shape values z b after a transient period will return to the original values. As a result, the dynamical responses to variations around the mean of a stationary stochastic environment, that is, around an adaptive peak, will be independent of the environmental location of the adaptive peak. This is demonstrated in Figures 2 and 3 in Section 3. In practice, however, complete genetic assimilation to any environment must be limited by biological constraints, plasticity costs etc.
As an alternative to the modeling of the reference environment As mentioned in the Introduction, a main difficulty appears to be to find estimates of G cc from data. With linear reaction norms, it is theoretically impossible to find G cc from data collected at stationarity, but as discussed in Ergon and Ergon (2017), signs of G cc ≠ 0 will show up in transient situations. For the simple example in Ergon and Ergon (2017), it is in fact possible to find G cc from dynamical experiments, as used in engineering control system identification (Appendix 2). A more general application of such methods on evolutionary problems is an interesting area for future research.
It is also interesting to note that there may exist mean reaction norm slope and shape parameters that are different from zero only in dynamical situations, as demonstrated in Figure 3. Although the corresponding individual parameters will be different from zero also in a stationary stochastic environment, this may make it difficult to find the corresponding covariance parameters from data collected at stationarity. In such cases, these parameters may possibly be found using dynamical identification experiments as introduced in Appendix 2.
Under the assumption that all individuals develop in the same environment, genetic assimilation leads to good tracking properties, and thus good adaptation to slowly changing environments. This may reduce the need for nonlinear reaction norms, and also the details of this is an interesting area for future research.

ERGON
Finally note that modeling of a constant reference cue as an undriven discrete-time integrator as in equations (7) and (10), with G cc = 0, has an interesting parallel in engineering control applications. Such modeling of a constant system input is there used to achieve model based integral control through state feedback, which assures that the stationary control error is zero also if the constant input is unknown (Friedland, 1986). It also makes it possible to follow an input ramp function without an ever increasing error. The similar effects of the state-space models (6, 7) and (9, 10), with G cc ≠ 0, is the genetic assimilation in any stationary stochastic environment, as described above, and good tracking properties when the environment changes slowly.

ACK N OWLED G M ENTS
I thank Torbjørn Ergon for comments on an early version of the manuscript, and an anonymous reviewer for comments on a recent version. I also thank University College of Southeast Norway for support and funding.

CO N FLI C T O F I NTE R E S T
None declared.

Modeling based on index environment phenotypes
In order to show that the plasticity reference environment needs to be modeled also in function-valued models, I here consider models based on index environment phenotypes. Such models lead to function-valued models through interpolation between the index environments (Kingsolver et al., 2001;Kirkpatrick & Heckman, 1989). I also describe two additional problems in such cases. For clarity, I assume a univariate individual phenotype y i,t and a univariate environmental cue u t .
In an index environment model, the phenotypic values y i,t in Equation (1) is defined as the individual phenotypic values at r discrete index environments, where γ is the in general nonlinear reaction norm function, and where z c,i,t is the individual reference trait (which is set to zero in traditional models). The individual phenotypic values are also used as individual traits, that is, z 1,i,t = y 1,i,t etc., and these traits have a phenotypic covariance matrix P r = cov(y i,t ,y i,t ), and a corresponding additive genetic covariance matrix G r . Setting z c,i,t = 0, the multivariate breeder's Equation (3) would thus lead to When z c,i,t ≠ 0, this state equation must be augmented into where G rc and P rc are the covariance matrices of the vector y i,t augmented with z c,i,t . This raises two problems. First, with z c,i,t ≠ 0 the traits z 1,i,t = y 1,i,t etc. will not be normally distributed, even if the reaction norm has underlying normally distributed parameters, which is in conflict with the assumptions behind the multivariate breeder's equation (Lande, 1979). Equation (A3) will therefore be more of an approximation than it otherwise would be. Second, the state variable z c,i,t does not fit into a function found through interpolation between z 1,i,t to z r,i,t . A similar problem in a life-history trait setting is discussed in Irwin and Carter (2013). A possible solution is to assume that z c,i,t is independent of z 1,i,t to z r,i,t , and model the evolution of z c,t independently, that is, to use Equation (A2) combined with

APPENDIX 2
Preliminary example of evolutionary system identification System identification is a mature discipline in the engineering control community, with prediction error methods developed during the 1980's (Ljung 1999), and subspace methods from the 1990's and later (Qin 2006). For evolutionary system identification, predictor error methods of the output error (OE) type is a straightforward choice.
Here is an example of the OE prediction error method applied on an evolutionary system identification problem. Assume a system essentially as in Ergon and Ergon (2017) For identification of several unknown parameters, a better search method is needed. This requires experimental data that are informative enough, but it also requires a theoretical identifiability analysis (it may not be theoretically possible to identify all parameters). Also note that we must assume a model, i.e., a linear or nonlinear reaction norm, a fitness function, and a covariance structure.
The applicability of dynamical system identification methods in an evolutionary setting remains to be investigated.

F I G U R E A 1
The input function θ t = μ Θ + v θ,t (upper panel), and the model output for G cc = 0 (magenta), G cc = 0.5 (green) and G cc = 1 (blue). The green curve is also the output from the assumed evolutionary system (A5) to (A7) itself F I G U R E A 2 Identification of G cc = 0.5 based on model with population size N = 10000