Predator–prey power laws: trophic interactions give rise to scale-invariant ecosystems

Scaling laws and power-law distributions are ubiquitous in ecological systems. However, it is not clear what factors give rise to such universal regularities. Here, I show scaling laws are a simple consequence of scale-invariant distributions, and both result from simple commonalities of diverse ecosystems. I introduce a simple model of predator–prey interactions in which predators and prey move on a two-dimensional space in search of resources that they use to survive and reproduce. As primary resources increase, the food web exhibits a series of transitions to phases with equilibrium dynamics and top-down control of the food web, non-equilibrium dynamics with bottom-up control, and unstable dynamics exhibiting the paradox of enrichment. The model shows resource heterogeneity can solve the paradox of enrichment and ensure the stability of ecosystems. Scale-invariant spatial distribution of prey and predators and a surprisingly rich set of scaling laws, including predator–prey and Taylor’s power laws, appear in the non-equilibrium phase. The model predicts both Taylor’s power law and predator–prey power law can be extended to a rich set of fluctuation scaling laws governing the fluctuation of predator’s and prey’s densities and growth. A mathematical theory suggests scaling laws result from the scale-invariance of the spatial distribution of prey and predators.


Introduction
Scaling laws and power laws are ubiquitous in many biological and ecological systems [1][2][3].These laws arise across different levels, from the individual level, such as metabolic scaling laws pertaining to individuals' attributes [3,4], to the population level, such as Taylor's power law relating populations' fluctuation and abundance [5][6][7], or even in ecosystems and evolutionary level [8][9][10][11].Among the scaling laws present at the ecosystem level, recently, a curious scaling law according to which both predator abundance and prey production show a sub-linear power-law dependence on prey abundance with an exponent close to 3/4 is argued to be at work across different ecosystems from terrestrial to marine communities [12,13].Given the stark differences in composition and underlying interactions in these ecosystems, the emergence of such universal patterns may suggest some underlying scale-invariant properties are at work in ecosystems that originate from simple commonalities shared by diverse ecosystems.Indeed signatures of scale-invariance, such as scale-invariant group size distributions, have been observed in many biological populations [14][15][16][17][18]. Yet, it is not clear how such distributions relate to scaling laws and what factors give rise to these scale-invariant properties.To address this question, here I introduce a simple agent-based model of predator-prey interaction in which preys and predators move on a two-dimensional space in search of resources, which they use to pay for their metabolic costs and reproduce.Preys use primary resources, and predators reach resources by catching prey.
As the primary resources increase, the model shows different phases in which predators do not survive, prey and predators coexist in equilibrium, and the system exhibits top-down control of the food web, in which predators control the food web and play a more important role in deriving the systems dynamics, and a non-equilibrium state in which prey and predators exhibit highly aggregated dynamics in the form of traveling waves and the food-web shows bottom-up control, and the food web is derived by prey population.These observations shed theoretical light on empirical arguments for the existence of such different phases in predator-prey systems [19][20][21], and the conditions under which they may show top-down or bottom-up control [22][23][24].The predator-prey dynamics become unstable in too-rich environments, resulting from the too-fast exponential growth of predators acting faster than the regulatory feedback from prey to predator abundance.The instability of the predator-prey dynamics in resource-rich environments coined the paradox of enrichment, has been argued theoretically [25,26], and observed experimentally [19,20].As the model shows, in contrast to what has been suggested before [27], although spatial structure increases stability, it fails to restore the stability of the system fully.However, as I show here, resource heterogeneity-a nonuniform distribution of resources in space-restores the stability of the system.
A rich set of scaling laws emerges in large densities, where the system exhibits non-equilibrium behavior and bottom-up control.Several variables, including prey and predator densities, density fluctuation, birth and death rates, and the maximum group size of prey and predators, obey a power-law relation with the resource regeneration rate.Furthermore, the temporal fluctuation of these variables shows scaling with primary resources.The primary scaling laws imply any of these variables obey a power-law relation with respect to each other.This finding gives rise to a surprisingly rich set of scaling laws, including several fluctuation scaling laws, and suggests the scaling laws empirically discovered to hold in ecosystems are just the tip of an iceberg of scale-invariant ecosystems.For instance, predator and prey spatial density fluctuations scale as a function of each other with an exponent close to 1 in large densities, which recovers Taylor's law in space with highly aggregated dynamics [5].In addition, predator density and prey production obey a scaling relation as a function of prey density with a sublinear exponent, which recovers the predator-prey power laws [12,13].Furthermore, the model predicts the spatial distribution of prey and predators, which can be considered as group size distribution obeys a power law with a sharp cut-off.The cut-off, the largest group size, shows a scaling law with respect to resource regeneration rate.I will develop a general theory extending beyond predator-prey interactions to argue that these two conditions can give rise to a scale-invariant ecosystem exhibiting power-law dependence on population densities and highly aggregated dynamics.Given that the scale-invariance of group size is argued to hold in many biological populations [14][15][16][17][18], the theory developed here suggests similar scaling laws may be at work in other contexts where populations coexist and interact.
The model also sheds light on one of the most fundamental questions in ecology, the direction of control in food webs [22][23][24], by establishing a link between the nature of control, prey and predator life histories, their collective dynamics, and behavioral responses to each other.In large densities with bottom-up control, prey's lifespan increases, and their per-capita reproduction decreases with density, while predators' lifespan decreases and their per-capita reproduction slightly increases with density.As the birth rate of predators and prey scales almost linearly, a higher lifespan of prey in rich and dense environments implies a sublinear scaling of predator abundance with respect to prey abundance.This observation suggests increasing the resources, and thus increasing the population, favors prey and is consistent with the fact that the prey population derives the ecosystem in a bottom-up controlled food web.On the other hand, in small densities with top-down control, predators' lifespan shows little density dependence, and their per-capita reproduction increases with density.Prey's lifespan decreases, and their per-capita reproduction increases with density.This suggests increasing the density favors predators, and it is consistent with the fact that in a top-down controlled food web, the dynamic is governed by predator population.These observations link the top-down or bottom-up nature of control in food webs with the life histories of prey and predators.
I will argue that a sublinear scaling in large densities and the nature of the control in the food web results from behavioral responses of prey and predators to each other.In large densities, a shelter effect-the benefit of living in large groups which can result from preys' or predators' behavioral responses to each other-plays a role, and the system exhibits highly aggregated non-equilibrium dynamics.In contrast to a simple mass-action law underlying Lotka-Volterra models [28][29][30][31], the shelter effect makes the resource flow to predators non-invariant to, and a decreasing function of density fluctuations.Due to higher density fluctuations in higher densities, resulting from Taylor's law with highly aggregated dynamics, prey per capita production in large densities decreases.This leads to a longer lifespan and lower population turnover of prey in large densities and a sublinear scaling of predator-prey densities.On the other hand, in small densities, the shelter effect does not play a role, mass-action law holds, preys exhibit weak density dependence, and the system shows top-down control.These findings point toward the importance of preys' and predators' behavioral responses to each other in determining the nature of the control and dynamics of predator-prey interactions.

The model
The model and its phenomenology are illustrated in figure 1.In the model, prey and predators move on a grid of linear size L.Each agent α, prey or predator, has an internal resource, ϵ α , which increases by gaining resources, κ α , and decreases by a constant rate, η, due to paying for metabolic costs.Thus, the internal resource of agent α at time t evolves according to: Preys acquire resources by consuming primary resources, which are regenerated on each site r, with a rate λ(r).Here λ is the primary resource regeneration rate.Basal resources available on a site are divided among all the prey visiting a site.Predators capture prey to gain resources.When a predator visits a site where prey lives, the predator captures a prey with probability q.When capturing prey, the internal resource of the predator increases by an amount γ, which can be considered as the resource gain by predators per prey or predator conversion efficiency.Preys and predators die if their internal resource reaches zero and reproduce if their internal resource reaches a threshold d.When reproducing, the individual produces an offspring, and its resources are divided between the individual and its offspring.After reproduction, prey and predators move to a randomly chosen neighboring site on the lattice.If it happens that one of the species goes to extinction, assuming there is migration, one individual of that species is added to a random location on the lattice.I will consider both a homogeneous environment, in which the resources on all the sites are regenerated with a constant rate λ, and a heterogeneous environment with a linear distribution of resources in which the resource regeneration rate λ(r) on a site r = (x, y) is given by λ(r) = 2λx/L.The case of binary resource distribution is considered in the supplemental material, and shows a similar phenomenology (see supplementary note 7).For the spatial structure, both a nearest neighbor square lattice with linear size L, fixed boundaries and von Neumann connectivity (structured population) and a fully connected network (mixed population) are considered.
In the following, I will show densities (number of individuals per site by ρ).The birth rate is shown using the symbol b, and the death rate is shown by d.Prey production, defined as the number of prey caught by predators per time step, is shown by P p .I use the symbols K and K 0 to refer to an analytical expression for prey production and its approximation, respectively, derived in the section 5. δρ p and δρ π show the spatial density fluctuations of prey and predators, respectively, defined as the spatial standard deviation of densities, δρ = √ ⟨s(r) − ρ) 2 ⟩, where s(r) shows the number of individuals (prey or predators) on a lattice site and ρ = ⟨s⟩ is the average density.⟨.⟩ indicates an average, and ⟨.⟩ t indicates a time average.A subindex p refers to prey and π to predators.We use a super index c to refer to per capita variables.A per capita variable, such as per capita birth rate, is the birth rate per individual per timestep and is calculated by dividing the per area variable by density, b c p = b p /ρ p .A super index ' * ' on a variable shows the value of the corresponding variable at the transition point (see below).

Results
I begin by showing that the model predicts predator-prey systems show different phases in section 3.1.Furthermore, I show that while unstable predator-prey dynamics result in too rich environments in a homogeneous environment, resource heterogeneity restores the stability of the ecosystem.In section 3.2, I show that the model exhibits scaling laws for the mean and fluctuation of several variables.These laws can explain and extend the predator-prey power law and Taylor's power laws.In section 3.3, I provide a mathematical theory to argue that scaling laws result from the scale-invariant spatial distribution of individuals.Finally, in section 3.4, I discuss how the model sheds light on predator and prey life histories and relate food web control with the life histories of prey and predators.

The phase diagram and phase transitions in predator-prey systems
The phase diagram of the model in the λ-γ plane and the densities of preys and predators as a function of λ for γ = 0.5 are plotted in, respectively, figures 2(a) and (b).For too low resource regeneration rates, predators do not survive, and prey density is determined by equating resource production λ and resource consumption ρ p η, ρ p = λ/η.By increasing λ, a phase transition above which predators survive occurs (blue circle).A naive mean-field argument predicts the transition occurs at λ * = η 2 /γq (by equating resource flow to predators, γqρ p , and predator resource consumption η, and using ρ p = λ/η, see Methods), in good agreement with The illustration of the model and its phenomenology.(a) Prey and predators move on a lattice.Preys use primary resources, regenerated by a rate λ per site, and predators capture prey (successfully with probability q) to reach an amount of γ of resources.Individuals pay a metabolic cost of η per time step, reproduce when their energy reaches a threshold d, and die out if they go out of energy.(b) The model predicts a top-heavy food web in small densities, where mass action law holds, a weak density dependence follows, predator-favoring life histories emerge, and the food web shows top-down control and higher resource leads to a higher density of predators.On the other hand, in large densities, the shelter effect gives rise to a scale-invariant ecosystem with a scale-invariant group size distribution and scaling laws, bottom-heavy food web, bottom-up control, and prey-favoring life histories.simulation results (figure 2(a)).The predator-prey phase transition is a continuous phase transition and exhibits critical scaling [32].Close to the phase transition, predator density, prey and predator birth rates, and prey production, when subtracted from their value at the critical point, show critical scaling with respect to the reduced resource regeneration rate, λ − λ * , where λ * is the critical value of the resource regeneration rate.These critical scaling relations are presented in figure 2(c).We note that the critical exponents of all these variables are universal and are the same for different parameter values.
Above the predator-prey transition, the system reaches an equilibrium state where preys and predators coexist, and no collective movement emerges (see the supplementary videos and supplementary notes 3 and 4 for the dynamics of the model).In this region, the density of prey shows little variation and slightly decreases with increasing λ.Instead, higher resources available in the system are transferred to predators, and their density increases by increasing λ.A simple mean-field argument predicts the density of preys in the equilibrium to be equal to η/(γq), with a fair agreement with simulation results (see the supplementary note 1).
By increasing λ beyond a second threshold (red squares), non-equilibrium fluctuations are set in a structured population.These fluctuations are derived when the density of predators becomes high enough such that the speed with which predators consume preys is higher than the speed with which preys arrive in that neighborhood.As a result, the prey's density field becomes inhomogeneous, and the dynamic is derived by traveling waves of prey followed by fronts of predators (see the supplementary videos and the supplementary notes 3 and 4).A simple argument suggests this phenomenon occurs when predators become abundant enough to be able to deplete the prey density field locally (see the supplementary note 1).
Above a third threshold (orange diamond), the density fluctuations increase, and the dynamics become unstable due to the fact that high prey abundance elicits a too-fast exponential growth of predators, which acts faster than the regulatory feedback from prey to predator density.Adding prey or predators into the system after the extinction of their population leads to cycles through which the density of the prey population increases when there are low or no predators, leading to the growth of predators.This, in turn, results in the decline and eventual extinction of prey.The dynamics of the system thus go through cycles of coexistence of prey and predators and their extinction.The time average of the densities of prey and predators in this regime in a homogeneous environment shows a linear scaling with λ.For too low resource regeneration rates, only prey survive (marked by the letter (p).By increasing λ, the system shows a series of successive phase transitions to a phase where preys and predators coexist in equilibrium, E, a non-equilibrium phase marked by traveling waves of preys and predators, N, and a phase where the dynamics become unstable due to too large density fluctuations resulting in predators driving preys to extinction, U.The dashed line, λ * = η 2 /γq, where η is the resource consumption per time step, and q is the probability of successful predation, is a naive mean-field prediction for the phase transition above which predators survive, which is in good agreement with simulation results in a mixed population (black stars).(b) The densities of prey, ρp, and predators ρπ, as a function of λ for γ = 0.5.Different phases are shown in the system by letters.For too high resource regeneration rates, the system is unstable, and the time average of prey and predator densities in the presence of immigration shows a linear relation as a function of λ in a homogeneous environment.(c) The predator-prey transition, above which predators survive is a critical phase transition.Consequently, the reduced predator density, ρπ − ρ * π , reduced prey birth rate bp − b * p , reduced predator birth rate, bπ − b * π , and reduced prey production, Pp − P * p , show scaling as a function of reduced resource regeneration rate, λ − λ * close to the phase transition.Here, a superindex c indicates the value of the corresponding variable at the transition point.A reduced variable is defined as the value of the variable subtracted by its value at the transition.Prey production is defined as the number of prey captured by predators per time step.Parameter values: η = 0.05, d = 2, γ = 0.5, q = 0.5.A homogeneous environment is considered.Error bars are standard deviation (divided by 2 in (b) to increase visibility).
The instability of the predator-prey dynamic in resource-rich environments, coined the paradox of enrichment, is argued to result in simple models of predator-prey dynamics [26] and has found some empirical support [19,20,25].However, it is argued natural predator-prey systems can escape the paradox of enrichment using different mechanisms [25].For instance, using theoretical models, it is argued that spatial structure can stabilize predator-prey dynamics [27].The model suggests that while spatial structure can improve the stability of the system by shifting the instability to richer environments, the instability of the dynamics endures in a spatial structure (see the supplementary note 2).However, resource heterogeneity provides a simple mechanism to restore the stability of the dynamics due to the fact that in a heterogeneous environment, resource-poor regions can provide a niche for preys in which they can survive due to the low abundance of predators and from time to time immigrate to resource-rich environments (see the supplementary note 2).

Scaling laws in predator-prey systems
By stabilizing the predator-prey dynamics, resource heterogeneity allows inspection of scaling in the stable regime and reveals a rich set of mean and fluctuation scaling laws.Several variables, including prey density, predator density, preys' density fluctuation, predators' density fluctuation, preys birth rate, predator's birth rate, and prey production, show scaling with primary resources, λ, in large densities (see the supplementary note 6).Furthermore, the temporal fluctuation of these variables, defined as their standard deviation in time, exhibits a scaling relation with primary resources (see the supplementary note 6).The existence of scaling with primary resources suggests these variables show scaling as a function of each other (as for two generic variables, v and w, with scaling relations, v ∝ λ a λ,v and w ∝ λ a λ,w , we have v ∝ w a λ,v /a λ,w ).While many of the scaling laws predicted by the model, such as the scaling of the birth rates of prey and predators or primary scaling relations, are not yet empirically discovered in ecosystems, some of the consequences of the scaling relations are empirically observed but not yet theoretically understood.An example is the predator-prey Figure 3. Predator-prey power laws.Predator density, ρπ, as a function of prey density, ρp (a), prey production, Pp, defined as the number of prey caught by predators, as a function of prey density, ρp, (b), for three different system sizes, L, and a linear distribution of resources are plotted.The system shows a sublinear predator-prey and prey production-prey density scaling in large densities in both stable and unstable regimes.(c) and (d) show prey and predator spatial density fluctuation, δρp and δρπ, respectively, as a function of their density.The system is unstable for large densities for small system sizes (L = 25 and L = 50).Both prey and predators show Taylor's power law with an exponent close to 0.5 for small densities, indicating weak density dependence, and 1 for large densities, indicating a highly aggregated dynamic.Parameter values: η = 0.05, d = 2, γ = 0.5, q = 0.5.Error bars are standard deviation divided by 5 to increase visibility.power law: in figures 3(a) and (b), I plot the density of predators and prey production as a function of prey density for three different system sizes.As the system size increases, the dynamics stabilize up to larger values of λ.For L = 25, both predator density and prey production show sublinear scaling with respect to prey density.For the unstable dynamics, the scaling is closer to linear (with a slope of ≈0.86), and in the stable regime, the scaling is close to a three-quarter scaling observed in predator-prey systems [12].A sublinear scaling holds consistently for other parameter values and resource distribution (see the supplementary notes 6 and 7).Another consequence of scaling relations is Taylor's power law.Both prey and predator spatial fluctuations show Taylor's power law with an exponent of 0.5 in the equilibrium regime, indicating weak density dependence, and an exponent of 1 in the non-equilibrium regime, indicating highly aggregated dynamics (figure 3(c) and (d)).
The model also exhibits a rich set of temporal fluctuation scaling relations.The temporal fluctuation of densities exhibits Taylor's power law in time, with an exponent which increases going to higher levels of the food web (4(a)).Taylor's power law in time also holds for other quantities: both the birth rate and spatial density fluctuations show fluctuation scaling (see 4(b) for fluctuation scaling of birth rates).In small densities, an exponent of 0.5 is observed.This exponent can result when the system lacks strong correlations, and birth rates can be regarded as a sum of independent variables.On the other hand, in the scale-invariant regime, a super-linear exponent is observed, which signals the emergence of strong dependencies and synchrony between the elements of the system [7], as manifested by the onset of traveling waves.The fact that the mean and fluctuation of different variables scale with primary production in the ecosystem implies Taylor's power laws can be extended to heterogeneous power laws, which relate the fluctuations of a variable to the mean of another variable.Two examples of these power laws, presented in figure 4(c), are the scaling of the temporal fluctuation of prey production as a function of prey density (left) and temporal fluctuation of prey birth rate as a function of their density (left).Finally, primary fluctuation scaling relations imply temporal fluctuations of variables pertaining to different species also scale with each other.This observation implies predator-prey fluctuation power laws, which can be considered as a corollary to the predator-prey power law, according to which the temporal fluctuations of both predator and prey production show a power law as a function of prey temporal fluctuation (4(d)).Examination of other parameter values suggests predator-prey fluctuation power law is a sublinear power law.However, it shows a larger exponent compared to the predator-prey power law.Furthermore, in contrast to the predator-prey power laws which shows the same predator-prey density exponent, a ρp,ρπ , and prey production-prey density exponent, a ρp,Pp , predator-prey fluctuation power law shows a slightly larger predator-prey temporal density fluctuation exponent, a σρ p ,σρ π , than prey production-prey density fluctuation exponent, a σρ p ,σP p , (see the supplementary note 6).

Scale-invariant spatial distribution and the origin of scaling laws
An important question is the origin of scaling.In the following, I argue that the scaling of averages can be understood in terms of the scale-invariance of the ecosystems.As can be seen in figures 5(a) and (b), the spatial or the group size distribution of both prey and predator, s x , defined as the number of prey or predators observed in a patch, for different resource regeneration rates, shows scale-invariance with a lower cut-off for small group sizes, s 0 x , and an upper cut-off for large group sizes, s F x , where x stands for prey, p, or predator π.Postulating a simple scaling form f(s x , λ) ∼ λ −β f 0 (s x λ −α ), the distribution for different λ collapse into a universal power-law f 0 (the insets in figures 5(a) and (b)).Besides, the upper cut-off, which is the largest group size, obeys a simple power-law relation with respect to λ, with an exponent, ν x (figure 5(c) and (d)).In the following, I show that the scale-invariance of the group size distribution, together with the scaling of the cut-off, leads to the scaling of the species population densities (see Methods for details).Noting that the density of a species x can be written as ρ x (λ) = ´∞ 0 s x f(s x , λ)ds, using the simple scaling form, taking the lower cut-off to be equal to zero, approximating the infinity in the upper bound of integral by s F x = λ νx , and assuming f 0 (s x ) has a power law form, f 0 (s x ) = s −δx x , it is easy to see: Thus, a λ,ρx = −β x + α x δ x + (2 − δ x )ν x .Using similar calculations, it is possible to show that the exponent of the nth moment of population density of species x, ρ Finally, noting that in the stationary state the prey production should balance predator consumption, K(λ) = ηρ π (λ), we get the same exponent for prey production, a λ,K = a λ,ρπ = −β p + α p δ p + (2 − δ p )ν p .
Finally, using the normalization of the probability density function, it is possible to derive the approximation −β x + α x δ x = −(1 − δ x )ν x , which gives a λ,ρx = ν x .This predicts a Taylor's power-law with an exponent of one, δρ x = (ρ x − ρ 2 x ) 1/2 ∼ ρ x ∼ λ νx , that is a highly aggregated dynamics observed in the non-equilibrium regime.
Although the model exhibits a sub-linear predator-prey power law, the preceding mathematical theory regarding the origin of scaling is consistent with both superlinear and sublinear scaling.A curious question is, what mechanism drives the sublinear scaling of prey and predator densities?It is possible to derive an analytical form for prey production based on their spatial distribution.We denote the analytical expression for prey production by K, which can be approximated by K 0 = q⟨min(s p , s π )⟩, and is calculated in the section 5, equations ( 4) and ( 5), respectively.Here, s p and s π are the numbers of prey and predators on a lattice site.Underlying these functional forms of prey production is a key ingredient of the model: predators get only one chance to attempt predation.Realistically, this feature can result from the prey's behavioral responses to an unsuccessful predation attempt as prey can better avoid predation following an unsuccessful predation attempt, for instance, by scattering around [33,34], showing a higher vigilance [33][34][35], alarming the peers [33,36,37], or due to predator's fatigue, and can be called the shelter effect.Importantly, this functional form is in stark contrast with the mass-action law functional response typically considered in Lotka-Volterra models [28][29][30][31], which can be written as the correlation function of prey and predators densities, ⟨s p s π ⟩, times the probability of successful predation, q.In figure 6(b), I plot K and K 0 together with a mass action law functional response q⟨s p s π ⟩.I note that these quantities are calculated based on the spatial distribution of prey and predators derived empirically from simulations (that is, the average involved in the calculation of these quantities is performed using the empirical distribution of prey and predators observed in simulations).A mass action law functional response is a valid approximation in low densities and correctly predicts an invariant prey density.This can be seen by noting that in equilibrium, resource gain per predator, which according to the mass action law is qγρ p , should balance resource consumption, η.This gives, ρ p = η/(γq).Thus, mass action law holds in small densities where density fluctuations can be ignored and predicts an invariant prey density with respect to resource regeneration rate.For larger densities, the shelter effect plays a role, and mass action law no longer holds.The minimization procedure involved in K 0 , or the shelter effect reflected in K, makes prey production sensitive to density fluctuations due to the fact that lower p , and natural death rate, b c p − P c p , defined as the death rate of preys due to going out of resource, as a function of prey density.K exactly describes prey production.K and K0 are derived in the methods, equations ( 4) and ( 5), and are the analytical expression for prey production, and an approximation for prey production, respectively.(b) While K and K0 show scaling with prey density with the same exponent as prey production, the correlation function between prey and predator densities times q, q⟨spsπ⟩, which determines prey production if mass action law holds, deviates from prey production in the non-equilibrium regime.The population resides on a first nearest neighbor lattice with von Neumann connectivity and fixed boundaries, with linear size L = 200.Parameter values: η = 0.05, d = 2, γ = 0.5, q = 0.5.resource flow in places where prey are found in a small density is not compensated by a higher resource flow in prey rich places.As density fluctuation increases with density due to Taylor's law, prey production decreases with density.This leads to a sublinear scaling of prey production and predator density.We note that the deviation of predation functional response from a mass action law and its sub-linearity, suggested by the model, has been recently argued for empirically [38].

Life histories and food web control
The model also sheds light on life history strategies (see the supplementary note 5).While in small densities most of the prey die fast and increasing the resources decreases the mean lifespan of prey and increases their reproduction, in large densities, the prey lifespan distribution develops a peak for finite lifespan due to the shelter effect.As a result, most of the prey can live for a finite period of time and increasing the density increases the mean lifespan of prey and decreases their per-capita birth rate, while predator's mean lifespan decreases with increasing density.The difference in the lifespan of prey and predators, in turn, leads to the sublinear scaling of predator density as a function of prey density: as prey and predator per area birth rates scale almost linearly as a function of each other, longer lifespan of preys by increasing density leads to higher prey density in the system.These observations suggests that food-webs show different causality and control directions in small and large densities: in small densities, where equilibrium holds, food-webs exhibit top-down control and are controlled by predators, and in large densities, they exhibit bottom-up control and are derived mainly by prey.An examination of direction of causality based on flow of information [39] between prey and predator densities shows this is indeed the case (see the supplementary note 5).This finding suggests that the nature of the control in the food web results from prey and predators behavioral responses to each other.In small densities, predators are the winners of competition in this regime and control the food web: preys get caught rapidly and higher resources leads to higher predation, higher predator density, and shorter lifespan of preys.In large densities, the shelter effect hands the control of the food web to preys, and the benefit of living in large groups give rise to the emergence of a characteristic lifespan of preys, which is further strengthened by increasing the density.

Discussion
As we have seen here, a simple model of predator-prey interactions gives rise to a scale-invariant ecosystem, with a scale-invariant spatial distribution of individuals, which can be thought of as group size distribution and a rich set of scaling laws governing means or fluctuations of ecological variables or relating means to fluctuations.These scaling laws recover some of the empirically observed scaling laws and suggest these power laws are just the tip of an iceberg of scale-invariant ecosystems.The model suggests that the predator-prey power law and Taylor's power law can be explained based on a set of primary scaling laws and can be extended to scaling laws governing several ecological variables and their fluctuations.By developing a general mathematical theory, this study suggests scaling laws can result from the scale-invariance of the spatial distribution of individuals.Notably, the analytical theory relating the scale-invariance of spatial distributions and the scaling laws developed here is rather general.Its only assumption is that the spatial distribution obeys a power-law relation with a cut-off proportional to the power of resource availability.These assumptions naturally emerge in the agent-based model introduced here and appear to hold in many species and ecosystems.For instance, it is argued that group size distributions are truncated power laws [14,15,17,18], or power-laws with exponential decay, and the universality of these distributions are argued [16].In addition, an increasing relation between density and group size is often observed in animal groups [40][41][42].As we have argued using a general mathematical argument, these observations can lead to a rich set of scaling laws as long as different species coexist and interact in an ecosystem.Importantly, this argument is independent of the nature of the interactions and extends well beyond predator-prey interactions.
The model also reveals that trophic interactions can result in different physical phases, and sheds light on the phase transitions in these systems and how these phases can relate to ecosystems' function and structure.While the observation of different trophic phases and phase transitions may be difficult in the field due to factors such as immigration, or environmental noise, the existence of different phases in predator-prey systems depending on resource availability are observed in controlled lab experiments [19][20][21].A curious question in this regard which can be subject to future empirical work, is whether a critical predator-prey phase transition and universal critical scaling in predator-prey systems can be observed empirically.
The existence of different phases in food webs also has important implications for one of the most fundamental questions in ecology: the question of whether food-webs exhibit top-down or bottom-up control.[22][23][24].The model shows a top-down control results in small densities where the system shows equilibrium and non-aggregated dynamics, and bottom-up control results in large densities with the onset of collective motion of prey and predators with highly aggregated dynamics, where the shelter effect starts to play a role, and a characteristic lifespan of preys emerges.This finding shows prey's collective responses to predators have important consequences for fundamental ecological patterns and determine the nature of control in food webs.This observation shows the importance of going beyond simple Lotka-Volterra models by taking the behavioral responses of prey and predators into account [43] and how these behavioral factors can underlie fundamental ecological patterns [44].
While we have argued that ecological scaling laws can be explained in terms of the scale-invariance of trophic interactions, an interesting question is, what factors give rise to the scale-invariance of trophic interactions in the first place?Similar power law distributions have been observed in different contexts, from group size distribution in animals [14][15][16][17][18] to city size distribution [45][46][47].It is generally accepted that power laws hint at fundamental and universal principles governing systems' dynamics and, as such, may admit universal explanations, even when they appear in apparently different systems [45].The origin of power law distributions in biological and social systems has been subject to much research [45][46][47].Yule process, or 'rich gets richer effect' , also known as the Matthew effect and self-organized criticality, are two of the most appealed explanations of power laws in diverse systems [45,47,48].
A curious observation in this regard is the similarity of the predator-prey model introduced here and self-organized critical models [48,49].The dynamic of the predator-prey system is derived by the constant flow of energy into the system in the form of primary production and its dissipation along the food web gradient.Self-organized critical models have shown that such a flow of energy in dissipative physical systems can lead to the self-organization of the system in a critical state marked by the emergence of scale-invariant distributions and fractal structures [48,49].Similarly to such systems, ecological systems are derived by the constant flow of free energy into the system from the bottom-most level, and its dissipation through trophic interactions as it goes up the food web [50].Based on this observation, recent decades have witnessed the emergence of new insights according to which flow and dissipation of energy in ecosystems can provide a basis for a new ecological theory in which ecosystems can be regarded as self-organizing complex systems [51][52][53][54].The present study, by showing that a similar phenomenology to that observed in self-organized critical dissipative physical systems can hold in trophic systems and explains fundamental and ubiquitous ecological scaling laws, can take a step toward the conception of such an ecological theory.Future studies, by elaborating on the connection between trophic interactions and self-organized criticality, can take essential steps in this regard.
The research described here provides valuable insights that can be relevant to the field of social physics, particularly in the context of city size distribution and the Mathew effect.It discusses how scaling laws and scale-invariant distributions observed in ecological systems are analogous to patterns found in social and urban contexts.By revealing that these power laws are related to the scale-invariant spatial distribution of individuals, the study highlights the potential application of these ecological principles to understand phenomena like city size distribution and the Mathew effect in social systems.This suggests that the dynamics and laws governing ecological systems may offer valuable perspectives and analogies for understanding social phenomena, particularly those related to the distribution of resources and power.
Finally, the ecological model introduced here can provide a simple framework for evolution by allowing individuals to possess inheritable traits that determine their capability to obtain resources.As it will be shown elsewhere, such an evolutionary generalization of the model can provide a simple and powerful ecological account of evolution that can shed light on other fundamental ecological and evolutionary laws.Besides, while we have argued a sublinear prey production and predator-prey power law can result from the shelter effect, an interesting question is how other factors, such as cooperative defending, cooperative hunting, limited mobility, prey, and predator competition, or higher level food webs affect prey and predator growth and the nature of the predator-prey power law.While the model introduced here does not consider such complications for the sake of simplicity, it can be generalized to study these effects in future studies.Furthermore, as it is recently shown in the context of evolutionary games [55], such an evolutionary generalization of the model can provide a novel framework to study evolutionary games, which provides essential insights into the evolution of cooperation and life history traits.

The simulations
At each time step of the simulation, at first, preys consume resources.For this purpose, the resources available on each site are equally divided among the prey visiting that site.Then, predators are given a chance to capture prey.For this purpose, each predator on the site is given a chance to attempt to capture a prey.The predator is successful with probability q.This is continued until either all the prey on the site are captured or all the predators on the site attempt capturing a prey.In the next stage, preys and predators reproduce and move.For this purpose, those with an internal resource above d reproduce an offspring, and the resources are divided equally among the agent and its offspring.After reproduction, agents move to one of the neighboring sites chosen at random.Agents pay for the metabolic cost, η, and those with internal resources below 0 die out.Finally, the resources available on each site r increases by a rate λ(r).If it happens that the population of one of the species goes to zero, assuming there is migration, one individual of that species with internal resource 1 is added to a random location on the lattice.I note that immigration occurs only in the unstable regime (where extinction occurs) and does not occur, and thus, does not affect the dynamics of the system in the stable regime where no extinction occurs.At the beginning of the simulations, the individuals are distributed in random locations on the lattice.After checking that initial conditions does not affect the equilibrium dynamics, all the simulations start with 100 prey, 10 predators, and initial energy equal to 1.

Resource distribution
In the case of homogeneous resource distribution, the primary resource regeneration rate is taken to be a site-independent constant λ(r) = λ.For a linear resource regeneration rate, the primary resource distribution is given by λ(r) = 2λx/L, where x is the x-component of the position, r, r = (x, y).The factor of 2 ensures the average resources regeneration rate in the system to be equal to λ, and thus comparable with the homogeneous case.
In the supplementary information, I have also considered a two-patch environment or binary resource distribution.In the binary resource distribution case, the environment is divided into two patches, a poor patch, extending from x = 0 to x = L/2, with a constant resource regeneration rate λ 0 , and a rich patch, extending from x = L/2 to x = L with a resource regeneration rate λ.

Variable definition
Prey density, ρ p (predator density, ρ π ), is defined as the number density of preys (predators) defined as the total number of preys (predators) in the population divided by the area.In order to define prey and predator density fluctuation (δρ p and δρ π ) we first define the density field as the number of preys living on a lattice site (x, y), s p (x, y).The density fluctuation of preys is the standard deviation of this quantity δρ p = √ ⟨(s p (x, y) − ⟨s p (x, y)⟩) 2 ⟩, calculated over all the lattice sites.The density fluctuation of predators is defined in a similar way.
Prey and predator per area birth rate (b p and b π ) are defined as the average number of births per lattice site, derived by dividing the total number of births in a time step by lattice size.Per capita birth rates (b c p and b c π ) are defined as the average birth rate of a prey (predator) per time step.These are derived by dividing the per area birth rate by density.Prey and predator per area death rates (d p and d π ) are defined as the average number of deaths per lattice site.To calculate these, we divide the total number of deaths in a time step by lattice size.Per capita death rates (d c p and d c π ) are defined as the average death rate of a prey (predator) per time step.These are derived by dividing the per area birth rate by density.Prey per area production (P p ) is the number of preys captured by predators per lattice site.Prey per capita production, P c p is defined as per area prey production divided by prey density.
Group size of preys, s p and predator, s π is defined as the number of individuals in a lattice site.Maximum groups size is defined as the maximum group size observed in the stationary state of the simulation ran for T time steps (T = 20 000 in figure 5 in the main text and the fist 1000 time steps are discarded as non-stationary period).The lifespan of preys and predators is defined as the number of time steps that a prey or predator lives.Lifespan distributions are calculated based on a simulation run for T time steps and after the system reaches stationarity.
The temporal fluctuation of prey and predator density, birth rate, and spatial density fluctuation is defined as the standard deviation of these quantities in the stationary state.The fluctuation of prey and predator lifespan is the standard deviation of the lifespan distribution of all the prey or predators who born and died in the stationary state of the simulation.

Mean field approximation of the phase transition
To derive an approximate relation for the phase transition, note that assuming mass action law holds (which is a valid assumption in the regime of low densities), the average flow of resources into predators is equal to γqρ p , which should become larger than η at the transition point.Assuming that at the transition point, the prey's density remains close to its stationary value in the absence of predators, ρ p = λ/η, (based on the fact that predators' density is close to zero at the transition), we arrive at the naive mean field prediction, λ * = η 2 /γq.This prediction is plotted by a dashed black line in figure 2, which provides a lower bound for the estimated transition line from simulations in a structured population (blue circles), and is in high agreement with simulation results in a mixed population (black stars).The deviation from the mean-field prediction in a structured population can be understood by noting that predators deplete the prey density in their neighborhood.Consequently, a larger prey density than predicted by the naive mean filed argument is necessary to sustain predators in the system.

The shelter effect gives rise to sublinear growth
It is possible to derive an analytical form for prey production.Assuming there are s p preys and s π predators on a site, the expected number of preys caught by predators can be written as: where the theta function, Θ(x) is one if x > 0 and zero otherwise.To write down this expression, note that ( s π k ) (1 − q) sπ−k q k is the probability that k preys are caught.When k is smaller than the number of prey on the site, predators catch k prey.When k is larger than the number of preys, all the preys on the site are caught.Summing over k, we have the average number of preys caught on the site.Prey production is the average of K(s π , s p ) over all the sites: where, P(s π , s p ) is the probability of observing s p preys and s π predators on a site.It is possible to approximate this expression by a simple production function: This expression is exact when the number of prey on a site is larger than the number of predators, and an approximation is involved when the number of prey is smaller than the number of predators.As can be seen in figure 6, K (solid blue) exactly describes prey production, and K 0 (dashed red) shows a good agreement with simulation results.

The origin of scaling of averages
In the following I show the scale-invariance of the group size distribution together with a scaling of the cut-off with resource regeneration rate lead to scaling of the species populations.Consider the average density of species x in a given λ, ρ x (λ) = ´∞ 0 s x f(s x , λ)ds.(For discrete group size, s, we have ρ x (λ) = ∑ ∞ 0 s x f(s x , λ).One can use Euler-Maclaurin theorem to approximate this sum by an integral).Assuming the lower cut-off is zero, and using the simple scaling form we have ρ x (λ) = ´∞ 0 λ −β s x f 0 (λ −α s x )ds x .Using the fact that f 0 shows a sharp cut-off, the infinity in the upper bound of integral can be replaced by the cut-off, s F x = λ νx and assuming f 0 (s x ) has a power law form, f 0 (s x ) = s −δx x , we have: Thus, a λ,ρx = −β x + α x δ x + (2 − δ x )ν x .Using similar calculations, it is possible to shows that the exponent of the nth moment of population density of species x, ρ (n) x (λ) ≡ ⟨s n x ⟩ = ´∞ 0 x n f(x, λ), is equal to a λ,ρ (n) x = a λ,ρx + (n − 1)ν x .Finally, noting that in the stationary state the prey production should balance predator consumption, K(λ) = ηρ π (λ), we get the same exponent for prey production, a λ,K = a λ,ρπ = −β p + α p δ p + (2 − δ p )ν p .
Using the normalization of the probability density function, it is possible to derive approximate equalities which the exponents satisfy.Neglecting the mass of the probability density function above the upper cut-off and taking the lower cut-off equal to zero, we have ´sF 0 λ −β f 0 (sλ −α )ds = 1, from which one can derive −β x + α x δ x = −(1 − δ x )ν x .This implies a λ,ρx = ν x , which predicts a Taylor's power-law with an exponent of one, δρ x = (ρ x − ρ 2 x ) 1/2 ∼ ρ x ∼ λ νx , that is a highly aggregated dynamics observed in the non-equilibrium regime.
The best estimates of a λ,ρp and a λ,ρπ based on a log-log fit are respectively ≈1.19 and ≈0.91 (see the supplementary figure 10), and the best estimates for ν p and ν π are, respectively, 1.09 and 0.79 (figure 5).A linear fit in a log-log plot suggests δ p = 0.57 and δ π = 0.80 (figure 5).Upon inspection, α p = 1.1 and β p = 0.95 leads to a reasonable data collapse, which gives a λ,ρp = 1.23.For predators, α π = 0.9 and β π = 0.95, which gives a λ,ρπ = 0.72.The exponents continue to satisfy these equalities with an error of about 10% for other parameter values, which can be attributed to taking the lower cut-off equal to zero, and neglecting the mass of integral above a upper cut-off, as well as errors involved in a determination of exponents using a log-log plot.

Figure 1 .
Figure 1.The illustration of the model and its phenomenology.(a) Prey and predators move on a lattice.Preys use primary resources, regenerated by a rate λ per site, and predators capture prey (successfully with probability q) to reach an amount of γ of resources.Individuals pay a metabolic cost of η per time step, reproduce when their energy reaches a threshold d, and die out if they go out of energy.(b) The model predicts a top-heavy food web in small densities, where mass action law holds, a weak density dependence follows, predator-favoring life histories emerge, and the food web shows top-down control and higher resource leads to a higher density of predators.On the other hand, in large densities, the shelter effect gives rise to a scale-invariant ecosystem with a scale-invariant group size distribution and scaling laws, bottom-heavy food web, bottom-up control, and prey-favoring life histories.

Figure 2 .
Figure 2. Phase diagram and phase transitions in a predator-prey system.(a)The phase diagram of the model in the plane defined by the resource regeneration rate, λ, and resource gain of predators per prey, γ.For too low resource regeneration rates, only prey survive (marked by the letter (p).By increasing λ, the system shows a series of successive phase transitions to a phase where preys and predators coexist in equilibrium, E, a non-equilibrium phase marked by traveling waves of preys and predators, N, and a phase where the dynamics become unstable due to too large density fluctuations resulting in predators driving preys to extinction, U.The dashed line, λ * = η 2 /γq, where η is the resource consumption per time step, and q is the probability of successful predation, is a naive mean-field prediction for the phase transition above which predators survive, which is in good agreement with simulation results in a mixed population (black stars).(b) The densities of prey, ρp, and predators ρπ, as a function of λ for γ = 0.5.Different phases are shown in the system by letters.For too high resource regeneration rates, the system is unstable, and the time average of prey and predator densities in the presence of immigration shows a linear relation as a function of λ in a homogeneous environment.(c) The predator-prey transition, above which predators survive is a critical phase transition.Consequently, the reduced predator density, ρπ − ρ * π , reduced prey birth rate bp − b * p , reduced predator birth rate, bπ − b * π , and reduced prey production, Pp − P * p , show scaling as a function of reduced resource regeneration rate, λ − λ * close to the phase transition.Here, a superindex c indicates the value of the corresponding variable at the transition point.A reduced variable is defined as the value of the variable subtracted by its value at the transition.Prey production is defined as the number of prey captured by predators per time step.Parameter values: η = 0.05, d = 2, γ = 0.5, q = 0.5.A homogeneous environment is considered.Error bars are standard deviation (divided by 2 in (b) to increase visibility).

Figure 4 .
Figure 4. Fluctuation scaling laws.Taylor's power law can be extended to a rich set of fluctuation scaling laws, relating the fluctuations of densities and birth rates of species across the food web.(a) Temporal fluctuation (shown by σ) of prey (left) and predator density (right) as a function of their density shows Taylor's power law in time in the regime of large densities with a higher exponent for predators.(b) The temporal fluctuation of prey per area production (left) and predator per area birth rate (right) shows Taylor's power law in time with a higher exponent for prey production.In small densities, the trivial 0.5 exponent is observed.(c) The temporal fluctuations of prey per area production as a function of prey density (left) and predator per area birth rate as a function of predator density shows a power law, which can be considered a heterogeneous Taylor's power law as it relates fluctuations of different species characteristics.The birth rate fluctuations show a stronger density dependence for predators.(d) The temporal fluctuation of predator density as a function of the temporal fluctuations of prey density (left) and the temporal fluctuation of prey production as a function of the temporal fluctuation of prey density (right) shows a power law in the scale-invariant regime (large densities), which can be coined predator-prey fluctuation power laws.Similarly to the predator-prey power law, both power laws are sublinear.Parameter values: η = 0.05, d = 2, γ = 0.5, q = 0.5.Temporal fluctuations are defined as the standard deviation of the respective time series.

Figure
Figure Scaling of the spatial distribution.(a) and (b) The distribution of the number of prey, sp (a), and the number of predators, sπ, (b) on a lattice site show a power-law with a cut-off at large group sizes.Insets show data collapse for the re-scaled group size distributions.(c) and (d) The maximum number of prey, s F p , and predator, s F π observed on a lattice site shows scaling with λ at large densities.The scale invariant spatial distribution, together with the scaling of the maximum number of individuals observed on a lattice site, can explain scaling laws of population abundances.The population resides on a first nearest neighbor lattice with linear size L = 200.Parameter values: η = 0.05, d = 2, γ = 0.5, q = 0.5.

Figure 6 .
Figure 6.The shelter effect leads to a sublinear scaling.(a) K c = K/ρp, K c 0 = K0/ρp, per capita prey production P c p , prey birth rate b cp , and natural death rate, b c p − P c p , defined as the death rate of preys due to going out of resource, as a function of prey density.K exactly describes prey production.K and K0 are derived in the methods, equations (4) and (5), and are the analytical expression for prey production, and an approximation for prey production, respectively.(b) While K and K0 show scaling with prey density with the same exponent as prey production, the correlation function between prey and predator densities times q, q⟨spsπ⟩, which determines prey production if mass action law holds, deviates from prey production in the non-equilibrium regime.The population resides on a first nearest neighbor lattice with von Neumann connectivity and fixed boundaries, with linear size L = 200.Parameter values: η = 0.05, d = 2, γ = 0.5, q = 0.5.