Nutritional state-structured model

We begin by defining the nutritional state-structured population model, where the consumer population is partitioned into two states: (a) an energetically replete (full) state F, where the consumer reproduces at a constant rate λ and does not die from starvation, and (b) an energetically deficient (hungry) state H, where the consumer does not reproduce but dies by starvation at rate μ. The dynamics of the underlying resource R are governed by logistic growth with an intrinsic growth rate α and a carrying capacity C. The rate at which consumers transition between states and consume resources is dependent on their number, the abundance of resources, the efficiency of converting resources into metabolism, and how that metabolism is partitioned between maintenance and growth purposes. We provide a physiologically and energetically mechanistic model for each of these dynamics and constants (Methods), and show that the system produces a simple nondimensional form that we describe below.

Consumers transition from the full state F to the hungry state H at a rate σ—the starvation rate—and also in proportion to the absence of resources (1−R) (the maximum resource density has been nondimensionalized to 1; see Methods). Conversely, consumers recover from state H to state F at rate ξρ and in proportion to R, where ξ represents a ratio between maximal resource consumption and the carrying capacity of the resource. The resources that are eaten by hungry consumers (at rate ρR+δ) account for their somatic growth (ρR) and maintenance (δ). Full consumers eat resources at a constant rate β that accounts for maximal maintenance and somatic growth (see Methods for mechanistic derivations of these rates from resource energetics). The NSM represents an ecologically motivated fundamental extension of the idealized starving random walk model of foraging, which focuses on resource depletion, to include reproduction and resource replenishment33,34,35, and is a more general formulation than previous models that incorporate starvation36,37.

In the mean-field approximation, in which the consumers and resources are perfectly mixed, their densities are governed by the rate equations

$$\begin{array}{*{20}{l}} {\dot F} \hfill & { = \lambda F + \xi \rho RH - \sigma \left( {1 - R} \right)F,} \hfill \\ {\dot H} \hfill & { = \sigma \left( {1 - R} \right)F - \xi \rho RH - \mu H,} \hfill \\ {\dot R} \hfill & { = \alpha \left( {1 - R} \right)R - \left( {\rho R + \delta } \right)H - \beta F.} \hfill \end{array}$$ (1)

This system of nondimensional equations follows from a set of first-principle relationships for resource consumption and growth (see Methods for a full derivation and the dimensional form). Notice that the total consumer density F+H follows \(\dot F + \dot H = \lambda F - \mu H\). This resembles the equation of motion for the predator density in the LV model38, except that the resource density does not appear in the growth term. The rate of reproduction is independent of resource density because the full consumer partitions a constant amount of energy toward reproduction, whereas a hungry consumer partitions no energy toward reproduction. Similarly, the consumer maintenance terms (δH and βF) are also independent of resource density because they represent a minimal energetic requirement for consumers in the H and F state, respectively.

Steady states of the NSM

From the single internal fixed point (Eq. 5, see Methods), an obvious constraint on the NSM is that the reproduction rate λ must be less than the starvation rate σ, so that the consumer and resource densities are positive. The condition σ = λ represents a transcritical (TC) bifurcation39 that demarcates a physical from an unphysical (negative steady-state densities) regime. The biological implication of the constraint λ < σ has a simple interpretation—the rate at which a macroscopic organism loses mass due to lack of resources is generally much faster than the rate of reproduction. As we will discuss below, this inequality is also a natural consequence of allometric constraints for organisms within empirically observed body size ranges.

In the physical regime of λ < σ, the fixed point (Eq. 5) may either be a stable node or a limit cycle (Fig. 1). In continuous-time systems, a limit cycle arises when a pair of complex conjugate eigenvalues crosses the imaginary axis to attain positive real parts40. This Hopf bifurcation is defined by Det(S) = 0, where S is the Sylvester matrix, which is composed of the coefficients of the characteristic polynomial of the Jacobian matrix41. As the system parameters are tuned to be within the stable regime, but close to the Hopf bifurcation, the amplitude of the transient cycles becomes large. Given that ecological systems are constantly being perturbed42, the onset of transient cycles, even though they decay with time in the mean-field description, can increase extinction risk43,44,45.

Fig. 1 The transcritical (TC; dashed line) and Hopf bifurcation (solid line) as a function of the starvation rate σ and recovery rate ρ for a 100-g consumer. These bifurcation conditions separate parameter space into unphysical (left of the TC), cyclic, and steady-state dynamic regimes. The colors show the steady-state densities for the energetically replete consumer F* Full size image

When the starvation rate σ >> λ, a substantial fraction of the consumers is driven to the hungry nonreproducing state. Because reproduction is inhibited, there is a low steady-state consumer density and a high steady-state resource density. However, if σ/λ→1 from above, the population is overloaded with energetically replete (reproducing) individuals, thereby promoting transient oscillations between the consumer and resource densities (Fig. 1). If the starvation rate is low enough that the Hopf bifurcation is crossed, these oscillations become stable. This threshold occurs at higher values of the starvation rate as the recovery rate ρ increases, such that the range of parameter space giving rise to cyclic dynamics also increases with higher recovery rates.

The allometry of extinction risk

While there are no a priori constraints on the parameters in the NSM, we expect that each species should be restricted to a distinct portion of the parameter space. We use allometric scaling relations to constrain the covariation of rates in a principled and biologically meaningful manner (Methods). Allometric scaling relations highlight common constraints and average trends across large ranges in body size and species diversity. Many of these relations can be derived from a small set of assumptions. In Methods, we describe our framework to determine the covariation of timescales and rates across a range of body sizes for each of the key parameters of our model (cf. ref. 46).

Nearly all of the rates described in the NSM are determined by consumer metabolism, which can be used to describe a variety of organismal features47. We derive, from first principles, the relationships for the rates of reproduction, starvation, recovery, and mortality as a function of an organism’s body size and metabolic rate (Methods). Because we aim to explore the starvation-recovery dynamics as a function of an organism’s body mass M, we parameterize these rates in terms of the percent gain and loss of the asymptotic (maximum) body mass, εM, where different values of ε define different states of the consumer (Fig. 2; see Methods for derivations of allometrically constrained rate equations). Although the rate equations (1) are general and can in principle be used to explore the starvation-recovery dynamics for most organisms, here, we focus on allometric relationships for terrestrial-bound lower-trophic-level endotherms (Table 1), specifically herbivorous mammals, which range from a minimum of M ≈ 1 g (the Etruscan shrew Suncus etruscus) to a maximum of M ≈ 107 g (the early Oligocene Indricotheriinae and the Miocene Deinotheriinae). Investigating other classes of organisms would simply involve altering the metabolic exponents and scalings associated with ε. Moreover, we emphasize that our allometric equations describe mean relationships and do not account for the (sometimes considerable) variance associated with individual species. We note that including additional allometrically scaled mortality terms to both F and H does not change the form of our model nor impact our quantitative findings (Supplementary Note 1).

Fig. 2 The growth trajectory over absolute time of an individual organism as a function of body mass. Initial growth follows the black trajectory to an energetically replete reproductive adult mass of m = ε λ M (Methods). Starvation follows the red trajectory to m = ε σ ε λ M. Recovery follows the green curve to the replete adult mass, where this trajectory differs from the original growth because only fat is being regrown that requires a longer time to reach ε λ M. Alternatively, death from starvation follows the blue trajectory to m = ε μ ε λ M Full size image

Table 1 Parameter values for mammals Full size table

As the allometric derivations of the NSM rate laws reveal (Methods), starvation and recovery rates are not independent parameters, and the biologically relevant portion of the phase space shown in Fig. 1 is constrained via covarying parameters. Given the parameters of terrestrial endotherms, we find that the starvation rate σ and the recovery rate ρ are constrained to lie within a small region of potential values for the known range of body size M. Indeed, starvation and recovery rates across all values of M fall squarely in the steady-state region at some distance from the Hopf bifurcation. This suggests that cyclic population dynamics should be rare, particularly in resource-limited environments.

Higher rates of starvation result in a larger flux of the population to the hungry state. In this state, reproduction is absent, thus increasing the likelihood of extinction. From the perspective of population survival, it is the rate of starvation relative to the rate of recovery that determines the long-term dynamics of the various species (Fig. 1). We therefore examine the competing effects of cyclic dynamics vs. changes in steady-state density on extinction risk, both as functions of σ and ρ. To this end, we computed the probability of extinction, where we define extinction as a population trajectory falling below one-fifth of the allometrically constrained steady state at any time between t = 108 and t = 1010. This procedure was repeated for 50 replicates of the continuous-time system shown in Eq. 1 for organisms with mass ranging from 102 to 106 g. In each replicate, the initial densities were chosen to be (XF*, XH*, R*), with X a random variable uniformly distributed in [0, 2]. By allowing the rate of starvation to vary, we assessed extinction risk across a range of values for σ and ρ between ca. 10−8 and 10−3. Higher rates of extinction correspond to both large σ if ρ is small, and large ρ if σ is small. In the former case, increased extinction risk arises because of the decrease in the steady-state consumer population density (Figs. 1, 3). In the latter case, the increased extinction risk results from higher-amplitude transient cycles as the system nears the Hopf bifurcation (Fig. 3). This interplay creates an “extinction refuge”, such that for a constrained range of σ and ρ, extinction probabilities are minimized.

Fig. 3 Probability of extinction for a consumer with a M = 102 g and b M = 106 g as a function of the starvation rate σ and recovery rate ρ, where the initial density is given as (XF*, XH*, R*), where X is a random uniform variable in [0, 2]. Note the change in scale in b. Extinction is defined as the population trajectory falling below 0.2× the allometrically constrained steady state. The white points denote the allometrically constrained starvation and recovery rates for consumers of each body size Full size image

We find that the allometrically constrained values of σ and ρ, each representing different trajectories along the ontogenetic curve (Fig. 2), fall squarely within the extinction refuge across a range of Μ (Fig. 3a, b, white points). These values are close enough to the Hopf bifurcation to avoid low steady-state densities, yet distant enough to avoid large-amplitude transient cycles. Allometric values of σ and ρ fall within this relatively small window, which supports the possibility that a selective mechanism has constrained the physiological conditions driving starvation and recovery rates within populations. Such a mechanism would select for organism physiology that generates appropriate σ and ρ values that minimize extinction risk. This selection could occur via the tuning of body fat percentages, metabolic rates, and/or biomass maintenance efficiencies. We also find that as body size increases, the size of the low extinction-risk parameter space shrinks (Fig. 3b), suggesting that the population dynamics for larger organisms are more sensitive to variability in physiological rates. This finding is in accordance with, and may serve as contributing support for, observations of increased extinction risk among larger mammals48.

Damuth’s law and body size limits

The NSM correctly predicts that smaller species have larger steady-state population densities (Fig. 4). Similar predictions have been made for carnivore populations using alternative consumer-resource models49. Moreover, we show that the NSM provides independent theoretical support for Damuth’s law25,26,27,28. Damuth’s law shows that species abundances, Ν*, follow Ν* = 0.01Μ−0.78 (g m−2). Figure 4 shows that both F* and Η* scale as M−η, with η ≈ 3/4, over a wide range of organismal sizes and that F*+ Η* closely matches the best fit to Damuth’s data. Remarkably, this result illustrates that the steady-state values of the NSM combined with the derived timescales naturally give rise to Damuth’s law. While the initial metabolic studies supporting Damuth’s law provide arguments for the value of the exponent25,26 (Supplementary Note 2), our model predicts not only the exponent but also the normalization constant dependencies by explicitly including the resource dynamics and the parameters that determine growth and consumption. These predictions are complementary to recent work that also predicts the exponent and normalization constant of density relationships from the detailed allometries of reproduction, capture area, conversion efficiency, and mortality within predator–prey dynamic models49,50. It should be noted that density relationships of individual clades follow a more shallow scaling relationship than that predicted by Damuth’s law28. In the context of our model, this finding suggests that future work may be able to anticipate these shifts by accounting for differences in the physiological parameters associated with each clade.

Fig. 4 Consumer steady states F* (green) and H* (orange) as a function of body mass along with the data from Damuth25. Inset: resource steady state R* as a function of consumer body mass Full size image

With respect to predicted steady-state densities, the population sizes of both F and H go to zero at a finite body size, where the steady-state resources also vanish (Fig. 4). This behavior is governed by the body size at which f 0 Mγ−1+u 0 Mζ−1 = 1 (Supplementary Note 3) that causes the death rate to vanish, μ = 0, and corresponds to (F*, H*, R*) = (0, 0, 0). Just before this point H* becomes large, in accordance with the asymptotic behavior in Fig. 4. This point predicts an upper bound on mammalian body size at M max = 6.54 × 107 g. Moreover, M max , which is entirely determined by the population-level consequences of energetic constraints, is within an order of magnitude of the maximum body size observed in the North American mammalian fossil record29, as well as the mass predicted from an evolutionary model of body size evolution30. We emphasize that the asymptotic behavior and predicted upper bound depend only on the scaling of body composition and are independent of the resource parameters. The prediction of an asymptotic limit on mammalian size parallels work on microbial life where an upper and lower bound on bacterial size, and an upper bound on single-cell eukaryotic size, is predicted from similar growth and energetic scaling relationships3,51. It has also been shown that models that incorporate the allometry of hunting and resting combined with foraging time predicts a maximum carnivore size between 7 × 105 and 1.1 × 106 g52,53. Similarly, the maximum body size within a particular lineage has been shown to scale with the metabolic normalization constant54. This complementary approach is based on the balance between growth and mortality, and suggests that future connections between the scaling of fat and muscle mass should systematically be connected with B 0 when comparing lineages.

A mechanism for Cope’s rule

Metabolite transport constraints are widely thought to place limits on biological scaling47,55 and thereby lead to specific predictions on the minimum possible body size for organisms56. Above this bound, a number of energetic and evolutionary mechanisms have been explored to assess the costs and benefits associated with larger body masses, particularly for mammals. One important such example is the “fasting endurance hypothesis”, which contends that larger body size, with consequent lower metabolic rates and increased ability to maintain more endogenous energetic reserves, may buffer organisms against environmental fluctuations in resource availability57. Over evolutionary time, terrestrial mammalian lineages show a significant trend toward larger body size—Cope’s rule29,30,31,32. It is thought that within-lineage drivers generate selection toward an optimal upper bound of roughly 107 g29, a value that is likely limited by higher extinction risk for large taxa over longer timescales30. These trends are thought to be driven by a combination of climate change and niche availability32; however, the underpinning energetic costs and benefits of larger body sizes, and how they influence dynamics over ecological timescales, have not been explored.

The NSM predicts that the steady-state resource density R* decreases with increasing body size of the consumer population (Fig. 4, inset), and classic resource competition theory predicts that the species surviving on the lowest resource abundance will outcompete others58,59,60. Thus, the combined NSM steady-state dynamics and allometric timescales (Eq. 7) predict that larger mammals have an intrinsic competitive advantage given a common resource.

However, the above resource relationships do not offer a mechanism for how body size is selected. We directly assess competitive outcome between two closely related species: a resident species of mass M, and a competing species (denoted by ′) where individuals have a different proportion of body fat such that M′ = M(1+χ). For χ < 0, the competing individuals have fewer metabolic reserves than the resident species and vice versa for χ > 0. For the allowable values of χ (Methods), the mass of the competitor M′ should exceed the minimal amount of body fat, 1 + χ > ε σ , and the adjusted time to reproduce must be positive, which, given Eq. 7, implies that \(1 - \varepsilon _{\mathrm{\lambda }}^{1 - \eta }\left( {1 + \chi } \right)^{1 - \eta } \, > \, 0\). These conditions imply that \(\chi \in ( - f_0M^{\gamma - 1},1/\varepsilon _{\mathrm{\lambda }} - 1)\) where the upper bound approximately equals 0.05 and the lower bound is mass-dependent. The modified mass of the competitor leads to altered rates of starvation σ(M′), recovery ρ(M′), and the maintenance of both starving δ(M′) and full consumers β(M′) (see Methods for derivations of competitor rates). Importantly, ε σ , which determines the point along the growth curve that defines the body composition of starved foragers, is assumed to remain unchanged for the competing population.

To assess the susceptibility of the resident species to competitive exclusion, we determine which consumer pushes the steady-state resource density R* to lower values for a given value of χ, with the expectation that a population capable of surviving on lower-resource densities has a competitive advantage58. We find that for M ≤ 1.748 × 107 g, having additional body fat (χ > 0) results in a lower steady-state resource density (R′* < R*), such that the competitor has an intrinsic advantage over the resident species (Fig. 5). However, for M > 1.748 × 107 g, leaner individuals (χ < 0) have lower-resource steady-state densities.

Fig. 5 Competitive outcomes for a resident species with body mass M vs. a closely related competing species with modified body mass M′ = M(1 + χ). The blue region denotes proportions of modified mass χ resulting in exclusion of the resident species. The red region denotes values of χ that result in a mass that is below the starvation threshold and are thus infeasible. Arrows point to the predicted optimal mass from our model M opt = 1.748 × 107, which may serve as an evolutionary attractor for body mass. The black wedge points to the largest body mass known for terrestrial mammals (Deinotherium spp.) at 1.74 × 107 g31 Full size image

The observed switch in susceptibility as a function of χ at M opt = 1.748 × 107 g thus serves as an attractor, such that the NSM predicts organismal mass to increase if M < M opt and decrease if M > M opt . This value is close to but smaller than the asymptotic upper bound for terrestrial mammal body size predicted by the NSM, and is remarkably close to independent estimates of the largest land mammals, the early Oligocene Indricotherium at ≈1.5 × 107 g and the late Miocene Deinotherium at ≈1.74 × 107 g31. Additionally, our calculation of M opt as a function of mass-dependent physiological rates is similar to theoretical estimates of maximum body size30, and provides independent theoretical support for the observation of a “maximum body size attractor” explored by Alroy29.