Phenotypic plasticity as a cause and consequence of population dynamics.

Predicting complex species-environment interactions is crucial for guiding conservation and mitigation strategies in a dynamically changing world. Phenotypic plasticity is a mechanism of trait variation that determines how individuals and populations adapt to changing and novel environments. For individuals, the effects of phenotypic plasticity can be quantified by measuring environment-trait relationships, but it is often difficult to predict how phenotypic plasticity affects populations. The assumption that environment-trait relationships validated for individuals indicate how populations respond to environmental change is commonly made without sufficient justification. Here we derive a novel general mathematical framework linking trait variation due to phenotypic plasticity to population dynamics. Applying the framework to the classical example of Nicholson's blowflies, we show how seemingly sensible predictions made from environment-trait relationships do not generalise to population responses. As a consequence, trait-based analyses that do not incorporate population feedbacks risk mischaracterising the effect of environmental change on populations.


I N T RODUC T ION
Understanding how mechanisms of individual variation act upon populations is key to predicting how changes in the biotic and abiotic environment alter population processes. Phenotypic plasticity is the ability of individual genotypes to produce different phenotypes when exposed to different environmental conditions (Fusco & Minelli, 2010) and has been shown to be a mechanism by which species respond to climate change (Boutin & Lane, 2014;Crozier & Hutchings, 2014;Seebacher et al., 2015;Stoks et al., 2014). There is evidence that species exhibiting high levels of phenotypic plasticity are more successful at spreading across environmental gradients (Hahn et al., 2012;Szabó et al., 2018), and it is predicted that phenotypic plasticity contributes to determining the outcome of interspecific competition (Buskirk & Mccollum, 2016;Palkovacs & Post, 2009). Quantifying phenotypic plasticity in individuals is generally straightforward, but it is often more difficult to measure the effects on populations (Merilä & Hendry, 2014;Valladares et al., 2006). It is theorised that phenotypic plasticity contributes to the occurrence of seemingly paradoxical population dynamical behaviours such as the paradox of enrichment, whereby an increase in available resources causes a destabilisation of population's dynamics (Miner et al., 2005), and the hydra effect where an increase in per capita mortality results in increased population density (Cameron & Benton, 2004). Disentangling the complex network of inter-dependent individual and population processes necessary to demonstrate how phenotypic plasticity contributes to species responses to environmental change is inherently difficult using existing model frameworks (Forsman, 2015). In particular, the population dynamic consequences of phenotypic plasticity often manifest as delayed-density dependence (Beckerman et al., 2002) which is known to cause cryptic dynamical behaviours (Lima et al., 1999;Pedraza-Garcia & Cubillos, 2008).
Despite the potentially complex relationship between individual variation and population response, environment-trait relationships observed in individuals are routinely employed to predict the outcome of population processes (Figure 1). For example, in epidemiology environmental-trait relationships are used in parameter-based approaches for calculating the basic reproduction number, R 0 (Brand et al., 2016;Mordecai et al., 2017;Parham & Michael, 2010) (expected number of secondary cases produced by a single infection in a completely susceptible population). This implicitly assumes that variation observed in a population's trait distribution is independent of environmental stressors and population dynamics, such that an averaged trait value suitably represents the population at any given time and location (Liu-Helmersson et al., 2016). This is the mean-field approach and there is an increasing body of evidence that this approach under-represents the importance of variation between individuals and community structure in population ecology (Cator et al., 2019;Fox & Kendall, 2002;McGill et al., 2006;Violle et al., 2012). Consideration of purely stochastic forms of variation has demonstrated that the outcome of population processes such as species persistence do not always follow meanfield predictions (Hart et al., 2016;Morozov et al., 2013). In contrast to environmental or demographic noise, individual variation caused by phenotypic plasticity has a strong mechanistic component, and so can and should be suitably accounted for (Nylin & Gotthard, 2002).
To describe the effect phenotypic plasticity has on population dynamics it is key to link trait and effect mechanistically combining empirically derived relationships with theoretical methods. For populations with distinct generations this can be addressed using methods such as integral projection models (IPMs) which are widely used to represent trait variation within populations (Childs et al., 2003;Kuss et al., 2008). IPMs impose observed trait distributions upon populations and map changes in these distributions forwards in time to predict long-term and general trends. Due to their simplistic F I G U R E 1 Current predictive frameworks typically use environment-trait relationships, such as reaction norms, to predict population responses without consideration of how population processes may alter the traits individuals express. Our framework incorporates environment-trait relationships that interact with population dynamics and trait distributions. This allows the framework to account for the effect of interaction between environment, trait, and population as experienced by many organisms in our predictions of population processes representation of plasticity and ease of parameterisation, IPMs are used to make cross-taxa syntheses on global species trends (Salguero-Gómez et al., 2018). However, there is a recognised need for modelling approaches that can consider the effect of detailed intra-generation dynamics alongside inter-generational dynamics on the expression of phenotypic plasticity and so can predict how these feedback to alter population processes (Bolnick et al., 2011;Hendry, 2016;Johnston et al., 2019;Lion, 2018;Lipowsky et al., 2015;Lloyd-Smith et al., 2005;Sgrò et al., 2016;Turcotte & Levine, 2016;Violle et al., 2012).
Here, we propose a novel general mathematical framework that links experimentally derived environmenttrait relationships to well-parametrised stage-structured population models that allow trait distributions to emerge from population-trait-environment interactions ( Figure 1). We utilise a continuous-time stage-structured modelling approach, widely used to model organisms with multiple distinct life stages (Murdoch et al., 2003), adapted to represent the persistent and delayed effects of phenotypic plasticity across multiple developmental stages. By using our framework to represent mechanisms of individual variation in response to environmental change, we show that even simple forms of phenotypic plasticity can lead to complex population dynamical responses that previous approaches overlook. This is demonstrated by an application of our framework to a classical population ecology study, Nicholson's blowflies (Nicholson, 1957), where it has been hypothesised that previously unexplained population dynamics can be attributed to phenotypic plasticity. This application reveals a rich set of counter-intuitive population-dynamical behaviours caused by the interaction between phenotypic plasticity and population dynamics.

M AT ER I A L S A N D M ET HOD S
We present a modelling framework that dynamically links the expression of phenotypic plasticity in individuals to population dynamics. We combine a continuoustime stage-structured population model in the form described in Gurney et al., (1983) (as in Figure 2a), widely used to predict the population dynamics of interacting life-stages (e.g. Gurney et al., (1980)), with a set of empirically-derived reaction norms. Models created using our framework are systems of stage-phenotypically structured delay-differential equations, within which, we track cohorts of individuals based on their cumulative environmental experience. Each cohort is then associated with a unique phenotype. (Figure 2b). Within our framework an individual's phenotype may consist of multiple traits varying in response to multiple environmental factors, both current and historic. This creates a dynamic phenotypic structure that allows multiple phenotypes to be represented within a population simultaneously. By using this approach, the distribution of traits expressed within a population is not assumed, as it is in an IPM (Merow et al., 2014), but instead emerges as a feature of empirically verified mechanistic processes. This links individual-level variation in life-history traits to population-level response and so represents the effects of phenotypic plasticity on populations. Our approach is able to represent both intra-and inter-generational forms of phenotypic plasticity, in response to both instantaneous and delayed environmental conditions. Moreover, we can track the effects of multiple environmental cues on single or multiple traits, giving rise to a highly flexible modelling framework.

The general framework for representing phenotypic plasticity in stage-structured population models
Variation caused by phenotypic plasticity is often expressed in terms of a reaction norm, a function that describes how an individual's environmental experience alters the phenotype they express (Nylin & Gotthard, 1998) (e.g. food consumed as a juvenile predicts adult body mass). We assume that organisms with similar experiences of their environment express the same phenotype and to represent this we create multiple linked copies of the Gurney et al. framework (as in Figure 2b). Each copy corresponds to a unique set of environmental conditions and individuals move through this structure on a path determined by their current and historical experience of the environment. This allows us to track cohorts of individuals that share the same environmental history and so represent the effects of phenotypic plasticity on populations.
Consider a stage-structured population with n lifestages where phenotypic plasticity is expressed according to d reaction norms r 1 ( ) , . . . , r d ( ) in response to z environmental cues (t) = 1 (t) , . . . , z (t) . For computational tractability, we discretise each environmental cue, j (t), into m j subintervals and denote by jp the midpoint of the p th subinterval of the discretisation of j (t). We define an environmental class to be a vector of length z with entries that consist of one midpoint from each discretised environmental cue, that is, where l j ∈ 1, . . . , m j . We define g: ℝ z → Ω such that if j (t) takes values within the l th j subinterval of j (t). The function g defines a mapping of onto the discretisation of . The number of environmental classes is given by and we denote the set of all such vectors (environmental classes) by Ω . We assign an ordering to the m elements of the set of environmental classes Ω (see Supplementary Information 6 for an example ordering) and let k denote the k th element of the ordered set, i.e. the k th environmental class. Thus, g ( ) = k if g ( ) is the k th element of the ordered set Ω . Environmental classes define cohorts of individuals that have experienced a shared environmental history. We assume that an individual's current and historic experience of the environment completely determines the phenotype an individual acquires when it matures from life stage i to i + 1 or is born into stage life stage 1 at time t. This permits the effects of environmental variation to be deferred to future developmental stages or generations. As we discretise the environmental cues this means there are a discrete number of phenotypes that can arise in the framework. The traits that these phenotypes express are calculated using the reaction norms according to, r k (g ( )). This process pre-defines both the traits individuals express and the range of environmental conditions that give rise to those individuals.

Incorporating phenotypic plasticity into stagestructured models
Denote the number of individuals in life-stage i and environmental class j at time t by N i,j (t). Denote by R i,j (t) the rate of recruitment of individuals into life-stage i and environmental class j, M i,j (t) the rate of maturation out of life-stage i and environmental class j, and D i,j (t) the death rate in life-stage i and environmental class j. The population is described by the system of equations for i ∈ 1, . . . , n, and j ∈ 1, . . . , m. The death rate is is the mortality rate of individuals in life-stage i and environmental class j. The recruitment term in Equation 1 when i = 1 is given by for j = 1, . . . , m where w kj ( (t)) denotes the proportion of individuals from environmental class k that transition to environmental class j at time t and v,k (t) is the birth rate of individuals in life-stage v and environmental class k. The transition functions, w kj ( (t)), are the mechanism through which the environment acts to express phenotypic plasticity within the model. Equation 2 represents the birth of new individuals into the first life-stage and environmental class j by parents from across all environmental classes and life-stages. The birth term, v,k (t) N v,k (t) describes the number of new individuals produced by parents in lifestage v and environmental class k, and is summed across all life-stages and environmental-classes to account for all new individuals entering the population. The transition functions w kj ( (t)) then determine the proportion of the new births that are assigned to environmental class j dependent on the environmental state.
The number of individuals recruited into life-stage i and environmental class j is equal to the number of individuals maturing out of life-stage i − 1 that are assigned to environmental class j. Hence, we have that for i = 2, . . . , n (1) Schematics of the ways phenotypic plasticity in stage-structured populations can be described by the new model framework. The population being considered in all cases is stage-structured with n life-stages. The number of individuals in life-stage i, expressing phenotype j is denoted N i,j . (a), The Gurney et al. ) framework for stage-structured populations that is used as a basis for the novel framework. This framework represents a continuous age structure by a discrete number of developmental classes e.g. eggs, larvae, pupae, and adults. (b), The most general form of the novel framework, where an individual's experience of the environmental cues in each developmental stage determines the phenotype it expresses as partitioned by the environmental classes. (c), The new framework adapted to represent developmental plasticity in life-stage 2. It is assumed that individuals experience an environmental cue in life-stage 1 that does not effect individuals in life-stage 1 but results in the expression of phenotypic plasticity in subsequent life-stages. This allows the reduction of the phenotypic structure in life-stage 1 to just a single class, The new framework adapted to represent a maternal effect in response to an environmental cue experienced by parents in life-stage n that manifests as phenotypic plasticity in life-stage 1 which is then assumed to have no effect on subsequent life-stages Further, the number of individuals maturing out of lifestage i and environmental class j is equal to the number of individuals recruited into life-stage i and environmental class j one developmental period ago that survived. We denote the duration of life-stage i for individuals in environmental class j by i,j . Thus, we have that Although in this formulation of the framework the stage duration i,j is kept constant an extension to variable stage duration (Ewing et al., 2016;) is also possible.
The exact form of the transition function, w kj ( (t)), is left unspecified as the way individuals transition from one environmental class and life-stage to the next is case specific. However, the choice of w kj ( (t)) is subject to the constraints 0 ≤ w kj ( (t)) ≤ 1, ∀j, k ∈ 1, . . . , m and ∑ m j = 1 w kj ( (t)) = 1. Although the transition functions are stage-independent, the environmental vector (t) is able to refer the state of each environmental cue independently and so can consider the sequence of past environments that an individual has encountered.

Application of the novel modelling framework to Nicholson's blowflies
To demonstrate the insights that can be gained from our framework we applied it to Nicholson's classical blowfly study (Nicholson, 1957), which aimed to describe how populations adjust in response to changes in their abiotic environment. In this study, the population dynamics of Lucilia cuprina (Wiedemann, 1830) were examined under different competitive conditions. In each culture, food was supplied separately to larvae and adults and both food supplies were replenished daily. Cultures were maintained for over two years and the number of adults and eggs present was recorded every two days. The results of Nicholson's study have been extensively discussed in theoretical ecology (Bakker, 1963;Glyzin, 2018;Gurney et al., 1983;May, 1986;Wood, 2010).
In Nicholson's experiment, blowflies experienced competition for food their larval and adult stages. Competition for food between adult blowflies reduces fecundity if individuals cannot acquire enough protein to mature all their eggs (Vogt et al., 1985). Larval competition for food reduces adult body size and the probability of survival through the pupal stage (Jannicke Moe et al., 2002). In blowflies, body size is linearly related to the number of ovarioles an adult has (Vogt et al., 1985) which determines the maximum number of eggs the adult can produce (Jannicke Moe et al., 2002), and so larval competition alters the maximum potential fecundity of adults. We regard the intensity of larval competition for food resources as an environmental cue, altering maximum potential adult fecundity and through pupal stage survival. It is important to note that maximum potential adult fecundity is distinct from observed adult fecundity, the former representing the maximum number of eggs an individual could produce under ideal environmental conditions and the latter representing the actual number of eggs an individual produces under the environmental conditions that individual experiences.
In Nicholson's culture (reproduced in Figure 3a here), the daily larval food supply was kept constant, but the amount of adult food supplied was reduced from an "unlimited" amount to a more limiting 1000 mg after around 600 days. The reduction of adult food resulted in an increased average adult population density, and the stabilisation of the previously regular population cycles. This is somewhat counter-intuitive, since a decrease in available resource substantially increased the average number of individuals and stabilised the previously regular oscillations-an example of the paradox of enrichment (Roy & Chattopadhyay, 2007).
Nicholson hypothesises that the population dynamics observed in the blowfly culture can be explained by phenotypic plasticity induced by larval competition. Adults in the period of unlimited adult food were observed to produce many eggs. When these eggs hatched into larvae, they experienced high levels of competition for larval food. This caused very few larvae to gain sufficient mass to pupate successfully, resulting in increased pupal mortality and low adult numbers in the next generation. When adult food was limited, an increase in adult competition resulted in fewer eggs being produced. The lower number of eggs resulted in fewer larvae and a larger amount of food being available per larva, subsequently reducing larval competition and juvenile mortality causing an increase in average adult population density. We evaluate evidence for Nicholson's hypothesis and heuristic arguments using the modelling framework derived here to represent phenotypic plasticity induced by resource competition in blowfly populations.

Model description
To formulate a model that represents phenotypic plasticity in blowfly populations we extend a previously derived mean-field model from Gurney et al. (1983), detailed in Supplementary Information 3, that considered only the instantaneous effects of adult competition on blowfly population dynamics. We introduce reaction norms relating through pupal-stage survival (Jannicke Moe et al., 2002) and maximum potential fecundity (maximum number of eggs an individual could produce in conditions of excess adult food) (Webber, 1954), to the availability of larval food. As this is an example of developmental plasticity the model takes the form of Figure 2c, for examples of models where we consider maternal effects or multiple environmental cues with a cumulative effect over multiple stages see Supplementary Information 5.
We assume eggs are laid into a single egg class, within which all individuals express the same phenotype (i.e. we assume no maternal effects). After a fixed developmental period, eggs hatch into a single larval class where again all individuals express the same phenotype. When a larva matures into a pupa, the amount of food that it obtained in the larval stage is determined by dividing the total food provided over the larval period by the number of individuals present in the culture over that period assuming scramble competition. The food obtained by an individual in the larval period is subsequently used to determine the traits that the individual expresses as a pupa and as an adult.
We abbreviate the previously introduced N i,j notation dropping the i and replacing the N by a more descriptive letter reflecting the life-stage, for example, L is used for larvae and A for adults. Similarly, as adults are the only explicitly modelled life-stage that expresses phenotypic plasticity (pupae are implicitly modelled due to a lack of density dependence ) the j subscript is dropped completely for terms relating to larvae. As we only consider a single environmental cue (larval food) the set of environmental classes, Ω , consists only of the midpoints of (t).
Denote by L (t) the number of larvae at time t and by A j (t) the number of adults at time t in environmental class j. Associated with each environmental class are the maximum fecundity of adults q j , and survival through the pupal and juvenile stages S J j . Recruitment into the larval stage is denoted R L (t) and recruitment of larvae to adults in environmental class j is denoted R A j (t). The environmental classes are parametrised by discretising an adapted reaction norm for through pupal-stage survival (Jannicke Moe et al., 2002) and a reaction norm for maximum adult fecundity is approximated from various sources (Webber, 1954). A detailed discussion of how the model is parametrised is detailed in Supplementary Information 7. As a proxy for the environmental cue, total protein obtained per larvae over the course of the larval period, we use the average protein available per larvae per day over the course of the larval period. We assume larvae divide the available food equally allowing the cue to be expressed as follows: (4)

F I G U R E 3
Simulation of the Nicholson blowfly culture data using the novel framework to represent phenotypic plasticity. In the culture adults blowflies were given unlimited food for 610 days, represented by K A = 2000 mg, which converts to 1800 mg of food supplied. After day 610 the amount of adult food supplied, K A , was then reduced to K A = 1200 mg, which converts to 1000mg of food supplied daily.

Time (Days) Pupal survival
where K L is the amount of larval food supplied daily, and L is the duration of the larval stage. This is converted into a derivative for ease of computation further detail of which is provided in Supplementary  Information 8. The model takes the form Recruitment terms are given by for j ∈ 1, . . . , m where I t − E is an inoculation term that begins the dynamics (Kot, 2001), and represents the introduction of larvae into the system at t = 0, and q j and S J j are determined by the reaction norms.
We assume that adults compete equally for the total available food regardless of phenotype, and so the instantaneous effects of adult competition are represented by e − A Tot (t − E )∕KA where A Tot = ∑ m j = 1 A j indicates competition across all phenotypes. The transition function w j ( (t)), determines the fraction of individuals entering environmental class A j (t) at time t, and is defined for j ∈ 1, . . . , m. This choice of w j ( (t)) restricts recruitment of individuals into a single environmental class based on that individual's experience of previous larval competition and indicates that maximum adult fecundity is uniquely determined by past experience of larval competition. This restriction is appropriate, as due to the assumption that all food is split equally, larvae being recruited at time t will have identical experiences of larval competition over the duration of the larval period, and so will express the same traits. We further assume that this developmental plasticity is irreversible. Although this choice of w j ( (t)) precludes microenvironmental variation this could be incorporated through a different choice of transition function.
The system is initialised with 9500 larvae at t = 0 with history for t ≤ 0 given by L (t) = 9500, (t) = K L ∕9500, A j (t) = 0∀j ∈ {1, . . . , m}. The stability analysis for this model is detailed in Supplementary Information 9. The model was simulated in R (R Core Team, 2020) using the package PBSddesolve (Couture-Beil et al., 2019).

Simulating population dynamics in the Nicholson blowfly culture under experimental conditions
The model simulated under replica experimental conditions shows good qualitative and quantitative resemblance to experimental data from Nicholson's blowfly culture (Figure 3a), capturing the culture's dynamical behaviour before and after food limitation. Initially, when the adult food supplied was unlimited, the model predicts the regular population cycles observed in Nicholson's data. After the food was restricted to 1000 mg at 610 days, the oscillations dampen, and the average population density substantially increases capturing the change in dynamical behaviour observed in Nicholson's culture. Although the population density of the simulated blowflies matches the experimental data the amplitude of the oscillations does not match. This mismatch can be explained by the high sensitivity of the model to food supply, a discussion of which, accompanied by food supplies which correctly predict both amplitudes, is provided in Supplementary Information 4. Initially, the population exhibits temporal cycles in the dominant phenotypes (Figure 3b,c). In the time periods where no new adults are being recruited, the phenotypic composition of the pupal and adult population does not change, resulting in the flat regions of Figure 3b,c. After food restriction, the range of phenotypes expressed within the population is greatly reduced. Pupae and adults in this period belong to a group of closely related environmental classes of individuals with relatively low trait values. As there is no difference between the distributions of maximum fecundity and through pupal stage survival, we hereafter only discuss fecundity.
As our model extends a previously derived non-plastic blowfly population model by Gurney et al. it is natural to question whether the population dynamics observed in Figure 3a can be attributed to the non-plastic population model. To test this, we simulate the non-plastic blowfly model, the formulation of which is provided in Supplementary Information 3. The non-plastic model overestimates the average adult density in both food conditions, predicts a decrease in adult density when w j ( (t)) = 1, if g ( (t)) = j 0, otherwise , resource availability decreases, and maintains the same population dynamics before and after the resource change as can be observed in Figure 4 (see Supplementary  Information 3 for further details).

Understanding the wider effects of phenotypic plasticity and population dynamic interactions
To explore how robust the population dynamics observed in the blowfly system are to conditions beyond those in Nicholson's experiment, we simulate the population trajectories for a wide range of possible combinations of adult and larval food supplies. For each food supply we record the average adult population density (Figure 5a), the average potential fecundity (maximum number of eggs an individual could produce in conditions of excess adult food, Figure 5b), the average observed fecundity (the number of eggs an individual actually produces in the context of competitive pressures within the population, Figure 5c), and the difference between the average potential and observed fecundities (Figure 5d). The model predicts that the blowfly population exhibits one of three dynamical behaviours. In the leftmost region of Figure 5a-d, increases in larval food supply do not change the abundance or fecundity (potential or observed) of adults. This suggests that in this region the population is limited most by the amount of adult food available. The population in this region consists of a small number of phenotypes with similar trait values (c.f. the small range of colours in the lines representing the abundance of individuals in each environmental class in Figure 5e and point (e) in Figure ). In the rightmost region of Figure 5a-d, increases in adult food supply do not change the abundance or fecundity of adults, suggesting that the availability of larval food is a limiting factor. The phenotypes expressed within the population are more diverse and the population's phenotypic composition changes during the course of a population cycle (c.f. the larger range of coloured lines in Figure 5g and point (g)). In the central region (the dark segment in Figure 5a), increases in either adult or larval food supply change the abundance and fecundity of adults. The adult population exhibits dampened oscillations and a small number of phenotypes with low trait values (Figure 5f and point (f)). This suggests that in this region the population is limited by the availability of both larval food and adult food. We conclude that the balance of resource availability between adult and larval blowflies governs the dynamical behaviour of the blowfly population. The population dynamics we observe are therefore characterised by the interaction between the two sources of density dependence: the instantaneous effects of adult competition and the delayed effects of larval competition through developmental plasticity.
Nicholson observed that when a culture initially supplied with 50g of larval food was supplied with 1g of adult food that 'the oscillation [of the blowfly population] was comparatively slight and had lost almost all evidence of periodicity, whereas any appreciable departure from the rate of one gram of ground liver per day in either direction resulted in the increase in the amplitude of oscillation'. The model predicts that when a population with a low adult food supply is supplied with increasingly more adult food that there is a sharp rise and then fall in average adult density as observed in Figure 5a. Similarly, when a population with a relatively low larval food supply is provided with increasingly more larval food, we observe a similar sharp rise and fall in average adult density. The behaviour Nicholson describes is precisely the behaviour that our model predicts, demonstrating that phenotypic plasticity is a mechanism by which the paradox of enrichment can be reconciled. The predictions the model makes about the link between traits expressed by individuals and population responses are somewhat counter-intuitive and would be difficult to anticipate from reaction norms alone. From consideration of only reaction norms it would be F I G U R E 4 Simulations of a previously derived non-plastic blowfly model from Gurney et al., (1983) under experimental conditions with adult blowflies initially supplied with K A = 2000 mg, which after day 610 is reduced to reduced to K A = 1200 mg. The solid black line indicates the total number of adults, while the dashed black line is the original data from Nicholson's culture reasonable to predict that individuals are most reproductively successful when potential fecundity is high. However, in Figure 5b-d we see that the food conditions that produce individuals with the highest average potential fecundity are also those that prevent this from being exploited and are associated with low average observed fecundity and consequently low reproductive success. By comparison, conditions that produce individuals with low average potential fecundity allow those individuals to achieve this potential, meaning that individuals in these conditions are on average more reproductively successful despite their lower trait value. Therefore, the relative contribution of high and low trait-valued individuals to the intensity of future larvae competition changes dynamically according to the environmental conditions the population is subject to. This demonstrates that the seemingly reasonable assumption we made from the reaction norm alone, that high trait value corresponds to high reproductive success, does not hold (Reed et al., 2013). We predict that individuals with traits that are indicative of high individual performance, such as average potential fecundity and average observed fecundity, arise from environmental conditions where the population is least abundant and most unstable. When average observed fecundity is highest (rightmost regions of Figure 5a-d) population density is lowest. Conversely, when the population density is highest, and the oscillations are damped (central regions of Figure 5a-d), the average observed fecundity is low, and the average potential fecundity is at a minimum. This shows that over most food conditions average potential fecundity and average observed fecundity are poor predictors of individual and population success. Even when average trait value is a good predictor of observed fecundity (rightmost region of Figure 5a-d) the population's dynamics are regulated by phenotypic plasticity and so would be misrepresented by an approach that uses averaged trait values. Conversely, when average trait value is a bad predictor of fecundity (leftmost regions of Figure 5a-d) an averaged trait approach correctly predicts the population dynamics (as can be observed by comparing Figure 4 to Figure 5e).
The food amounts supplied in the simulations shown in Figure 5e-g were selected to produce populations with the same average potential adult fecundity. Despite each population sharing the same average trait value, each population exhibits distinct dynamics and trait distributions which would be overlooked by an approach using averaged trait values. Only by accounting for trait variation between individuals arising from the cumulative effect of each individual's experience is it possible to capture the population-level effects of these three scenarios. This highlights the need to consider the individual and

DI SC US SION
We demonstrate that the interaction between phenotypic plasticity in individuals and population-level effects can be a source of rich population dynamical phenomena. The disconnect between individual and population performance demonstrated in the example of Nicholson's blowflies, although certainly not universal, provides a mechanistic explanation of how pressures that are maladaptive for individuals can be beneficial for populations and vice versa (Edelaar & Bolnick, 2019;Louthan et al., 2013;Reed et al., 2013;Weiner et al., 2017). By representing the mechanisms by which individual variation and population-level processes interact, we generate insights into how populations adapt to changing environments, which is crucial for understanding phenomena such as ecological tipping points (Dakos et al., 2019). Further, our findings support numerous previous studies proposing that failure to represent the effects of individual variation on populations is more consequential than simply mis-estimating demography (Bolnick et al., 2011;Lloyd-Smith et al., 2005;Sgrò et al., 2016;Violle et al., 2012), as we demonstrate that patterns in individual variation can drive complete changes in the dynamical behaviour of the system being considered. This is corroborated by observational studies where it has been found that the response of populations to interventions was influenced by individual variation (Cameron & Benton, 2004;Cameron et al., 2013). The framework is broadly applicable to systems where interaction between population dynamics and trait is important in determining the outcome of a process of interest. For example, when considering insect vectors of diseases or crop pests (e.g. mosquitoes or locusts) both abundance and trait to interact to determine the health or economic risk posed (Chandrasegaran et al., 2020;Sword et al., 2010).
Developing the model for Nicholson's blowflies was considerably simplified as we consider a well-studied model organism under controlled laboratory conditions. Although reaction norms are widely available across a broad range of taxa, outside of laboratory settings additional sources of environmental variation require the inclusion of reaction norms of higher dimension (i.e. a reaction norm considering the effect of temperature and con-specific density on development rate). However, our framework is designed to represent complex systems and so this should not pose an obstacle to implementation. For species where particular environment-trait relationships are not fully quantified or are missing entirely, due diligence must be observed in determining how sensitive the dynamics are to these uncertainties. For example, reaction norms are often most uncertain at environmental extremes (Brady et al., 2013) and so this uncertainty would need careful consideration when using our framework to predict dynamics at population range limits.
Although we demonstrate that individual variation can change and be changed by population processes, we do not predict when trait variation alters the outcome of these processes. In invasion biology, metrics derived from reaction norms are often used to predict the competitive viability of native and invasive species (Richards et al., 2006). Although the approach of using reaction norms directly accurately predicts the success of some invasive species (Knop & Reusser, 2012;Luo et al., 2019), it fails to explain the success of others (Muth & Pigliucci, 2007). This inconsistency limits the usefulness of reaction norms as a general predictor of a species invasiveness (Hulme, 2008;Palacio-López and Gianoli, 2011).
Here, we demonstrate that if one directly compares reaction norms without also considering a greater ecological context, one may arrive at erroneous conclusions. Therefore, it is important to determine more generally when reaction norms alone are sufficient to predict population dynamics and in doing so reconcile the role of phenotypic plasticity in biological invasions and population biology.