Ecological Stoichiometry Theory (EST) explores the connections between an organism’s elemental composition and their life history traits, interactions with other organisms, and impacts on ecosystem functioning (Urabe et al., 1997; Sterner and Elser, 2002; Anderson et al., 2005; Frost et al., 2005). Using elements as a currency for metabolic and ecological processes of interest is particularly useful when connecting mechanisms between organisms and/or ecological scales, such as the relationship between animal egestion and fertilized primary production. The ability to quantify how macroelements, such as N, flow throughout an ecosystem is an informative tool for describing consumer-driven nutrient cycling (Elser and Urabe, 1999, 1999; Sterner and Elser, 2002; Grover, 2004; Moody et al., 2018) and how feedbacks between consumers and their environment impact ecological and biogeochemical processes (Sitters et al., 2017; Filipiak and Filipiak, 2020; Elser et al., 2022; Koeve et al., 2024).

EST can also be used to predict the fitness consequences of elemental imbalance between an organism’s demands and what is available in their environment (i.e., stoichiometric mismatch) (Sterner and Elser, 2002). These consequences can occur in different ways. First, element deficiency (i.e., nutrient limitation) can suppress organism performance (e.g. growth rate, reproduction, etc.) relative to more elementally balanced diets (Sterner et al., 1992). Second, consumers face costs when dealing with an overabundance of some elements in digesting, egesting, excreting, and/or storing unneeded nutrients, which has been shown to induce trade-offs in organismal performance (Boersma and Elser, 2006). The effects of stoichiometric imbalances on individual performance, population dynamics, and whole ecosystems can often be isolated experimentally, but it can be difficult to interpret observational data where there are often multiple, interacting limiting nutrients. Here, EST can explicitly describe the dynamic shifts in elemental limitations through the concept of TERs. TERs are the resource elemental ratio at which an organism’s growth shifts from limitation by one element to another (Frost et al., 2006), allowing for a combination of elemental ratios to describe a “stoichiometric niche” of resources (González et al., 2017). If a consumer cannot meet these thresholds, it must adapt or face fitness consequences. The strategies for adapting to, and quantitative consequences of, elemental imbalances both lend powerful and complementary modeling concepts when used in conjunction with both NG and DEB frameworks (Sperfeld et al., 2017).

DEB theory is a mechanistic representation of an individual organism’s metabolism that directly links resource quantity and quality to its life-history traits by connecting resource composition to physiological processes (Kooijman, 2000). DEB models track how assimilated resources are stored in organisms (i.e., reserves) and then mobilized to pay metabolic costs for maintenance (i.e., overhead costs of maintaining current biomass), growth, and reproduction. Reserves provide a temporal buffer for consumers to sustain their metabolic function during resource fluctuations, effectively modeling individual-level nutrient hysteresis (Nisbet et al., 2000; Sousa et al., 2010). DEB theory also models concepts described by EST and NG, such as nutritional constraints on growth and reproduction (Muller et al., 2001, 2004; Kuijper et al., 2004) and consequences of resource switching (Kooijman, 2000). DEB models have already been integrated into ABMs to scale individual-level energetics to population dynamics (Martin et al., 2013), model zoonotic infection control scenarios (Malishev and Civitello, 2019), develop predictions for population management under climate change conditions (Yang et al., 2022), and better understand animal life history trajectories (Goedegebuure et al., 2018). For those unfamiliar with the intricacies of DEB, this section aims to present relevant assumptions, mechanisms, and modules available for DEB models that would be relevant to an ABM integrated with EST and NG concepts.

Two assumptions for DEB theory underlie its connection to biochemical composition of organisms: strong and weak homeostasis (Kooijman, 2000; Nisbet et al., 2000). Strong homeostasis dictates that metabolic pools within the organism do not change chemical composition (i.e., somatic structure, reserve pools, etc.). Weak homeostasis dictates that under constant resource density environments, the composition of the consumer as a whole remains static, implying that as resource quantity or quality varies, the reserve density (amount of reserve relative to the somatic structure, m E = M E M V ) drives the overall chemical composition of the consumer. For example, an individual facing extended nutrient-poor conditions will have depleted reserves, and thus its reserve density will be much lower than an organism fed abundant, high-quality food (whose reserves will be at capacity). The difference in composition between reserves and structure determines the range of whole-body composition of an organism (i.e., heterotrophic consumers have a stricter compositional range than primary producers that can separately uptake inorganic nutrients) (Sousa et al., 2010). These theoretical tenets of DEB provide a complementary approach to EST’s description of the elemental homeostatic flexibility in producers and consumers.

DEB models can also use multiple reserves to represent different generalized pools of biochemical compounds (e.g., C-rich and N-rich reserves) (Kooijman, 2000). To convert fluxes from reserves of different stoichiometries into biomass products of different stoichiometries (i.e., somatic tissue), DEB theory uses “Synthesizing Units” (SUs) (Kooijman, 1998). SUs are mathematically analogous to enzyme kinetics and scaled functional responses: substrates arrive at the SU (i.e., reserve-C and reserve-N), and are “bound” to the SU, which converts them to a new product at a certain rate. In the process of this product formation, the SU will also produce “waste”, which is critical in using DEB to consider consumer impacts on nutrient cycles, as this is what comprises nutrient release to the environment. Importantly, this is yet another point of connection between DEB models and EST. SUs in DEB models have been used to explore the effects of N-limitation on copepod reproduction (Kuijper, 2004), phosphorus limitation on Daphnia growth (Muller et al., 2001), C, N, and phosphorus limitation in marine cyanobacteria (Grossowicz et al., 2017), and multi-trophic nutrient limitation in chemostats (Kuijper et al., 2004).

Historically, both EST and DEB practitioners have created models that highlight how each framework individually can create stoichiometric-explicit models that predict life-history trajectories of consumers under various resource qualities (Sterner, 1997; Frost and Elser, 2002; Andersen et al., 2004; Kuijper et al., 2004; Logan et al., 2004; Anderson et al., 2005; Peace et al., 2013). We want to stress, however, that the use of both frameworks simultaneously can develop more realistic model representations of elemental flows, which are critical when models cross different mechanistic, temporal, and spatial scales, akin to real ecological systems. For example, the use of reserves in DEB models provides an internally consistent mathematical tracking of true nutrient limitation—a consumer has a buffer to withstand eating low-quality food for periods of time before it begins to experience fitness consequences (Kooijman, 2000). Thus, DEB models are ideal for understanding the effects of temporally variable resource quantity and quality. Additionally, DEB models more adequately model the individual-level trade-off between how assimilated nutrients eventually are utilized for growth, reproduction, or excreta, which can map onto different life history strategies (such as emergency reproduction, Kooijman, 2000). Alternatively, DEB models in and of themselves cannot adequately model or interpret the feedback of elements between consumers and their environment. EST is more well-suited to explain the ecological significance of consumer-driven nutrient cycling and can capitalize on the consistent currency of elemental ratios to scale from individuals to populations and biogeochemical cycles (Schindler and Eby, 1997; Hessen et al., 2004; Van De Waal et al., 2010; Moody et al., 2018; Elser et al., 2022; Rizzuto et al., 2024). Finally, DEB models have already been explicitly connected to NG (Kearney et al., 2010; Kearney, 2012), and many DEB models are able to predict movement and organismal performance outcomes that are relevant to many NG concepts (Arnall et al., 2019).

NG was developed as a framework for understanding how animals manage nutritional intake when faced with resources that are imbalanced in their macronutrients (e.g., carbohydrates, proteins) (Simpson et al., 1997). Additionally, NG can help understand driving forces of consumer-resource interactions, especially how these occur in spatially-explicit contexts. NG conceptualizes nutrient imbalances within a two-dimensional nutrient space. In this space, each animal has an intake target, reflecting its nutritional needs for processes such as maintenance and growth. As the animal consumes food, its position shifts within the nutrient space, depending on the nutrients it ingests (Anderson et al., 2020). In theory, an organism’s goal is to reach their intake target whilst limiting excess or deficient intake of any particular nutrient. Accordingly, NG can be used to make predictions regarding what behavioral strategies an organism may employ when faced with imbalanced diets. For example, populations of migratory Mormon crickets were found to cannibalize each other when they were deficient in protein and salt, but when these nutritional imbalances were relieved, cannibalism slowed (Simpson et al., 2006). Further, these specific imbalances were important drivers of the movement of these crickets through the landscape (Simpson et al., 2006), highlighting the important connections between individual nutrition and behavior.

While EST and NG are similar in that they both focus on nutrients, they use different measurements that have trade-offs in explaining ecological phenomena. NG directly examines macromolecules and connects biochemical nutrition to ecological dynamics (i.e., distinguishing between the carbon in lipids versus proteins), whereas EST uses elements as a proxy for nutrition and is thus able to directly compare data between different organisms and environments. The complementary concepts and distinctions between NG and EST have been further described elsewhere (Sperfeld et al., 2016, 2017; Anderson et al., 2020). Ultimately, NG’s focus on consumer feeding behaviors in response to resource variability is a critical aspect of modeling consumer-interactions in stoichiometrically imbalanced environments.

In our ABM, an individual consumer’s feeding decisions can be shaped by its elemental intake. Penalties for over- or under-consuming individual elements are implicitly prescribed by DEB and are dictated by an organism’s elemental reserves (i.e., C-rich and N-rich reserve). Thus, the ABM can be used to explore the fitness consequences of different feeding strategies—such as between diet-switching and compensatory feeding—as predicted by NG. For instance, under a selective feeding (diet-switching) strategy, a consumer may struggle to grow or reproduce if resource quality fluctuates among food choices (Anderson et al., 2020). Below, we outline the relevant concepts and important model decisions underlying the specification of the DEB model. All flux equations can be found in Table 1 , and a more complete specification of the model can be found in the Supporting Material (Appendices A and B), along with code for solving the equations in R (Appendix C). Moles of C, C-mols, are used as the common internal currency for the organism to streamline mass balances, as all internal state variables contain some amount of C. Moles of N, N-mols, can be easily calculated from the stoichiometric parameters.

Important DEB model decisions are listed briefly below.

The mass balance equations are represented as a system of ordinary differential equations (ODE) that model fluxes of nutrients ( J x y , mol C/day) (Equations 1–9). A diagram of these fluxes through the DEB state variables can be seen in Figure 2 . Users altering equations (or implementing other DEB models of interest) should take care that these fluxes are “absolute” fluxes as referred to in DEB theory, simply indicating that they are not scaled by consumer biomass or volume.

Somatic biomass increases due entirely to growth.

Reserve mass for both reserves (Equations 2, 3) increases through assimilation ( J A E i ) and recycling of reserves mobilized but not utilized in the growth and reproduction processes ( J G R E i ). Reserve mass is drained from the mobilization for growth and reproductive processes ( J E i C ).

The reproductive maturity level ( M H ) increases from the development/reproduction flux J R R until it reaches the threshold to adulthood ( M H P ). Once it has reached this threshold, J R R builds the reproductive buffer ( M R ) , which will produce offspring when it reaches the threshold for reproduction events ( M R 0 ).

Excretion of nutrients are released as a byproduct of the three metabolic processes that use multiple fluxes of different stoichiometries (somatic maintenance, growth, and development/reproduction). C-excretion fluxes are quantified in C-mols, and N-excretion fluxes are quantified in N-mols.

Egestion (“waste”) consists of ingested but unassimilated resources, and are also quantified in C-mols or N-mols, respectively.

All flux equations can be found in Table 1 , a guide for terms and units can be found in Table 2 , and all derivation steps/model choices can be found in detail within Appendix A in Supplementary Material.

NetLogo software (Wilensky, 1999) enables users to create spatially explicit ABMs and provides a user-friendly interactive coding environment to specify agent (i.e., consumer) behavior and metabolism, environmental conditions, with a well-established community of example models and code to users new to the program and ABMs as a whole. For longevity and new user ease, we have thus created our model in NetLogo, v6.4.

The ABM’s algorithms control the environment (patches) and consumers (agents). Algorithms are temporally divided into those that initialize the environment and consumers (“setup” procedures), and those that run each timestep (“tick”) of the model (“go” procedures), based on parameter specifications. For utility, many parameters that control (or impose) environmental and consumer mechanisms are displayed in the user interface. Below, we outline the procedures for the environment and consumers first in the initialization, then per timestep run in the model, and indicate how each of the three frameworks is either used to calculate values or contextualize assumptions for the procedures.

The NetLogo User Interface for this ABM contains three main sections: the parameter selection ( Figure 3A ), the environment (world) view ( Figure 3B ) and a set of plots and monitors that update each time step ( Figure 3C ).

The current model environment loads as a landscape with 25 × 25 patches, where the borders of the environment continuously wrap. The landscape is divided into patches that contain resources and those that do not. Resource patches can be further divided into a user-determined amount of high-quality clusters (e.g., heterogenous environment), or homogenous resource quality. The setup constrains the average C:N resource ratio over the whole world, within high-quality resource clusters, and lower-quality non-clustered resource patches, with user-determined target average, max, and min C:N ratios. We created spatial aggregation of resource patches and high-quality resource clusters by generating temporary agents at random patches, whose random walks determine which patches contain resources, and which resource patches are part of a high-quality cluster. For example, to assign high-quality resource clusters, the number of temporary agents generated is equal to the user-selected number of resource clusters. These agents take a total number of random walk steps that is equal to the user-selected area of each cluster. As these agents walk randomly, whichever patches they walk on get “converted” to the high-quality resource clusters. After performing these tasks, the agents are removed from the simulation. We also included a variable that incorporates random variation in all resource patches. Finally, we calculate the total C- and N-mol on each patch. For convenience, C-mol always is equivalent to the resource quantity parameter and N-mol is calculated from the C:N ratios of resources.

The total initial consumer population and proportion of adults vs juveniles are selected in the user interface. Once assigned a life stage, their initial DEB parameters are assigned, which the user specifies in the code. Then, initial state variable values are selected based on life stage, with user-specified variation or condition. Binary variables for consumers were also created to track important events, such as death (and cause of death), reproduction, and maturation.

A diagram summary of the ABM run procedures at each time step can be found in Figure 4 . Each timestep corresponds to a day, based on the parameter units of DEB terms. Thus, we have discretized the ODE within the ABM to calculate per day and implement these calculations at every timestep.

All consumer procedures are the first to run at each timestep. Within the consumer procedures, first, all timestep-specific parameters are reset. Next, consumers move, based on sensing how full their N-rich reserve is (scaled reserve density). If the reserve is below a user-specified threshold (“mismatch_threshold”), then the consumer seeks out the best quality patch within a user-specified radius (“sight”) of patches and moves towards this target patch at a user-specified walk distance (“distance-consumer-walks”). Otherwise, the consumers will walk either in a correlated random walk (Renshaw and Henderson, 1981) or a Brownian walk (Kac, 1947) motion, based on user selection. Once at a patch, consumers’ ingestion potential is calculated based on a Holling Type II Functional Response. To avoid a patch having a negative number of resources due to multiple consumers’ ingestion, we implemented a rule that when multiple consumers are on a patch and there are not enough resources for all of them, consumers eat in turn, from largest to smallest, until no resources remain. After moving, an individual consumer’s DEB sub-model will compute the change in their state variables at that time step. To minimize NetLogo code clutter and utilize more advanced numerical solving algorithms, we used the SimpleR Extension for NetLogo (Hovet et al., 2022) to send parameters and state variables to an R script which calculates the change in state variables using the DEB sub-model, then sends the results back to NetLogo. Within our R script, we used the rootSolve package (Soetaert and Herman, 2009) and the `multiroot()` function, which performs Newton-Raphson methods to numerically solve variables (see Appendix A in Supplementary Material for more details on which variables need numerical solving methods). The roots solved in these equations are stored for each individual consumer at each time point, so that they can be used as the initial guesses in the next timestep, improving stability and solving efficiency. Following these root-solving calculations, consumer variables and relevant parameters are updated in NetLogo (i.e., the ODE is implemented as discrete calculations each timestep, essentially the Euler solving method). Based on the updated consumer states, if an individual exceeded (triggered) their thresholds for reproduction or death, those events occurred. Consumers can die due to a background mortality rate (set in user interface), both reserves reaching less than 1% capacity, or shrinking their biomass too significantly (threshold set in user interface, but we used 40% shrinkage per other DEB ABMs (Martin et al., 2013).

Next, environmental patches run their procedures. First, patches update how many resources were consumed and how much N was added to the soil pool by consumers. Then, patches will grow resources at a set rate each timestep, so long as they have not reached carrying capacity (rate and carrying capacity set in user interface). Resources uptake N from the soil at a certain rate (set in the user interface). To ensure that the C:N ratio of resources is maintained within the user-specified minimum and maximum, resource growth is constrained so that the resource C:N ratio cannot rise above the maximum and soil-N uptake is constrained so that the resource C:N ratio cannot dip below the minimum ratio. Then, the total C- and N- mols for resources and soil are updated again. Importantly, all procedures are written so that there is a world conservation of N, to better model where a known set of N is distributed throughout the system. This also provides a method for checking simulation experiments for numerical instability (N is not conserved from beginning to end).

Finally, all miscellaneous procedures run, such as assigning patch color based on the resource quality and updating plots, before the next timestep begins.

We modeled our case study around snowshoe hares (L. americanus) in boreal forests of North America (Krebs, 2010). We chose boreal forests since they are known to be N-poor environments (Sponseller et al., 2016; Högberg et al., 2017). Resource N limitation is known to influence snowshoe hare ecology, behavior, and physiology, and snowshoe hares are capable of modifying resource availability through herbivory (Rizzuto et al., 2021). Our model was designed to examine how hare population dynamics and ecosystem N cycling are influenced by variability in resource quality and hare feeding behaviors. Accordingly, we ran four distinct model scenarios, varying both environmental and consumer dynamics: (1) homogeneous vs. heterogeneous C:N molar ratios across the environment (see Sections 4.2 and 5.2) and (2) hare populations with either random feeding or N-selective feeding under N limitation (see Section 5.3).

In all scenarios, the average resource C:N ratio initialized at 32, representing higher resource quality for relevant hare forage species (Rizzuto et al., 2019). Additionally, the maximum and minimum C:N resource ratios were arbitrarily selected to be 60 and 30, and the percent of world without any resources set to 30%. Under a heterogeneous resource conditions, the number of clusters was set to 20. Consumer C:N somatic biomass was set at 5, selected from a range of snowshoe hares at different ages (Rizzuto et al., 2019). For calculation convenience, the reproductive biomass C:N was set to be the same as the somatic biomass.

We ran simulations across a factorial cross of two parameters: sigma (σ), which controls the multiplicative additional somatic maintenance costs of movement and also the N-rich reserve density threshold level at which consumers switch from a random correlated walk movement (random feeding) to seeking high quality food under N limitation (selective feeding, or, diet switching). For simplicity, we selected high values for both parameters: maintenance penalties for movement of 60% (sigma = 0.6) and thresholds for selective movement of 70% full reserves (scaled reserve density = 0.7). Other ABM parameters that were held constant are found in Appendix B in Supplementary Material.

Primary parameters for the DEB model were taken from the DEB Add-My-Pet Portal and used in a mass context with C-mols as common currency (Appendix B in Supplementary Material). While parameters for specifically L. americanus were not available in the database, we selected parameters from another cold-adapted species within the genus, the mountain hare (L. timidus). From these primary parameters, we calculated any other secondary DEB parameters using the energy-mass relationships described in Table 2 and “Notation and Symbols” section in Kooijman (2000).

To account for the fact that the primary parameters were derived for a single-reserve DEB model, we split the total maximum reserve density between the maximum reserve density of both our reserves: 40% of the total maximum reserve density was arbitrarily chosen for the maximum C-rich reserve, while the rest was used for the N-rich reserve. Additionally, we lowered certain parameters to generate qualitatively realistic dynamics around maturation and reproduction: maturation threshold for puberty, the size of an individual at birth, and thus needed to lower the cost of reproduction. We did this in part to ensure that we were more closely matching reproduction yearly dynamics with field observations of snowshoe hares (Snowshoe Hare Species Profile, Alaska Department of Fish and Game, n.d.). We calculated the cost of reproduction as the needed amount of C and N exactly to the size of a newborn plus its precocial maturity level (which we kept the same as the AMP Portal parameter) and reserves. For simplicity, consumers were all born with these same traits and scaled reserve density (0.9 for C-reserve and calculated ~0.86 for N-reserve). Births occurred when an adult consumer’s reproduction buffer reached the level of having a litter of 5. A full table of parameters used is in Appendix B in Supplementary Material, and the code to run the ABM is available in Appendix C in Supplementary Material.

Within the model, N is cycled from resources to consumers through ingestion, then to a soil pool through consumer release (egestion and excretion), then back to the resources through uptake ( Figure 1 ). We aimed to have a closed N system, though it is not necessarily realistic, so that we could ensure the model was numerically stable and track how N is shifted in different biomass throughout the simulation.

We used BehaviorSpace in NetLogo to perform these simulation experiments. Each simulation ran for 1000 time steps (i.e., 1000 days), with 10 replicates per parameter set. Replicates were needed to identify trends even among identical parameter sets, as consumers’ movement and the initialization of resource aggregation introduces stochasticity. Simulations ended early if no consumers were alive. The selected outputs were collected at each timestep from the simulations: total individuals as adults and juveniles, total resource quantity, total soil N, total resource N, and total consumer N, as well as the mean and variance of resource C:N, resource quantity, consumer size, consumer scaled reserve densities, consumer reproduction and maturity. We removed 3 replicates from the full set of simulations, due to the simulation becoming numerically unstable mid-run, potentially due to hardware constraints.

Across all model simulations, hare populations exhibited rapid population growth following model initialization, presumably driven by initially high resource availability ( Figure 5 ). However, reproduction was slightly lower for hare populations that fed randomly in the environment with homogeneous resource quality ( Figures 5A, B ). Although adult hare numbers temporarily declined after the initial reproductive surge, populations rebounded as juvenile growth increased in all simulations ( Figure 5C, D ). While population growth (adult count) was similar between scenarios with homogenous and heterogenous resource quality, hare populations that fed selectively on N-rich resource patches during periods of N-limitation exhibited abundances ~2× higher than their randomly feeding counterparts ( Figures 5C, D ).

While the heterogeneity of resource quality had only minor effects on hare population dynamics, it had important influence on N cycling and distribution by hare populations. Total resource consumption was higher for selectively feeding populations which exhibited greater abundances ( Figures 5E, F ) but was specifically highest for selectively feeding populations in the homogenous environment. Resource C:N ratios—reflecting resource quality—rapidly increased following initial plant growth but returned to baseline values following plant consumption during reproductive surges. However, resource C:N ratios eventually increased from starting values across simulations, reflecting an overall decline in resource quality ( Figures 5G, H ). Interestingly, this was most pronounced for selectively feeding populations, presumably due to the fact that these hares were able to seek out and consume N-rich resource patches.

Hare populations played an important role in generating and maintaining heterogeneity in resource quality—represented by the variance of N across the landscape ( Figures 5I, J ). Notably, hares dramatically amplified N heterogeneity in the initially homogeneous environment, and only marginally decreased N variance in the starting heterogenous environment. These effects occurred irrespective of feeding behaviors, suggesting a general importance of N transport by hare populations. As hares traversed the homogeneous environment, their waste and carcass deposition generated new high-N patches, amplifying the contrast between N-rich and N-poor areas. In contrast, within the heterogeneous environment, hares redistributed N from initial hotspots—first increasing heterogeneity by creating additional high-N patches but ultimately reducing heterogeneity below initial levels as their movement led to N deposition across the landscape ( Figure 5I ).

Ultimately, our model results suggest that the ability for hares to seek out N-rich resources when faced with N limitation significantly enhanced population fitness. While the initial heterogeneity of resource quality had only modest effects on hare population dynamics, its interaction with hare behavior strongly influenced patterns of N cycling. In general, hares reduced average resource quality through N consumption and increased spatial heterogeneity of resource quality. This case study thus also highlights the ability for researchers to study importance of consumers in distributing nutrients throughout the landscape.

We demonstrate how a multi-reserve DEB model within an ABM framework can be used to simultaneously explore the dynamic feedbacks among variable nutrient conditions, consumer traits, and population structure. These processes occur while maintaining physiological, spatial, and temporal specificity, akin to the complexity of real ecological systems. By parameterizing the DEB model using stoichiometric data and allowing for simulated consumers to behaviorally respond to their internal state, this approach provides a new integrative exploration of EST, NG, and DEB concepts. While we recognize that the general model outputs are highly dependent on the selection of parameter values (e.g., hare structural shrinkage allowed, environment spatial structure) and simplified ABM procedures (e.g., not allowing hares to remain in one patch over multiple timesteps), our primary aim here is to demonstrate how the model can help develop hypotheses about how stoichiometric imbalance influences populations and ecosystems. Furthermore, we believe it is most helpful to provide a basic set of procedures and broad parameter manipulation abilities so that other researchers can tune the model for different systems.

The results from our case study demonstrate the functionality of the model yet only touch the surface of the model’s utility to examine a broad range of questions across DEB, EST, and NG based on real world scenarios. Many novel anthropogenic changes to ecosystems, such as climate change, may affect consumer nutrition, energetics, and behaviors simultaneously (Van De Waal et al., 2010; Toseland et al., 2013; Diehl et al., 2022). Additionally, consumers have concomitant impacts on ecosystem functioning, indicating the importance of studying the impacts of anthropogenic change on zoogeochemistry (Abraham et al., 2023). Theoretical models offer valuable opportunities for understanding these challenges by generating testable hypotheses and making novel predictions. Such models are particularly useful when empirical data is difficult to obtain due to time, spatial scale, or logistical constraints. In cases where resources are limited, models can also help identify which variables are most critical to measure. Additionally, theoretical models allow for the exploration of complex feedback mechanisms that may be challenging to detect through experimental or observational approaches.

Despite the intricacy and utility of our model, there remain many limitations that also highlight key areas of future integration. First, our model currently only reflects a two-trophic system, with one consumer population and a generalized primary producer community. Expanding the model to include predators and a more diverse consumer community would allow for a deeper understanding of how biodiversity and trophic complexity alter and are influenced by stoichiometric imbalances. For example, predator-prey interactions have been shown to substantially affect consumer population dynamics, movement and behavior, energetics, and nutrient cycling within ecosystems (Denno and Fagan, 2003; Schmidt et al., 2012; Kohl et al., 2015; Herzog et al., 2023). On the other hand, predator risk affects the physiology of prey which may result in higher metabolic rates as well as habitat shifts to fulfill their nutrient demands (Guariento et al., 2018). For example, grasshoppers under high predation risk increase requirements for energy from carbohydrates affecting the C of their food source as well as their C:N release by excretion (Schmitz et al., 2010). Additionally, while the model is already capable of measuring how consumers influence spatial availability of nutrients, it could be further developed to include specific spatial statistics (e.g., spatial autocorrelation, such as Morran’s I), an open-nutrient exchange landscape (e.g., nitrogen fixation, as opposed to the closed-nitrogen environment currently in the model), or environmental responses and feedbacks (e.g., primary productivity).

Second, our model currently only incorporates C:N ratios, reflecting prominent trade-offs between C-rich carbohydrates (energy) and N-rich proteins (nutrition). However, both the ABM and DEB sub-model are capable of modelling other elemental ratios, such as C:P, N:P, or micronutrients (i.e., Na, vitamins). For example, the reciprocal relationship between ribosome allocation and protein synthesis concurrently increases demands for P to maintain demands for growth (Hessen, 2008). Additionally, other micronutrients, such as Na, Ca, K, and Mg are gaining increasing attention for their contribution to ecophysiology and ecosystem functioning (Kaspari and Powers, 2016; Duvall et al., 2023). The immense amounts of elemental stoichiometry data including emerging data repositories, such as the STOICH database (STOICH Database Home, n.d), may be used to examine how a plethora of ecological outcomes differ under various simulated stoichiometric imbalance scenarios.

Third, our model can be used to examine theoretical outcomes under different behavioral scenarios and nutritional trade-offs, informed by NG approaches. Although our case study focused on consumer diet-switching responses, we developed the model so that a user can select between diet-switching and compensatory ingestion strategies. For example, shrimps can change their diet in low-P streams from C-rich detritus to P-rich insects and periphyton (Snyder et al., 2015), and many large vertebrates will constantly select for N, P, and Na-rich plants and other materials under nutrient stress (Duvall et al., 2023; Monk et al., 2024),. However, a consumer may also elevate their feeding rates and compensate for the insufficient food quality by higher food quantity (Jensen et al., 2012; Suzuki-Ohno et al., 2012; Flores et al., 2014; Anderson et al., 2020). In our ABM, users could compare the outcomes of compensatory feeding versus diet-switching behavior. For example, while diet switching may induce increased metabolic demands associated with food processing or increased movement, compensatory feeding may induce metabolic consequences associated with excess nutrient intake and processing. Adding these additional behavioral strategies can be used to examine which strategies are, in theory, more evolutionarily advantageous, and how these advantages may differ under various stoichiometric imbalance scenarios, especially those that play out in spatially heterogeneous landscapes (e.g., high vs low spatial heterogeneity of nutrients (Smith and Moore, 2003), high vs low N digestibility (Weiser et al., 1997).

Ecological systems are highly complex. Efforts to understand the intricacy of ecological interactions have ranged from generalized frameworks to complex ecological models. There are inherent trade-offs in each approach: general frameworks provide an intuitive way of describing nature, whereas complex models can attempt to provide realistic predictions that elucidate the importance of specific variables and mechanisms in ecological interactions (Sun et al., 2016). We recognize that our model lies on the latter end of this spectrum. Indeed, a general frustration with both ABMs and DEB models is their complexity and the number of parameters and individual-level processes included, creating compounding uncertainty in model outcomes (An et al., 2021). Regardless, the aim of this integrative model is to more accurately account for multiple processes when assessing the scaling effects of elemental imbalance—an outcome achievable only through increased model complexity. Since we aimed to make our model as user-friendly as possible to pull together researchers from different disciplines, we have tried to provide detailed descriptions of how we arrived at theoretical decisions, placed many important ABM parameters in the NetLogo User Interface, and designed the ABM and DEB sub-model to be mostly modular, such that researchers would be able to use the DEB model without needing to also use the ABM. Future research should compare our model outcomes with simpler models of DEB, EST, and NG to ascertain whether individual frameworks alone fail to produce the outcomes of the ABM as hypothesized.

The quality of model outcomes is also determined by the respective quality of empirical data needed for parameterization, particularly when trying to replicate real ecological systems. It is important that one must know enough about the biology of the organism to generate reasonable relationships between the parameters of two focal nutrients, or instead, have a tractable model system to collect enough of these parameters experimentally (e.g., Daphnia). DEB parameters can be estimated from empirical data (van der Meer, 2006), but are also available in a large database, Add-My-Pet portal (Marques et al., 2018). Users can install packages for downloading entries in the portal and work with their data through Matlab, R, and other programs. While many basic parameters can be extracted from the DEB Add-My-Pet portal for DEB models, it is still likely that one of the above strategies for parameter selection will need to be employed. Additionally, as more reserves, structures, and/or ontogenetic stages are added to the DEB model, the parameters increase, so parsimony in model design is needed for tractable calculations. Beyond consumer parameters, it can also be difficult to acquire sufficient environmental data for modeling abiotic and biotic nutrient cycling processes. For example, a challenge we faced in setting up our case study was determining the rates of nitrogen incorporation into soils and uptake by plants, for which data is lacking and can be highly context-dependent. Furthermore, users will need to appropriately tune the spatial and temporal scale of the model to match their respective system, as this will have important impacts on individual-level processes and their emergent effects as populations and communities. For example, in our case study, we use ‘days’ as a time unit based on the specific DEB parameters we chose, yet the ideal timescale will depend on focal systems, organisms, and questions.

In our model, consumers actively shape their environment by altering the distribution of elements and resources through space and time. A key advantage of the agent-based modeling (ABM) approach is its ability to simultaneously quantify the interactions between consumers and their environment. Variations in resource quantity and quality driven by organisms can significantly influence broader ecosystem functions, such as primary productivity, plant community composition, and biodiversity. These changes can create greater variability in the resources encountered at the individual level, generating long-term feedbacks that affect consumer populations.

Measuring the impacts of consumers on environmental heterogeneity is often challenging or even impossible in the real world at proper scale. For example, predicting consumer impacts on nutrient cycling must account for discernment between different nutrient release pathways (e.g., excrement, egestion, carcass deposition) including their variable rates and quantities. These measurements are largely unobtainable using field or observation methods, and thus, the ABM approach can provide a step towards a more realistic and holistic perspective on how consumers contribute to the heterogeneity of their environments, with bottom-up effects on ecosystem structure and function (e.g., primary productivity, habitat, biodiversity).

Ultimately, the integrative approach of our model presents new opportunities to explore concepts from EST, NG, and DEB, as well as other frameworks. The future directions we present here are only a selection of possibilities and opportunities. Continuation of model development by other users will ensure that modules are refined for specific landscapes (e.g., aquatic habitats), and how the wealth of stoichiometric, nutritional, and ecological data collected by our community can be combined to inform understanding and future study designs.