Predation is a biological interaction that impacts community structure (Marquet, 2002; Hatton et al., 2015; Brose et al., 2019) and due to the extinction of large carnivores (Ripple and Beschta, 2007; Atkins et al., 2019) it is reasonable to expect profound modifications to ecosystems (Ripple et al., 2014). Predator density is usually limited by metabolic rate and prey abundance (Carbone and Gittleman, 2002), and the coexistence between predators and prey requires mechanisms that allow for the maintenance of an adequate level of resources for each species to persist (Chesson, 2000). Thus, characterizing population fluctuations through demographic patterns (Barbier and Loreau, 2019) allows assessing the energy requirements used by a community (Damuth, 1987, 1993; Brown, 1995), and the indirect effects of top predators on the basal resource through interaction with some intermediate species, i.e., trophic cascade (Terborgh et al., 2001, 2010). Determining how these phenomena arise from energetic, evolutionary and physiological constraints and inter-species interactions is necessary to understand community structuring processes (Pawar et al., 2015; Momo, 2017; Piovia-Scott et al., 2017).

It has been recorded that body size is allometrically related with basal metabolism, reproduction rate, mortality rate and population abundance. Mathematically it can be expressed as: Y = b·m^α, where Y is one these biological attributes of interest, b is a species-dependent factor, m is the body size and α is the allometric exponent (Kleiber, 1932; Damuth, 1981; Hatton et al., 2019). For instance, a demographic pattern recorded in empirical research is the “size-density relationship,” which indicates a negative relationship between density and average body size (Damuth, 1987, 1993, 2007; Marquet et al., 1990). Interestingly, a recent study suggests that the relationship between body size and population density has changed over a relatively short period of time, a change that the allometric exponent can capture (Santini and Isaac, 2021). This could eventually change body size-dependent coexistence conditions, even impacting the energy requirements used by the species. Furthermore, the indirect phenomena produced in food webs could intensify (Delong et al., 2015). Barbier et al. (2021), using a macroecological approach performed a meta-analysis of predator-prey pairs (mammals, birds and reptiles), and showed that predation rates at the macro level differ from that presented in Lotka-Volterra-type models (the probability of encounters grows as the product of the density of the predator and its prey).

On the other hand, predator-prey interactions have been intensively studied through mathematical modeling (May, 1973; Yodzis and Innes, 1992). Weitz and Levin (2006) using a Rosenzweig-MacArthur type model, incorporating body masses in the parameters obtained the size-density relationship with allometric exponent 3/4. There are strong similarities between empirical patterns and those predicted by mathematical models for carnivorous mammals (Carbone and Gittleman, 2002; Weitz and Levin, 2006; DeLong and Vasseur, 2012). Extending the modeling approach to a three-level food chain (apex predator-consumer-prey) may allow to improve the understanding of fundamental aspects of community structuring, such as species coexistence and the effect of indirect interactions (Rinaldi and De Feo, 1999; Terborgh et al., 2010; Abdala-Roberts et al., 2019).

Size-density relationships are essential for understanding ecosystem functioning, and dynamic models provide us with an adequate description of natural systems. This study presents a theoretical exploration using mathematical models, considering the law of mass action, that integrates species body sizes into a dynamic representing a three-species food chain. The relationships between species' body sizes eventually determine species' coexistence and energetic requirements. In addition, this trait allows identifying the influence and control that species from higher trophic levels exert on species at lower trophic levels. Through mathematical modeling, it is possible to show relationships in simple terms that allow the notions of coexistence, energy and intensity of the trophic cascade to be studied. This contributes to a better understanding of the factors that determine the structuring of communities.

The article is structured in six sections. The second section deals with the presentation of the model. The third section focuses on the dynamic analysis of the model and the determination of coexistence conditions. The fourth section presents a type of invasion criterion and the trophic cascade intensity phenomenon. The fifth section contains the results of the coexistence analysis, energy requirements and trophic cascade intensity. Finally, the sixth section corresponds to a discussion of the results obtained.

We analyze a three-level system of equations of the Rosenzweig-MacArthur-type model. Consider a food chain in the following energy transfer sequence: Resource ↪ Consumer ↪ Predator, with respective body mass sizes m_R, m_C and m_P and, abundances functions in time represented by R(·), C(·) and P(·), respectively. We write the change in population densities as follows: dX/dt = G_X·X, withX∈{R, C, P}. Where G_X is the per capita growth rate and is described in detail in the following terms:

where L(mR,R)=r(mR)·{1-R/K(mR)} is the logistic type per capita resource rate. The parameters r(m_R) and K(m_R), also denoted by r and K, are considered proportional to mRαr-1 and mR-γ respectively (Damuth, 1993), where scale exponents α_r and γ are number in the interval [2/3, 1]. In addition, the parameter f_CR is linked to consumer predation of the resource but independent of their body sizes or abundances.

where parameters e_C (related to conversion efficiency) and δ_C (related to mortality) are not masses depend. Moreover, f_PC is associated with the top-predation of the consumer and does not depend on their body sizes or abundances. For example, the parameters f_CR and f_PC are linked to the behavior of the predator and the prey, the hunting strategies of the predator, the movement of the prey, the anti-predator behavior or the physiological stress that the predator exerts on the prey (Schmitz, 2017).

where e_C and e_P are not body sizes depend.

Equations (1–3) articulated in a three-level chain allow to address trophic interactions across different scales of body sizes. This approach contribute to the understanding of complex processes that take place in different ecosystems (Abdala-Roberts et al., 2019). For a better understanding of the results, let's consider the following direction of energy transfer from X to Y, i.e., X↪Y (R↪C or C↪P). The conversion efficiency of species Y (C or P) it is proportional to the quotient, m_X/m_Y, with e_Y as the constant of proportionality. In addition, the per capita mortality rate is proportional, with constant of proportionality δ_Y, to a mass-dependent term d(m_Y) equals to mYαd-1, with 2/3 ≤ α_d ≤ 1 (Loeuille and Loreau, 2005). These parametric relationships have been corroborated by theoretical and empirical investigations in the context of the metabolic theory of ecology (West et al., 2001; Brown et al., 2004; Savage et al., 2004; Hatton et al., 2019).

In a link X↪Y, the consumption is at a non-negative per capita rate proportional to ϕ(m_Y, m_X)·X called a functional response, that allow us to relate predation rate with the abundance of the consumer species, encapsulating physiological and behavioral aspects (Holling, 1959). In empirical research, it has been recorded that these functions also depend on the body sizes of the species (Weitz and Levin, 2006; Kalinkat et al., 2013; Ortiz and Arim, 2016; Pawar et al., 2019).

The parameter ϕ(m_Y, m_X) represents the mass-dependent component in the capture efficiency. It is constructed as the product between the probability, Π(m_Y/m_X), of successfully killing a prey of body size m_X consumed by a predator of body size m_Y (the probability increases as the predator has more body size) and I(m_Y, m_X) is the interaction rate per unit predator density (Weitz and Levin, 2006; Campillay-Llanos et al., 2021). So, it is

This product has been proposed by Weitz and Levin (2006) and described in detail in Campillay-Llanos et al. (2021). The interaction rate I(mY,mX)=mYβ·F(mY/mX), is a product where β represents the intensity of the predator mass in the interaction and the factor F(·), a positive value, is a scaling function that describes the behavior of the interaction rate relative to the quotient m_Y/m_X. In two extreme cases (domain boundary), F(·) has the necessary properties to ensure that: If m_Y≪m_X (read the ≪ sign as “much less than"), then we have I(mY,mX)~f0mXβ and on the other extreme, if m_X≪m_Y, we must have I(mY,mX)~f∞mYβ, where f₀ and f_∞ are positive constant. In this study, we consider F(ν)=f0/νβ+f∞ and Π(ν) = 1−e^−ν2, therefore the ϕ(m_Y, m_X) model that we will follow for graphics and simulations is:

On the other hand, in the absence of apex predators (P = 0), the model is expressed by Equations (1) and (2), which is the called Rosenzweig-MacArthur model (Rosenzweig and MacArthur, 1963), that has been proposed and studied by Weitz and Levin (2006). While, in the absence of intermediate consumers and apex predator (C = P = 0), the basal resource grows logistically, represented by dR/dt=L(mR,R)·R.

Table 1 provides a summary of variables and the parameters presented.

Let us introduce the reproductive number RY(X) as the average number of predators type Y that one predator leaves while it lives when X is the abundance of prey. Considering this conceptualization, we can mention that the population of apex predators Y (we refer to P or resp. C when P = 0) increases if RY(X)>1 , and decreases if RY(X)<1. More precisely, we have RY(X):=RY(1)·X, with

We also introduce RR(R,C):=L(mR,R)·[f CRϕ(mC,mR)C]-1, corresponds to the effective reproductive number of a resource unit when abundance is R and there are C consumers.

In order to have a better control over the signs of the derivatives (increase, balance or decrease) of the system represented by Equations (1–3) we perform algebraic operations to obtain an equivalent system, but in terms of the reproductive numbers:

where the state space for the variable (R, C, P) is the set Δ = [0, K(m_R)] × [0, ∞) × [0, ∞).

Our model is a straightforward representation of the community dynamics. However, it allowed us to investigate the conditions of coexistence characterized by the net reproductive rates, to study the energy requirements and the intensity of the trophic cascade. The empirical values of the allometric exponents play an essential role in our analysis. As in the existing literature, we assume that α_r = α_d = γ equals to 3/4 (Weitz and Levin, 2006; DeLong and Vasseur, 2012; Hatton et al., 2019) and β equal to 1/2 (Weitz and Levin, 2006). In our approach, we used the epidemiological notion of the basic reproductive number (Garrione and Rebelo, 2016) in the context of a three-level food chain. This notion, called net reproductive rate, represents the number of predators that a predator of average body size leaves during its lifetime. Coexistence and overall stability of the system are obtained by considering: RC(K)>1 and RP(C*)>1. These rates are not monotonic concerning the m_C/m_R and m_P/m_C ratios respectively. As a consequence, this morphological constraint is associated with the existence of feeding hierarchies (Arim et al., 2007). Moreover, as argued by Weitz and Levin (2006) using the invasion criterion (see Section 4.1), the process of substitution by invading predators of different body sizes through the mechanism of natural selection will converge to an optimal mass size (stable convergent strategy). Thus, we deduce that an invading predator cannot survive in the chain where the native predator settles with the same body size (stable evolutionarily strategy).

Under the condition α_r>β, our research proposes a space of species body sizes that can be explored in current or past natural environments as indicated in Figure 2 (Marshall et al., 2021). In this three dimensional space, the value of the allometric exponent α_r represents (a) α_r = 2/3, the surface law indicating the surface-volume constraints on heat dissipation over the surface of geometrically similar body planes, (b) α_r = 3/4, this value is known as Kleiber's law (Kleiber, 1932) and (c) α_r = 1, this value indicates a constant energy flux per unit of tissue mass and is called isometry. In this way, in Figures 2A,B, the mass space is similar to a parallelepiped. This space guarantees a variety of body sizes, indicating a greater diversity of species. Whereas in Figure 2C, it is observed that a more limited region of body masses is obtained for the first two trophic levels. Under these energetic assumptions, the coexistence region admits apex predators that exceed the size of one of the largest predators that have existed (average body size of Tyrannosaurus rex: 5,200 kilograms, Farlow, 1993). Therefore, our approach may serve to explore current trophic interactions and interactions of species that lived in other geological times.

In relation to demographic patterns, our model allows us to establish that the equilibrium density of intermediate consumers, C^*, is linked to body mass m_C with allometric exponent −[β−α_r+1] (In the absence of apex predators, a similar result is established for R_*). An analysis of a dataset reported by Damuth (1987) reveals that it is appropriate to consider the values of the allometric exponents α_r = 3/4 and β = 1/2 resulting in a negative relationship between density and body size at the second trophic level (herbivorous mammals) with allometric exponent –3/4 (Damuth, 1981). As a result, we obtain that the equilibrium abundance R^* is bounded below and above by body mass-dependent amounts (Damuth, 2001). It should be added that several studies have shown that densities of some species are not proportional to body size in −3/4 and even positive relationships between density and body size have been observed (Russo et al., 2003; Maxwell and Jennings, 2006; van Langevelde et al., 2020). Moreover, these relationships can be generated in a non-linear and polygonal fashion (Leaper and Raffaelli, 1999; Andrew and Hughes, 2008). However, Figure 3 shows that to keep the apex predator population in equilibrium, consumer (C) and resource (R) densities fall between a band that corresponds to a wide number of species.

The demographic patterns analyzed above, multiplied by the individual metabolic rate (Kleiber, 1932) allow estimating the energetic requirements, which represent the resources needed to maintain a population of a given density as proposed in macroecological theory (Brown, 1995). In the model, considering the prey-predator link X↪Y, the energetic requirements of a local population of species Y, denoted by E_Y, are mathematically presented by the product between the prey density X^* and individual metabolic rate B_X (EY=X*·BX). Setting the parameter values: α_r = α_d = γ = 3/4 and β = 1/2, it is possible to estimate the energy requirements of the predators, E_P, as shown in Figure 4. In the same way as Damuth (1981) has recorded, in our model and with the parameter values specified in the beginning, E_P is set to be independent of the body size of intermediate consumers (Damuth, 1981). This phenomenon is called the energy equivalence rule (Isaac et al., 2013). This implies that body size does not confer advantages in competition for energy use among populations (Isaac et al., 2013; Sewall et al., 2013). Recently, Hatton et al. (2019) have made it clear that the energy equivalence rule persists over a wide taxonomic and body size range. Furthermore, in our model it is possible to estimate the energy requirements of intermediate consumers. As shown in the Figure 4, we obtain that E_C is bounded by a lower and upper level of energy resources that does not depend on the body sizes of the basal resource (Damuth, 1987; Marquet et al., 1995; Charnov et al., 2001). The limiting of basal resource densities is linked to alternative assumptions to the energy equivalence rule that energy requirements will increase with body size, decrease with body size or reach a maximum level for a given body size (Sewall et al., 2013). This is consistent with Silva et al. (1997) who argue that species from different trophic groups exhibit different allometric exponents. The mathematical expressions that have been used to perform the analysis of the demographic patterns and energy requirements (Y = b·m^α), are also called Power Laws (Marquet, 2000; Schneider, 2001). These mathematical relationships describe the behavior of nonlinear interactions, where one dependent variable is expressed as the power of an independent variable and where the exponent value determine the scaling relationship. Power laws are scale invariant, a property that contribute to the understanding of general principles of ecological systems across different levels of organization (Marquet et al., 2005).

On the other hand, predators affect directly the abundance of their prey but this effect can also extend to lower trophic levels indirectly. Considering a three-level chain, the apex predator reduces the abundance of the second level, which in turn allows higher abundances at the lower trophic level. The strength of this indirect interaction measured in terms of abundance differences at equilibrium is called trophic cascade intensity. Piovia-Scott et al. (2017) argue that trophic cascade intensity depends on several species and ecosystem factors. Moreover, multiple models and metrics are used to quantify the intensity of trophic cascades (Leroux and Loreau, 2008). However, it is reasonable to assume that intensity depends on the species' body size (Delong et al., 2015). Our mathematical representation of trophic cascade intensity corresponds to a positive function with varying species body sizes. In agreement, with the coexistence conditions, we require that R*<R* so that T_C is greater than 1. In addition, the condition that the energy requirements of apex predators (E_P) are limited is necessary for this model to make ecological sense (i.e., that it is a positive function). This is because the model specifies that consumers are specialists and does not incorporate external food flows. These factors are fundamental to determining the intensity of the cascade and can be modified by incorporating food subsidies at the different trophic levels (Leroux and Loreau, 2008). On the other hand, by setting the mass of intermediate consumers at m_C = 10, we hypothesize that in communities larger predators can generate stronger trophic cascade intensities than smaller predators, shown in Figure 5. Our results are consistent with empirical research (Brose et al., 2006a; Terborgh et al., 2010) and is reported in Delong et al. (2015). Furthermore, our theoretical contribution by characterizing trophic cascade intensity with the net reproductive rates may eventually support current research on possible mechanisms to explain variation in trophic cascade intensity.

Our research is not without limitations, as our models do not incorporate spatial and/or environmental factors that may contribute to maintaining or modifying demographic patterns. Despite this, we propose that body sizes influence the strength, distribution, and characteristics of interactions: delimit species body size ranges, which determine coexistence, demographic patterns, energetic requirements and the intensity of the trophic cascade inherent to the three-species food chain. It should also be mentioned that body size is a climate-sensitive trait (Binzer et al., 2012; Ohlberger, 2013; McLean et al., 2020). Challenging these complications, we would like to integrate the knowledge of how different types of organisms function in their physical, chemical, and biological environments into these types of ecological-mathematical models, to understand how global changes will affect organisms and biodiversity. Potentially, our research could be used to guide future empirical studies from a macroecological point of view.

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

WC-L carried out both the mathematical model and analysis of this and wrote the first draft of the manuscript. All authors contributed to the conceptualization of the study, wrote, reviewed, and edited the manuscript.

This research was supported by the Doctorado en Modelamiento Matemático Aplicado (DM₂A-UCM), a program belonging to the Universidad Católica del Maule, Talca, Chile.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.