Integrating theoretical and empirical approaches for a robust understanding of endocrine flexibility

To maintain homeostasis, hormones mediate interactions between the external and internal environments. As examples, hormone levels can change in response to fluctuating environmental conditions, social interactions, and disease exposure, and these changes in hormone levels then stimulate modifications of behavior and physiology. Individuals and species vary in their endocrine responses to external cues in important ways that we can describe and quantify. We consider endocrine flexibility to be a sub-type of phenotypic flexibility, which is broadly defined as reversible variation in trait expression within an individual's lifespan (Hau et al., 2016). We use the term flexibility, rather than plasticity, to acknowledge its reversibility and ability to occur throughout an individual's lifetime, in contrast to phenotypic plasticity (see Glossary), which is sometimes defined as an irreversible change in phenotype induced during development (Pigliucci, 2005; Hau et al., 2016). Our broad definition of endocrine flexibility includes endocrine scope and speed (see Glossary; Fig. 1A; Taff and Vitousek, 2016), but also includes flexibility in baseline and induced hormone concentrations (Fig. 1B; Guindre-Parker, 2020). Endocrine scope and speed are measures of flexibility in the increases (and decreases) in hormone levels from baseline to induced levels. Individuals may also exhibit endocrine flexibility in either their baseline or induced hormone levels across environmental contexts (Guindre-Parker, 2020). Endocrine flexibility might be quantified by calculating endocrine scope and speed or by using reaction norms (see Glossary; Guindre-Parker, 2020; Malkoc et al., 2021). Using the reaction norm approach, endocrine flexibility is derived from the intercept and slope of hormone concentrations measured across contexts. The correlation between the intercept and the slope can also be quantified as a measure of the independence of these two metrics (e.g. Lendvai et al., 2014).

Endocrine flexibility may impact the ability of individuals, populations and species to respond to environmental change, fueling interest in its quantification. To date, empirical research has primarily documented that endocrine flexibility exists and has started to test its function and adaptive value. Integrating mathematical modeling with future empirical studies can help us understand the mechanisms, evolution and adaptive value of endocrine flexibility. We promote the integration of mathematical modeling with studies of endocrine flexibility, echoing other recent calls for increased interdisciplinary efforts that merge theoretical modeling with biological research (White et al., 2021). We outline additional avenues for empirical research, potential modeling approaches, requirements for utilizing models and questions for which different approaches are best suited. We recognize that not all researchers enjoy working with mathematical models and encourage empirical endocrinologists to examine published models of their system and to contact the authors of the models that best meet their accepted assumptions (Box 1). Collaborative approaches to asking and answering questions about endocrine flexibility with mathematical models will yield more robust and out-of-the-box studies than any narrow-field attempt.

Endocrine flexibility can be produced through a variety of mechanisms that allow some individuals to manufacture and clear hormones faster, or mount more robust responses than other individuals do (Fig. 1). For instance, endocrine flexibility may depend on upstream releasing or inhibiting hormones, such as those secreted by the hypothalamus. In the hypothalamic–pituitary–gonadal (HPG) axis, gonadotropin-releasing hormone (GnRH) from the hypothalamus stimulates secretion of luteinizing hormone (LH) and follicle-stimulating hormone (FSH), which stimulate reproductive activity and sex steroid hormone secretion. Individuals may be flexible in their hormone secretion along this axis, thereby affecting downstream responses (Aboul-Ela et al., 1983; Bluhm et al., 1991). Testosterone secretion in response to HPG activity can also vary seasonally, and there is evidence that individuals have repeatable differences in testosterone secretion (Jawor et al., 2006). Similarly, in the hypothalamic–pituitary–adrenal (HPA) axis, variation in stress reactivity may produce flexibility in secretion of corticotropin-releasing hormone (CRH) among individuals (Anisman et al., 1998). As a result, the speed and secretion of downstream hormones, including adrenocorticotropic hormone (ACTH) and glucocorticoids, can vary. Steroid hormones typically are bound to proteins while in circulation and may be released to diffuse into tissues and bind to intracellular receptors. Corticosteroid binding globulin (CBG), which binds glucocorticoids, may play a role in endocrine flexibility by modulating glucocorticoid levels (Breuner et al., 2013). However, CBG also binds other steroids with varying affinities, and CBG-bound glucocorticoids or CBG itself may be biologically active in the stress response (Schoech et al., 2013), complicating the role of CBG in endocrine flexibility.

Variation in the affinity and capacity of receptors on target cells can alter the strength of the response to hormones and influence downstream endocrine flexibility. Differences in steroid production may also arise through interactions of steroid receptors with heat shock proteins and chaperones that alter the availability of the receptor for hormone binding (Walker et al., 2017). A given receptor may have several isoforms (see Glossary), and these may also contribute to endocrine flexibility. For instance, isoforms of the GnRH receptor differ in expression based on life history stage (Joseph et al., 2009; Ciani et al., 2020) and season (Ciani et al., 2020).

Hormones can have pleiotropic effects and may influence multiple phenotypic traits, including physiology, morphology and behavior (Dantzer and Swanson, 2017; Mauro and Ghalambor, 2020). Hormones may then alter endocrine flexibility through their influence on other molecules that control cellular activity. Relatedly, hormones act as part of physiological regulatory networks, rather than linear pathways with a single response (Cohen et al., 2012; Mauro and Ghalambor, 2020). Regulatory molecules, such as hormones and receptors, are linked through their regulatory relationships (Cohen et al., 2012). Physiological regulatory network state reflects the concentrations of regulatory molecules within a context. The state of a physiological regulatory network can shift, altering multiple regulatory molecules to cope with internal or external changes and maintain homeostasis (Cohen et al., 2012; Mauro and Ghalambor, 2020). Variation in physiological regulatory network state can occur within and among individuals, and across time and contexts (Cohen et al., 2012). Flexibility in the expression or interactions of individual regulatory components as well as in physiological regulatory network state (Di Poi et al., 2016) may be important drivers of endocrine flexibility.

Endocrine flexibility may also arise through individual variation in perception and evaluation of stressors. In some instances, previous exposure can alter the perception of a stressor. In wood frog (Lithobates sylvaticus) tadpoles, exposure to predator cues caused increased glucocorticoid levels (Bennett et al., 2016). When tadpoles were exposed to predator cues continuously for 3 weeks, however, their glucocorticoid levels after cue exposure were lower than those in tadpoles raised without predator cues (Bennett et al., 2016). Thus, habituation to a stressor may influence perception and the scope of the endocrine stress response. Social rank and social instability can also affect stress responses. In bison (Bison bison), dominant bulls have higher glucocorticoid levels (Mooring et al., 2006), but in olive baboons (Papio anubis), higher glucocorticoid levels are found in subordinates (Virgin and Sapolsky, 1997). Regarding social stability, rhesus macaques (Macaca mulatta) in unstable social groups had lower glucocorticoid levels than individuals in stable social groups, which may be due to altered negative feedback sensitivity in macaques that are chronically stressed from social instability (Capitanio et al., 1998). Such alterations may dampen the ability to elevate glucocorticoids in response to a stressor.

There are a variety of ways in which individuals may display endocrine flexibility (empirical recommendations: Box 2; modeling recommendations: Box 3). First, individuals may exhibit flexibility in baseline hormone levels (see Glossary) under different environmental conditions (Guindre-Parker, 2020; Lendvai et al., 2014). As an example, food-restricted house sparrows (Passer domesticus), increase baseline corticosterone levels as their body mass decreases but differ in glucocorticoid responses to the restriction (Lendvai et al., 2014). Similarly, female North American red squirrels (Tamiasciurus hudsonicus) increase baseline glucocorticoid levels in higher population densities but differ in the degree to which glucocorticoid levels change in response to density (Guindre-Parker et al., 2019). A study of male chimpanzees (Pan troglodytes) that measured cortisol in urine samples (Box 2) found that cortisol levels were highest in the morning and declined across the day; however, initial cortisol levels and the magnitude of decline throughout the day varied consistently among males (Sonnweber et al., 2018). Second, individuals may exhibit flexibility in induced hormone levels (Glossary) in response to environmental (Lendvai et al., 2015) or pharmacological challenges. For example, older house sparrows elevate corticosterone less in response to capture and restraint than younger birds do, but individuals change their endocrine flexibility differently with age (Lendvai et al., 2015). Finally, individuals may differ in endocrine flexibility in the way their hormone levels change in response to environmental or pharmacological challenges (Ambardar and Grindstaff, 2017). One example comes from work on testosterone production in response to a GnRH challenge in Eastern bluebirds (Sialia sialis). Male and female bluebirds displayed consistent individual variation in testosterone levels but did not differ in potential endocrine scope in response to GnRH challenge (Ambardar and Grindstaff, 2017). This study highlights the use of releasing hormone challenges to quantify potential endocrine flexibility (Box 2; Taff and Vitousek, 2016). Two studies of different fish species have not detected consistent individual variation in endocrine flexibility. In three-spined stickleback (Gasterosteus aculeatus), there was a significant effect of individual identity on baseline cortisol levels, but individuals did not differ in their cortisol responses to a simulated predator (Box 2; Fürtbauer et al., 2015). Similarly, individual Trinidadian guppies (Poecilia reticulata), differ significantly in baseline cortisol levels, but do not display consistent differences in induced cortisol levels during habituation to a stressful environment, although there was some evidence that individuals with higher average cortisol levels also had higher slopes of cortisol levels across repeated exposures (i.e. positive correlation between the intercept and slope of reaction norms; Houslay et al., 2019).

For selection to act on endocrine flexibility, there must be heritable variation in endocrine flexibility traits. Little work has assessed the heritability of endocrine flexibility, but existing analyses suggest the potential for endocrine flexibility to evolve. Baseline and restraint stress-induced corticosterone levels have low heritability in nestling tree swallows (Tachycineta bicolor), whereas the difference in corticosterone levels between restraint stress-induced and baseline levels was moderately heritable (Stedman et al., 2017).

Considering the proximate and ultimate causes of endocrine flexibility, there are limits to endocrine flexibility that might reflect physiological limits or be the products of trade-offs that disfavor very rapid responses (Rich and Romero, 2005; Lipowska et al., 2020). Importantly, selection may shape physiological limits of hormone production (Clotfelter et al., 2004; Øverli et al., 2007). For instance, glucocorticoid levels are repeatable in various taxa (Rensel and Schoech, 2011; Narayan et al., 2013; Vuarin et al., 2019), and may interact with fitness (Comendant et al., 2003; Vitousek et al., 2018b), making them potential targets of natural selection. Thus, it is likely that both mechanistic and evolutionary forces play a role in the causes of variation in endocrine flexibility. Recent reviews have addressed the environmental conditions that might favor endocrine flexibility and some mechanisms that could contribute to variation in endocrine flexibility (Guindre-Parker, 2020; Hau et al., 2016; Malkoc et al., 2021; Taff and Vitousek, 2016; Wada and Sewall, 2014).

To determine if endocrine flexibility is limited by mechanistic constraints, releasing hormone challenges and pharmacological manipulations can be used to probe maximal potential endocrine scope or speed in relation to realized endocrine scope and speed (Taff and Vitousek, 2016). Hormone challenges and pharmacological manipulations can also be administered within different environmental and social contexts to test for individual variation in maximal potential endocrine flexibility (e.g. Ambardar and Grindstaff, 2017) and can be compared with non-challenge hormone levels in the same environmental and social contexts to determine if maximal and realized endocrine flexibility are the same or if they might differ owing to fitness trade-offs.

At the ultimate level, to determine if fitness trade-offs shape realized endocrine flexibility, investigators could test relationships between endocrine flexibility and reproductive success or survival. Realized endocrine flexibility may be less than maximal endocrine flexibility because selection acts against very rapid responses or maximal flexibility (Taff and Vitousek, 2016). Realized endocrine flexibility may be constrained by costs associated with mismatches between an optimal endocrine phenotype and current environmental conditions, particularly when endocrine speed (see Glossary) is slower than the rate of environmental change. In these circumstances, optimal endocrine responses and flexibility are shaped by current and likely future environmental conditions, and overall endocrine scope is reduced (Gabriel et al., 2005; Luttbeg et al., 2021).

At the proximate level, the contributions of different mechanisms to endocrine flexibility can be investigated using theoretical models of dynamic systems, such as dynamical models and virtual screening techniques (Fig. 2, Supplementary Materials and Methods). At the ultimate level, how endocrine flexibility is expected to optimally evolve can be predicted using state-dependent dynamic models (Taborsky et al., 2020) and these can be linked to evolutionary dynamics using adaptive dynamics or genetic algorithms (Fig. 2, Supplementary Materials and Methods). Mechanistic limitations can be incorporated using evolutionary simulations (Taborsky et al., 2020). Predictions about the direction and speed of evolution are best made using quantitative genetic models (see Glossary) or evolutionary simulations.

One way researchers could use modeling to better understand how biochemical processes influence endocrine flexibility is through construction and manipulation of differential equation-based models (i.e. dynamical models of biochemical pathways; see Glossary). The type of model selected by the researcher will depend on the assumptions made and the data available (dynamical models of the HPA axis, reviewed in Stanojević et al., 2018). Dynamical models integrate nicely with empirical work as they are constructed using empirically collected data on biochemical reactions. Additionally, the assumptions about the physiological and biochemical processes underlying an endocrine response used to construct these models and predictions made by the models can be tested with empirical studies, such as pharmacological manipulations. In this way, researchers can evaluate otherwise difficult to test mechanistic hypotheses.

Dynamical models of the HPA axis have been developed to examine potential causes of endocrine flexibility. Examples include models examining interactions between the HPA axis and the immune system (Malek et al., 2015), the regulatory network controlling glucocorticoid synthesis (Spiga et al., 2017) and regulatory mechanisms of the HPA axis under acute and chronic stress (Marković et al., 2011). The structure of these models varies widely. In one model of the HPA axis, equations were generated for the change in concentrations of CRH, ACTH and cortisol, and only negative feedback from cortisol on CRH and ACTH was included (Vinther et al., 2011), while other models incorporate the steroid precursor cholesterol, as well as aldosterone (Marković et al., 2016). One can see, then, how the predictions of – and relative support for – these models by empirical work using pharmacological or environmental manipulations can yield insight into the relative importance of different physiological mechanisms of endocrine flexibility. Multiple mathematical models also exist for the HPG axis (reviewed in Clément, 2016; Clément et al., 2020).

Virtual screening techniques (see Glossary) are used in drug discovery and can further our mechanistic understanding of endocrine flexibility in several ways (Gohlke and Klebe, 2002). Models using these techniques are divided into structure- and ligand-based methods. Structure-based methods (e.g. molecular docking) can be employed when information about the 3D structure of the receptor is known; if the 3D structure is not known, then ligand-based methods [e.g. pharmacophore modeling and quantitative structure activity relationship (QSAR)] can be used (Halperin et al., 2002; Hamzeh-Mivehroud et al., 2016). In essence, these models examine the interaction between small molecules and receptors at the subatomic level to determine the behavior of those molecules at binding sites.

Within the context of endocrine flexibility, virtual screening techniques might be used to explore how individual variation in binding site activity contributes to endocrine flexibility, investigate how environmental factors or other biochemical molecules interact with other endocrine systems to contribute to endocrine flexibility, or study compounds to experimentally manipulate endocrine flexibility. Like dynamical models, these models are built on assumptions that reflect hypothesized physiological and biochemical relationships. Empirical studies, especially pharmacological challenges, can effectively test the assumptions and predictions of these models to evaluate hypotheses. Theoretical exploration into the proximate mechanisms of endocrine flexibility can help to focus the effort and resources of future empirical work.

Because endocrine responses help to maintain homeostasis, it would seem adaptive for endocrine responses to be as flexible as mechanistic restraints allow. However, trade-offs might favor slower responses or more variation among individuals in response rates. One potential trade-off is balancing the benefits of quickly responding to environmental changes against the costs of mistakenly responding to false information. Individuals receive cues that contain information about the current environmental state, such as the presence of predators or lack of available food. Cues are rarely perfect, and individuals are at risk of perceiving false positives (e.g. cue of a stressor that is absent) or false negatives (e.g. failure to perceive the presence of a stressor). Given imperfect information, there is an optimal amount of information needed to trigger an individual to up- or downregulate hormone levels depending on the costs and benefits of quicker responses versus false responses. If stressor severity is best described along a continuum (e.g. the probability a predator is present), then optimal hormone regulation can be found using Bayesian updating and optimality models (see Glossary) that assume natural selection will ultimately favor phenotypes that produce the highest average reproductive success of individuals (Luttbeg and Trussell, 2013; Zimmer et al., 2022). If environmental states can be simplified to a binary presence or absence of stressors, then signal detection theory can be used to find the optimal threshold for triggering the up- or downregulation of hormones (Getty, 1985).

Often optimal responses depend on an individual's state, such as their age, sex, size, hunger, or estimate of the state of their world. For example, how much predation risk an individual is willing to encounter to find food should depend on hunger levels (Houston and McNamara, 1993). If there are resource or opportunity costs to synthesizing or maintaining hormones, receptors, enzymes or binding proteins that affect the flexibility of hormonal responses, then there could be a trade-off where faster regulation has higher costs. If individuals differ in the resources available for building and maintaining endocrine systems, then this state variable would affect the optimal trade-off of maintaining endocrine flexibility. One of the consequences of state dependency is that current hormonal responses can affect future response possibilities. For example, endocrine responses that affect foraging effort affect how much energy is available later. Thus, the optimality of different responses depends not only on how responses affect current success, but also how they affect the future options for the individual and the expected success of those options. This can be a tricky question to solve. Fortunately, dynamic state variable modeling methods were developed starting with the key insight that state-dependent questions are more tractable if one works backwards (Mangel and Clark, 1989).

Optimality models, including state-dependent models, can easily be integrated with empirical work on endocrine flexibility by using experimentally derived values to parameterize the models (e.g. hormone concentration ranges, survival rates, fecundity) or by designing empirical studies that test model predictions. Studies that determine the relationship between realized and maximal endocrine flexibility and fitness (e.g. survival and reproduction) or modulate the environment in predictable or unpredictable ways (as described in the empirical recommendations) and quantify its effect on realized and maximal endocrine flexibility and fitness could be used to test model predictions.

Mathematical models can be used to predict and give insight into the evolution of endocrine flexibility. Quantitative genetic models have long been used by evolutionary biologists to make predictions about the direction and speed of evolution. A core assumption of these models is that traits are influenced by many genes with small additive effects, and thus the traits and their changes are continuous (Hill, 2010). They can be used to ask if endocrine responses are likely to evolve given available genetic variability and selection strengths. Dantzer and Swanson (2017) used a quantitative genetics approach to study the relative strength and occurrence of hormonal pleiotropy. They found that whether hormonal pleiotropy of two traits facilitated or constrained each other's evolution depended on the relative directions of selection on the traits and the directions of correlations between each trait and a shared hormone. Quantitative genetic models could be used to answer questions about endocrine flexibility such as how varying selection strength affects the speed of evolution of flexibility, how the type and strength of selection on flexibility affects correlated traits, and how opposing selection pressures between the sexes or among age classes might affect the evolution of endocrine flexibility. Because evolutionary questions are difficult to address empirically, unless a researcher works with an organism in a controlled environment with fast generation times, quantitative genetic models provide insights into the evolution of endocrine flexibility that may be unobtainable via empirical work. Should a researcher with a suitable model system wish to design an experiment to test evolutionary hypotheses about endocrine flexibility, however, quantitative genetic models can provide predictions to which empirical results can be compared.

Genetic algorithms (see Glossary) are another approach that can be used to predict how endocrine systems will evolve. In these models, individuals with higher fitness pass more alleles into the next generation. Over generations, the population evolves with alleles that confer greater fitness increasing in frequency. The purpose of the approach can differ based on the degree of mechanistic specificity. It can be viewed as an algorithm to find solutions to optimality problems that are too difficult to solve analytically with no claim that the observed process of evolution is likely to occur (Mitchell, 1996). However, it is possible to include more mechanistic details, including the location and interaction of genes, recombination and mutation rates, and then present the results as a prediction of the likely course of evolution. There are few examples of this modeling approach in the field of endocrinology. However, Bourg et al. (2019) used genetic algorithms to examine how endocrine systems should evolve and shape the trade-off between investing in two different traits that affect fitness. In their model, individuals acquire energy through a meal, then genes coding for hormones and receptors determine how that energy is allocated to two different traits or stored for future allocation. They find that the shape of the trade-off in investing in the two traits depends on mutations in the regulatory and coding regions of genes for hormones and their receptors, selection on the endocrine system, and the efficiency costs of storing energy. Genetic algorithm models would allow the incorporation of more mechanistic details of the genetic architecture underlying endocrine flexibility, obtained through empirical studies, to predict its likely evolution. Like quantitative genetic models, genetic algorithm models allow researchers to examine trade-offs that are difficult to assess empirically.

Adaptive dynamic models (see Glossary) are also used to find optimal trait combinations in populations inhabiting different environments. However, unlike genetic algorithm models, adaptive dynamic models assume that all individuals in a population use a set trait or proportion of traits, and then a mutation in the trait is introduced to assess if it can outcompete the current trait in the population. The expectation is that after repeated iterations of mutational invasion, ultimately a trait will be found that cannot be invaded, and that is the optimal trait. By ignoring variation within the population, this approach ignores the effects of frequency dependence on the success of traits but can be an effective modeling approach if that is not a major consideration. Quantitative genetic models, genetic algorithm models and adaptive dynamic models are all useful for understanding evolutionary questions of endocrine flexibility, and allow researchers to examine trade-offs ultimately underlying endocrine flexibility that would be difficult to test empirically. These approaches differ in the level of mechanistic detail that can be incorporated and the starting assumptions.

Although there are many theoretical options that could be used to investigate the causes of endocrine flexibility, only one study (to our knowledge) compares the predictions of two types of models. Taborsky et al. (2020) used two different types of models to understand why the pattern of low glucocorticoid levels in the absence of stressors, a quick rise to a peak in response to a stressor, and a decline back to baseline is ubiquitous across vertebrates (Breuner et al., 2008; Romero and Wingfield, 2016). Their first model, a dynamic state variable model in which they altered the autocorrelation between current and future environmental states, demonstrated that environmental predictability may be a primary factor generating the conserved shape of the glucocorticoid response to stress. Their second model, a genetic algorithm model, incorporated mechanistic constraints on how rapidly hormone levels could be upregulated in response to stressors with a single hormone clearance rate. This mechanistic model produced similar baseline and stress-induced hormone levels to a model without constraints on rates of hormonal regulation, but much slower clearance rates. Their models highlight that assumptions about the mechanisms of hormone regulation can exert large effects on the predicted dynamics of hormone regulation, and that the process of modeling leads to assumptions being clearly stated and thus subject to discussion. The next step might be to collect empirical data to determine which model is most accurate by manipulating environmental predictability or by determining whether there is a single clearance rate, and if so, what the consequences might be for how endocrine systems work.

To assess the consequences of variation in endocrine flexibility in response to environmental variability, investigators could manipulate environmental predictability and then measure endocrine flexibility as well as fitness. Alternatively, researchers could compare populations or species inhabiting environments that differ substantially in predictability. Even when circulating hormone levels are similar across populations that differ in environmental predictability, there may be differences in other components of HPA axis regulation such as CBG and receptor affinity and capacity that confer greater endocrine flexibility (Breuner et al., 2003). However, we would predict that environmental unpredictability would enhance endocrine flexibility when either the rate of environmental change is slow enough or hormone responses are fast enough to minimize mismatch between the environmental state and endocrine responses. Studies have not yet simultaneously quantified endocrine flexibility and fitness to assess how differences in components of flexibility might impact fitness.

Other approaches that would provide valuable insight into the consequences of variation in endocrine flexibility within species include selection line experiments, studies of pedigreed populations, cross-fostering experiments and common garden experiments. Artificial selection on performance or behavior may have correlated effects on endocrine flexibility that provide insight into the consequences of variation in endocrine flexibility. Bank voles (Myodes glareolus) selected for aerobic exercise metabolism while swimming, predatory responses or body mass maintenance on a poor-quality diet differed in endocrine scope and speed, but did not differ in baseline or stress-induced corticosterone levels (Lipowska et al., 2020). In particular, voles selected for physical activity and alertness on the swimming task or predatory assay had decreased endocrine scope, whereas animals selected for responses to prolonged dietary restriction had decreased endocrine speed. Moreover, voles in the dietary restriction line had reduced glucocorticoid production in response to restraint stress relative to the maximal potential response (determined by ACTH challenge; Lipowska et al., 2020). These results demonstrate how stressor type and duration may alter different components of endocrine flexibility. Rainbow trout (Oncorhynchus mykiss) lines have been developed that differ in growth on an herbivorous diet early in life. As in the vole study, trout with higher early growth rates on the plant-based diet have reduced cortisol release rates (Sadoul et al., 2015). Fitness can be difficult to quantify in the captive conditions required for artificial selection studies. Accordingly, analyses of pedigreed wild populations can provide valuable information linking variation in endocrine flexibility with fitness consequences. In tree swallows, glucocorticoid responses are heritable (Jenkins et al., 2014) and baseline and acute stress-induced levels are not genetically correlated (Stedman et al., 2017; but see Béziers et al., 2019). If baseline and induced hormone levels are not linked, then this would allow for independent evolution and more rapid responses to environmental change (Dingemanse and Wolf, 2013). Common garden experiments in which organisms from different habitats that differ in endocrine responses are tested in a shared environment could be used to determine if the previously observed differences were the results of microevolutionary differences or phenotypic plasticity (Bókony et al., 2021).

Individual variation in endocrine speed can have consequences for baseline and acute stress-induced glucocorticoid levels. Recent modeling work predicted that slower rates of glucocorticoid release lead individuals to have elevated baseline hormone levels, while slower rates of downregulation should lead individuals to reduce acute stress-induced glucocorticoid levels (Luttbeg et al., 2021). For individuals with slower endocrine speeds, there is a longer lag time for glucocorticoid levels to either increase or decrease after stressor exposure, and during those periods of time, hormone levels are not optimally regulated to maximize fitness. Therefore, individuals with slower upregulation rates might be expected to have higher baseline glucocorticoid levels to keep them closer to the appropriate stress-induced level needed during periods of acute stress. The approach to modeling this question started with a dynamic state variable model to determine how optimal hormone levels are affected by an individual's current glucocorticoid levels (a state variable), the probability of changes in the environmental state, and the specified rates of up- and downregulation. The existence of finite rates of endocrine regulation was made as an assumption, and the consequences were explored. Then, a simulation was used to view how individuals would respond given the model assumptions and the resulting solutions of the optimality model, with a focus on hormone regulation. One could also investigate the effects of endocrine flexibility on health, performance, allostatic load or reactive scope and test model predictions with an empirical study in which endocrine flexibility is environmentally or pharmacologically manipulated and the relevant consequences documented. For example, one could model how differences among individuals in endocrine flexibility impact the transition from reactive homeostasis to homeostatic overload (Romero et al., 2009) and then design an experiment to test model predictions. Furthermore, environmental variation and its effect on endocrine flexibility can also be incorporated using these modeling approaches, allowing researchers to generate testable predictions for empirical studies that manipulate environmental predictability and examine the resulting consequences for endocrine flexibility and fitness.

If one makes assumptions about individual endocrine responses or if those responses are predicted from optimality models, then population-level models such as integral projection models (IPMs; see Glossary), Euler–Lotka models (see Glossary) or individual-based simulation models can be used to make predictions about the population- or community-level consequences of variation and limitation in endocrine flexibility (Fig. 2, Supplementary Materials and Methods). IPMs can be applied to any stage-structured population – a population that has distinct life stages that vary in survival, reproductive, and growth rates – or size-structured population, and are built on regression models that link individual state data (e.g. endocrine responses) to vital rates, so empirical data linking these two measures are central to the use of this modeling approach (Easterling et al., 2000). The regression models explaining the relationship between individual state variables and vital rates provide opportunities for researchers to augment the validity of their model by incorporating abiotic and biotic covariates that might explain additional variation in population vital rates. Although IPMs primarily generate estimates of population dynamics, these population estimates can be expanded to explore the consequences of variation in individual state on emergent biological patterns such as range limits (Merow et al., 2014). For example, a researcher could conduct experimental or observational studies to document how individual or population variation in an organism's ability to exhibit endocrine flexibility relates to their probability of survival. Those data could then be used to construct an IPM to determine how population dynamics might differ between populations that exhibit variation in endocrine flexibility.

Euler–Lotka models find population growth rates using the Euler–Lotka equation, which is based on age-specific survival and per capita reproductive rates. With this approach, researchers can define equations that describe how survival and reproduction change with individual-level parameters (e.g. growth rate and age at maturity; Mangel and Stamps, 2001) and examine the effect that variation in individual parameters has on population growth rates. Like IPMs, however, a researcher interested in exploring the population-level effects of endocrine flexibility either must have empirical data to parameterize the relationship between endocrine flexibility and population vital rates, or they must provide robust theoretical support for the parameters they select. Regardless, integrating empirical research and Euler–Lotka models can provide insights into the population-level consequences of endocrine flexibility that cannot be obtained independently through empirical or theoretical work.

Individual-based simulation models (IBMs; see Glossary), also called agent-based simulation models, simulate populations and communities based on the actions and interactions of individuals that vary in state variables (Railsback and Grimm, 2019). Unlike other population and community modeling approaches that are based on initially generating population vital rates, IBMs are built from the bottom-up – that is, the behavior of the population-level parameters results from the interactions of individuals with their simulated environment. Preliminary data on the relationship between state variables and population vital rates are not necessary. For example, previous models have used IBMs to investigate how stressor exposure and the costs of coping with stressors can impact population size (Fefferman and Romero, 2013). An added benefit of employing IBMs to investigate the population- or community-level consequences of variation in physiological traits, such as endocrine flexibility, is that any number of individual state variables can be incorporated. So, a researcher may be able to simultaneously explore the effect that several aspects of endocrine flexibility have on population or community dynamics. IBMs can incorporate data from any experiment aimed at quantifying intraspecific variation in endocrine flexibility (e.g. selection line experiments, pedigreed populations, common garden) and generate testable predictions of population size. As a result, the integration of empirical data with IBMs can provide a robust understanding of the consequences of endocrine flexibility.

We have outlined empirical approaches that can further understanding of the causes and consequences of endocrine flexibility and advocate for future integration of empirical and theoretical techniques. Modeling can significantly advance our understanding of endocrine flexibility by filling in gaps where empirical studies are intractable, identifying gaps in our understanding of the ‘how’ and ‘why’ of endocrine flexibility, and generating additional predictions for empirical testing. Furthermore, although there are many appropriate approaches for modeling endocrine flexibility, model choice should ultimately depend on the purpose of the research.

The authors would like to acknowledge the helpful comments of two anonymous reviewers and the invitation from Charlotte Rutledge to contribute to the Special Issue.