A scale-invariant model of marine population dynamics

# A scale-invariant model of marine population dynamics

José A. Capitán Grupo Interdisciplinar de Sistemas Complejos (GISC), Departamento de Matemáticas,
Gustav W. Delius Department of Mathematics, University of York, York, UK
###### Abstract

A striking feature of the marine ecosystem is the regularity in its size spectrum: the abundance of organisms as a function of their weight approximately follows a power law over almost ten orders of magnitude. We interpret this as evidence that the population dynamics in the ocean is approximately scale-invariant. We use this invariance in the construction and solution of a size-structured dynamical population model.

Starting from a Markov model encoding the basic processes of predation, reproduction, maintenance respiration and intrinsic mortality, we derive a partial integro-differential equation describing the dependence of abundance on weight and time. Our model represents an extension of the jump-growth model and hence also of earlier models based on the McKendrick–von Foerster equation. The model is scale-invariant provided the rate functions of the stochastic processes have certain scaling properties.

We determine the steady-state power law solution, whose exponent is determined by the relative scaling between the rates of the density-dependent processes (predation) and the rates of the density-independent processes (reproduction, maintenance, mortality). We study the stability of the steady-state against small perturbations and find that inclusion of maintenance respiration and reproduction in the model has a strong stabilising effect. Furthermore, the steady state is unstable against a change in the overall population density unless the reproduction rate exceeds a certain threshold.

###### pacs:
87.23.-n, 89.75.Da, 87.10.-e

## I Introduction

The population dynamics in the pelagic marine ecosystem (the open sea) are particularly amenable to mathematical modelling and analytic understanding. Because size is the most important factor in determining who eats who, rather than species Jennings2001 (), it is possible to work with a model in which a single function describes the total population density of organisms of weight at time , aggregated over all species. We can build on many previous works using partial integro-differential equations to describe the time evolution of this total population density, in particular SilvertPlatt1978 (); SilvertPlatt1980 (); Cushing1992 (); CamachoSole2001 (); BenoitRochet2004 (); Arino2004 (); Datta2008 (); Blanchard2009 ().

In this paper we exploit an additional special property of the pelagic zone: its approximate physical scale invariance. The open sea looks similar at a range of scales. That is due to the fact that over many orders of magnitude there are no physical features and no strong physical principles that would single out a particular intermediate scale. We expect that this approximate scale invariance of the environment breaks down only at small scales at which molecular diffusion processes dominate and at large scale at which geography affects ocean currents. This kind of scale-invariance is not present in terrestrial environments, where not only the physical structure sets a scale, but where in addition the effects of gravity quickly become important for larger organisms.

It is not a priori guaranteed that the scale invariance of the physical environment will also lead to scale invariance of the ecosystem. The organisms populating the environment and their interactions could break the scale invariance. This paper is however based on the assumption that evolution, in its drive to make use of all available ecological niches, has made optimal use of the scale invariant environment by filling it with an ecosystem that roughly preserves scale invariance. It is beyond the scope of the current paper to investigate the evolutionary mechanisms that would lead to the self-organisation of such a scale invariant ecosystem. We will be content with deriving the consequences of scale invariance for the dynamics of the population density.

The main observational evidence for approximate scale invariance is that the equilibrium size distribution of organisms in the open ocean is approximately given by a power law , valid over almost ten orders of magnitude Sheldon1972 (); Boudreau1992 (); KerrDickie2001 (); JenningsMackinson2003 (); Quinones2003 (); Marquet2005 (). Several theoretical models have been proposed in the past to derive this steady-state size spectrum PlattDenman1977 (); SilvertPlatt1978 (); SilvertPlatt1980 (); CamachoSole2001 (); BenoitRochet2004 (); Arino2004 (); Law2009 (); Datta2008 () but, as far as we know, in this paper we are the first to point out that this can be understood as a consequence of the scale invariance of the underlying dynamical model.

The main processes that affect the abundance of organisms as a function of weight are predation, reproduction, maintenance respiration and intrinsic death. We will in the next section write down a dynamical model that incorporates all these processes. In general such a model will not reproduce the observed power-law size-spectrum in the steady state. This is the reason why existing work either only models predation SilvertPlatt1978 (); SilvertPlatt1980 (); BenoitRochet2004 (); Arino2004 (); Datta2008 () or makes the simplifying assumption that reproduction, maintenance respiration and death are exactly proportional to predation CamachoSole2001 (). We are able to avoid those difficulties by exploiting scale invariance.

While scale invariance leads to a power-law steady-state solution, it does not guarantee that this steady state is a stable equilibrium solution. Indeed, a linear stability analysis of the jump-growth model Datta2010 (), which is also scale invariant, showed that, for realistic choices of the parameters, the steady state in that model is unstable against small perturbations. The steady state observed in nature is stable, so the pure jump-growth model is missing some important stabilising effect. This was the motivation for the investigations of the more general model in this paper. We will see that the inclusion of maintenance respiration and reproduction has a strong stabilising effect.

We do not model spatial variation in the population density. Nevertheless our model formally resembles a model in one space and one time dimension, where the space dimension represents the weight of the organisms. This formal analogy between space and weight will become particularly clear when we transform from the weight variable to the logarithmic weight variable in Section IV. A scale transformation in then corresponds to a translation in . In the special case where the scale transformations do not affect the time variable we end up with a model that is translationally invariant in both space and time, and many of the usual techniques can be applied, like for example the use of the standard Fourier transform in the solution of the linearised model in Section VI.

Another difference between our model and those usually studied by physicists derives from the fact that feeding is a non-local interaction in weight space. Most fish and plankton do not feed on other individuals that are close to their own weight. Instead they prefer prey that are substantially smaller. Similarly they produce offspring at a weight far below their own. This is the reason why the population density is modelled not by a partial differential equation but by the partial integro-differential equation (8) involving integral terms that encode the feeding behaviour (9) and reproductive behaviour (10). In the linear stability analysis this leads to the non-local dispersion relation (49).

The organisation of this paper is as follows. In Section II we formulate a stochastic model encoding feeding, reproduction, maintenance, and death. We do not make any assumptions about the rates for these processes but keep them as general parameters of the model. By taking the macroscopic limit of the stochastic model we derive the deterministic evolution equation (8). The derivation of a Fokker–Planck equation for the stochastic fluctuations away from the macroscopic model is left to Appendix A. In Section III we impose scale-invariance and derive sufficient scaling conditions (15) for the rate functions. Section IV is devoted to transforming the evolution equation into the more convenient form (27). We then use scale invariance to find a power-law steady state solution of the model in Section V, where we discuss as well its main properties and restrictions. The overall population level is predicted by the model and we show in Appendix B that is positive under biologically relevant assumptions. A time-damped power-law solution can also be found, as we show in Appendix C. In order to investigate the stability of the steady state we determine the spectrum of small perturbations around the steady state in Section VI. We summarise our results in the Section VII. In Appendix D we describe how variability in the absorption efficiency can be approximated by the addition of a diffusion term to the model.

## Ii Size-structured population model

We start by constructing a stochastic model for the dependence of abundance of individuals on weight and time, taking into account the basic processes of predation, reproduction, maintenance and intrinsic mortality.

As in previous work Datta2008 (), instead of keeping track of the weight of each individual, we aggregate individuals of similar weight into discretised weight brackets for . Weight is the only attribute of individuals that is used in this model. Species identity and life stage are ignored. The weight distribution of individuals in a large fixed volume is described by the vector whose entries give the number of organisms in each weight bracket. We will later let the size of these brackets go to zero to obtain the continuum model. Working with discrete weight brackets allows us to avoid the mathematical complexities involved when trying to describe a stochastic model directly in continuous space, see for example Kondratiev2008 ().

The primary stochastic processes involved in the model are illustrated in Figure 1, which shows the various ways in which an event can affect the number of organisms in bracket . Below we will show how the deterministic equation can be read off directly from this figure. In Appendix A we start from a master equation for the probability distribution and carry out the systematic expansion of van Kampen vanKampen () in powers of the inverse system volume . The same method was used in Datta2008 () to derive a jump-growth equation. This method is based on splitting each variable into a deterministic, macroscopic component describing the density of individuals in weight bracket , and a fluctuation component as

 ni(t)=Ωϕi(t)+Ω1/2ηi(t). (1)

The powers of volume are chosen so that the new variables and no longer scale with the system volume. This method not only gives the macroscopic behaviour for the densities at leading order in the expansion, at higher orders it describes the stochastic fluctuations around the macroscopic solution as well, giving at next-to-leading order a linear Fokker–Planck equation. However, because the system volume is so large, in this paper we will concentrate on the macroscopic, deterministic equation.

From Figure 1 we can obtain the contributions to the time evolution of from each of the processes we consider: predation (), reproduction ( for Birth), maintenance respiration ( for Respiration) and intrinsic mortality ( for Death),

 dϕidt=(dϕidt)P+(dϕidt)B+(dϕidt)R+(dϕidt)D. (2)

We will describe each of them in turn below.

### ii.1 Predation

A predation event moves a predator from a weight bracket before feeding to a higher weight bracket after feeding and removes a prey organism from a weight bracket . Let be the rate constant for such a predation event. As illustrated in Figure 1a, there are three ways in which a predation event can affect the number of individuals in a weight bracket : 1) an organism in bracket can eat another organism and grow; 2) an organism in bracket can be eaten by another organism; and 3) an individual in a lower bracket can absorb enough prey weight and grow into bracket . Because the rate at which a particular individual will encounter prey will be proportional to the density of prey, the probability that one of the individuals in bracket eats any of the individuals of bracket and increases its size to bracket in the time interval from to is . Hence the contribution of the predation events to the deterministic time evolution of the density of organisms in bracket is

 (dϕidt)P=∑j,k(−Pijkϕiϕj−Pjikϕjϕi+Pjkiϕjϕk). (3)

Our model for predation can be viewed as a generalisation of the Smoluchowski coagulation equation Smoluchowski1916 (); Drake1972 () which is obtained in the special case when the resulting weight is equal to the sum of the weights and of predator and prey.

### ii.2 Reproduction

Most fish reproduce by laying a large number of eggs that are subject to heavy predation. Only a small fraction of the eggs survives to hatch and join the consumer spectrum. In a size-spectrum model however, in which size is the only attribute of an organism and its life stage is ignored, all organisms are assumed to be prey and predator simultaneously from the moment they are born. We can therefore not model the egg life stage and therefore can not provide an entirely realistic model of the reproductive processes. Our model is thus designed to only capture two features of reproduction that we deem essential to the size-spectrum dynamics, namely that it moves biomass up the size spectrum from large weight to small weight and that it replenishes the population numbers at smaller weight.

We assume that a reproduction event moves a parent organism from a weight bracket to a lower weight bracket and produces a large number of smaller offspring in weight bracket . We set the number of offspring to , where denotes the largest integer smaller than , so that their combined weight is approximately equal to the weight lost by the parent. Let be the rate constant for such an event. The probability that one of the individuals in bracket reproduces in such an event in is then (note that reproduction, unlike predation, is a density-independent process).

As depicted in Figure 1b, there are three ways in which reproduction changes the number of the organisms in bracket : 1) the parent belongs to the -th bracket, and after spawning moves to a lower bracket; 2) a parent looses weight during spawning and moves to bracket ; and 3) the bracket receives the offspring from a reproduction event. The contribution of these three possibilities to the deterministic equation is therefore given by

 (dϕidt)B=∑j,k(−Bijkϕi+Bjikϕj+mjkiBjkiϕj). (4)

### ii.3 Maintenance and mortality

In between feeding events, organisms continuously draw upon their reserves to maintain themselves. The weight loss due to maintenance respiration is modelled by events that move individuals to the next-lower weight bracket, assuming that the width of each interval is small enough (this is not a restriction, because at the end we will take the continuum limit where all these widths tend to zero). The probability that any of the individuals in the bracket undergoes such an event in is , where is the maintenance respiration rate for the bracket . Figure 1c shows the two ways that these primary processes change the number of individuals in bracket , giving the contribution

 (dϕidt)R=Ri+1ϕi+1−Riϕi. (5)

Organisms can also die for reasons other than being eaten. This introduces a fourth process in the model that accounts for intrinsic mortality. With a probability , a single individual in bracket is removed from the system in (see Figure 1d). This gives the contribution

 (dϕidt)D=−Diϕi. (6)

Note that, like reproduction, both maintenance and intrinsic mortality are modelled as density-independent processes.

### ii.4 Continuum limit

We now take the continuum limit of each contribution to the macroscopic equation by writing and then taking the limit uniformly for all . The variables are combined into a density of individuals per unit weight per unit volume as a function of weight and time so that . The sum over weight brackets is replaced by an integral, . Continuum rate functions are introduced as

 P(wi,wj,wk)=Pijk/Δk,B(wi,wj,wk)=Bijk/(ΔjΔk),R(wi)=ΔiRi,D(wi)=Di. (7)

After taking the continuum limit, the ordinary differential equations (2) for the assemble into a partial differential equation for ,

 ∂ϕ(w,t)∂t=(∂ϕ(w,t)∂t)P+(∂ϕ(w,t)∂t)B+(∂ϕ(w,t)∂t)R+(∂ϕ(w,t)∂t)D, (8)

which will be the object of study in the remainder of the paper. The contribution from predation (3) in this limit becomes

 (∂ϕ(w,t)∂t)P=∫dw′∫dw′′(−[P(w,w′,w′′)+P(w′,w,w′′)]ϕ(w,t)ϕ(w′,t)+P(w′,w′′,w)ϕ(w′,t)ϕ(w′′,t)). (9)

Integrals run over all positive weights. Similarly, the continuum limit of (4) is

 (∂ϕ(w,t)∂t)B=∫dw′∫dw′′(−B(w,w′,w′′)ϕ(w,t)+B(w′,w,w′′)ϕ(w′,t)+w′−w′′wB(w′,w′′,w)ϕ(w′,t)). (10)

Maintenance and intrinsic mortality are described by

 (∂ϕ(w,t)∂t)R=∂∂w[R(w)ϕ(w,t)] (11)

and

 (∂ϕ(w,t)∂t)D=−D(w)ϕ(w,t), (12)

respectively.

We note that the model with only the predation term reduces to the jump-growth equation derived in Datta2008 () when the feeding rate is chosen as , for some feeding preference function and a fraction accounting for the feeding efficiency. As pointed out in Datta2008 (), when the typical predator:prey mass ratio is sufficiently large and the population density is sufficiently smooth that model in turn can be approximated by the McKendrick-von Foerster equation McKendrick1925 (); vonFoerster1959 (). It is this simplified model in terms of the McKendrick-von Foerster equation that forms the basis of most previous analytic studies of the size distribution SilvertPlatt1978 (); SilvertPlatt1980 (); Cushing1992 (); CamachoSole2001 (); BenoitRochet2004 (); Arino2004 (). We will see that, thanks to scale invariance, the much more general evolution equation (8) is amenable to similar analytical investigations.

## Iii Scale invariance

We now derive restrictions on the feeding rate , the reproductive rate , the maintenance rate and the mortality rate that guarantee the scale invariance of our model.

As we have a model describing the dependence of the abundance on two variables, namely weight and time, we can consider separate scaling in weight and in time. However, we expect our model to be invariant only under simultaneous scaling of both weight and time, because we expect life processes to run faster for smaller organisms. Hence, when we change the weight scale by a factor we should also change the time scale by a factor , where the constant expresses how the speed of the dynamics scales with weight. So we will consider the transformations

 (w,t)↦(cw,cξt) with c>0. (13)

Under such scale transformations the density transforms as

 ϕ(w,t)↦cγϕ(cw,cξt) (14)

with some, so far undetermined, exponent , also called the scaling dimension of DiFrancesco1997 ().

Requiring the evolution equation (8) to remain unchanged under this transformation imposes the conditions

 cγ−ξ−2P(w,w′,w′′)=P(cw,cw′,cw′′),c−ξ−2B(w,w′,w′′)=B(cw,cw′,cw′′),c1−ξR(w)=R(cw),c−ξD(w)=D(cw). (15)

We choose to factorise the feeding rate as

 P(w,w′,w′′)=S(w,w′)A(w,w′,w′′), (16)

where determines the rate at which a predator of weigh eats a prey of weight and gives the probability density that such a feeding event makes the predator grow to weight . On average only a certain proportion of the prey biomass will be absorbed by the predator, so that on average , but there will be some variability, due to differences in predator digestion and prey composition. This will be modelled by the probability density , so that

 ∫dw′′ w′′A(w,w′,w′′)=w+Qw′. (17)

We should make it clear that the model, and therefore the results, depend only on the combined rate . We factorise it into and only to make the ecological origin of the rate clearer BenoitRochet2004 (); Datta2008 (), but how we choose this factorisation has no influence on the results.

Because the model only depends on the product of and , we are free to choose the relative scaling between these factors, as long as the product scales as in (15). We want to scale as a probability density in and hence we choose

 A(w,w′,w′′)=w′′−1A0(w/w′,w′′/w),S(w,w′)=(w/w0)γ−ξ−1S0(w/w′), (18)

Here is an arbitrarily chosen reference weight. We have also introduced the scaling functions and that are invariant under scale transformations.

Similarly, we require that the reproduction rate scale as a density in and , so according to the scaling rules given in (15) we can write

 B(w,w′,w′′)=w′−1w′′−1(w/w0)−ξB0(w/w′,w′′/w) (19)

for some scaling function . This behaviour for the reproduction rate implies that the average weight of an offspring is proportional to , where is the weight of its parent. This should however not be taken as a prediction of the scaling of egg sizes because, as stressed previously, our model of reproduction is not formulated at a sufficient level of detail for that purpose.

The scale transformations for the maintenance and death rates given in (15) allow us to express them as

 R(w)=(w/w0)1−ξR(w0),D(w)=(w/w0)−ξD(w0), (20)

The functions , and , the constants and and the exponents and are not fixed by the requirement of scale invariance and need to be determined from observations or separate theoretical arguments.

The restrictions that we impose to achieve scale invariance predict the relative scaling of these rates. In particular, the scaling of all the rates contains the exponent . It is widely believed that the maintenance metabolic rate scales as ClarkeJohnston1999 (); Brown2004 (); Peters1986 (). Comparing that to the scaling in (20) gives . This then predicts that the mortality rate scales as , which is in line with observations Peters1986 (); Lorenzen1996 (). Moreover, the exponent of the scaling transformation (14) of the density is determined once the scaling of the feeding rate and are known. If we assume, following Benoît and Rochet BenoitRochet2004 (), that the feeding rate is proportional to the volume searched per unit time, we can identify the scaling exponent of as the exponent of the search volume, which is around 0.8 Ware1978 (). Assuming it follows that , which is in line with observations as well Blanchard2009 (). We will see in Section V that the exponent coincides with the exponent of the power-law steady-state solution. Therefore, the relative scaling between density-dependent processes (feeding) and density-independent processes (reproduction, maintenance and death) determines the steady state exponent.

Besides scale invariance, our model also has time translation invariance. Time translations and scale transformations together generate the non-abelian symmetry group , the group of orientation-preserving linear transformations of the real line, also known as the group. In the special case the symmetry group becomes abelian, a point which we will exploit in Section VI, because in that case the equation is translationally invariant and hence the standard Fourier transform can be applied to solve the equation for small perturbations of the steady-state.

The symmetry is enhanced when the model only contains the predation terms, which are quadratic in . In that case there is also an invariance under scaling in time alone,

 ϕ(w,t)↦λϕ(w,λt) (21)

for any .

## Iv Change to logarithmic weight

For the upcoming analysis it is convenient to make a change of variable to , where is the arbitrarily chosen reference weight. We refer to as the logweight. A scale transformation corresponds to a translation .

We introduce the density so that is the number of individuals in volume with a logweight between and . Thus and we can easily translate our results for back into results for , if desired. Under a scale transformation the function transforms to . We will apply this change to the various terms of the evolution equation (8).

We first introduce the predation rate such that . The factorisation of into and in equation (16) together with equation (18) leads to

 p(x,x′,x′′)=e(ρ−ξ)xs(x−x′)a(x−x′,x′′−x), (22)

where we have defined the functions and and we have introduced the exponent for latter convenience. We now transform to logweights, substitute this form for into (9) and perform a change of variables in each of the feeding terms so that the integration variable coincides with the argument of . The result is

 (∂u(x)∂t)P=e(ρ−ξ)x∫dy s(y)(−u(x)u(x−y)−e(ρ−ξ)yu(x)u(x+y)+∫dz e−(ρ−ξ)za(y,z)u(x−z)u(x−y−z)), (23)

where we have taken into account the fact that is a probability density and hence it is normalised. The integrals all run over the whole real line. In what follows, we will often not indicate the time-dependence of explicitly but write just instead, as in the above equation.

We can perform the same changes for the reproduction term. Using the scaling form (19) for in (10), then transforming to logweights and defining , we obtain

 (∂u(x)∂t)B=e−ξx∫dy∫dz b(y,z)(−u(x)+e−ξyu(x+y)+e(ξ−1)z(1−e−y)u(x−z)). (24)

Finally, the contribution of maintenance and intrinsic mortality can be expressed in logarithmic weights as

 (∂u(x)∂t)R=r∂∂x[e−ξxu(x)] (25)

and

 (∂u(x)∂t)D=−de−ξxu(x), (26)

respectively. Here we have introduced the constants and .

The full dynamical equation in terms of logweights is

 ∂u∂t=(∂u∂t)P+(∂u∂t)B+(∂u∂t)R+(∂u∂t)D. (27)

We have chosen not to fully non-dimensionalise this evolution equation: still has dimension of time, has the dimension of inverse volume, , and have dimension of inverse time, has dimension of volume over time and is dimensionless.

## V Power-law steady state solution

Solving the integro-differential evolution equation (27) is difficult in general. However we can simplify the task by looking for solutions that are invariant under symmetry transformations.

In this section we will study the solution that is invariant under both time-translations and scale transformations. We leave the study of the general scale-invariant solutions to Appendix C. Invariance under time-translations means that we are looking for a steady-state solution that has no dependence on time. Invariance under scale transformations (13) and (14) then implies

 ^ϕ(w)=^ϕ(w0)(ww0)−γ, (28)

where is the arbitrarily chosen reference weight. After transforming to logweights as in Section IV this takes the form

 ^u(x)=u0e−ρx, (29)

where and . Substituting this form for the solution into the evolution equation (27) gives an equation for the overall population level ,

 cPu0=−r(ρ+ξ)−d+cB, (30)

where the constants and are given by

 cP=∫dy s(y)(eρy+e−ξy−eρy∫dz e(ρ+ξ)za(y,z))cB=∫dy∫dz b(y,z)(−1+e−(ρ+ξ)y+e(ρ+ξ−1)z(1−e−y)). (31)

When this uniquely determines and hence the steady state solution. We will show in Appendix B that is positive under some ecologically reasonable assumptions about the parameter functions.

Note that scale invariance fixes the power-law form of the steady-state size spectrum and the steady-state exponent is determined entirely by the scaling behaviour of the parameter functions, see (15). It is not dependent on any other details of the interactions in the model.

A special situation arises in the case where maintenance, reproduction and intrinsic mortality are absent from the model. In this case only the first scaling relation in (15) remains and it is not enough to determine both and . However an equation for is obtained by noticing that in this case the right-hand side of (30) is zero and for this implies that . This constraint should then be used to determine the scaling exponent given a particular choice for and . The overall population level is not determined by the model in this case. This special situation was considered in most previous work SilvertPlatt1978 (); SilvertPlatt1980 (); BenoitRochet2004 (); Datta2008 ().

### v.1 Conservation of number of individuals

There is a continual flux of individuals from lower weight to larger weight to make up for the losses due to predation and intrinsic death. In previous models that considered only the predation process BenoitRochet2004 (); Datta2008 () there was no source for this influx of small individuals. Instead they appeared from . Now that we are modelling the reproduction process, we do have a source of individuals and can impose that in the steady state this source should exactly balance the losses.

The easiest way to impose this balance is to impose for each weight bracket that the number of individuals entering the bracket from the left due to predation exactly equals the number of individuals leaving that bracket to the left, either as offspring or through weight-loss. This gives

 ∑j,kPjki^ϕj^ϕk=∑j,k(mijk+1)Bijk^ϕi+Ri^ϕi. (32)

Thanks to scale invariance, all these conditions for different are equivalent. In the continuum, after substituting the steady state solution and changing to logweight notation, this condition reads

 fPu0=fB+r, (33)

where we have defined the constants

 fP=∫dy∫dzs(y)a(y,z)e(ρ+ξ)z+ρy,fB=∫dy∫dzb(y,z)(e−z(1−e−y)+1). (34)

We can use this constraint, together with the steady-state condition (30) to fix the maintenance rate in terms of the other parameters of the model,

 r=(cB−d)fP−cPfBcP+(ρ+ξ)fP. (35)

If we use this to eliminate from the equation (30) for we get

 u0=cB−d+(ρ+ξ)fBcP+(ρ+ξ)fP. (36)

Obviously, for the model to make sense we need to be positive. This is not possible for all choices of the other parameters. In particular, this requirement defines a range of allowed exponents . To investigate this further, we will now make specific choices for the functions that appear in the feeding and reproduction rates.

### v.2 Choice of parameter functions

For the prey selection function we choose a Gaussian that expresses that there is a preferred value for the log of the predator:prey mass ratio and a certain variance around this mean Ursin1971 (). So we set

 s(y)=s0gσβ(y−β) (37)

with

 gσ(x)=1√2πσe−x2/2σ2. (38)

The parameter has dimension of volume over time and sets the overall feeding rate.

For the absorption probability density the simplest assumption would be that a fixed proportion of prey mass is absorbed in all feeding events, i.e., that in terms of the predator mass and the prey mass the mass after feeding is always . This corresponds to a choice where

 ψ(y)=log(1+Qe−y). (39)

This was used in Datta2008 (). However the proportion of the prey mass that is absorbed by the predator is not exactly the same in each feeding event. Variability arises for example from the difference in digestion between predator species and also from the difference in organic composition of prey organisms. In this paper we will allow variation by replacing the delta function by a Gaussian. So we will set

 a(y,z)=gσψ(z−ψ(y)). (40)

For the reproduction function we will use the product of Gaussians

 b(y,z)=b0gσν(y−ν)gσμ(z−μ). (41)

This gives a mean offspring:parent mass ratio of and a mean mass ratio between parent before reproduction and parent after reproduction of .

Finding the correct values for the parameters requires a close investigation of the data and is outside the scope of this paper. However, for the purpose of the plots, we have chosen parameters that appear at least reasonable from a biological point of view. For example, the preferred predator:prey body mass ratio is believed to be around or JenningsMackinson2003 (), so we have chosen . In order to estimate , we need the average weight loss caused by reproduction processes. The average gonadosomatic index (ratio between the gonadal weight and the body weight) is actually measured for fish and is rather variable. It is found to be on average around 0.1 or 0.2 Roff1992 (), thus we have chosen so that the average fraction of weight loss due to reproduction is around 20%. We have set the logweight difference between offspring and parent to a small value around -8 or -9. For the standard deviations in the parameters we use , and , although a careful analysis of the data will be necessary to determine them properly. We have set the absorption efficiency to PandianMarian1985 () because respiration and other metabolic processes have been modelled separately. In previous work Datta2010 (), where these processes were not separated, the net absorption efficiency was replaced by a conversion efficiency of around . In most plots we will set the variability in equal to zero, as well as the mortality rate.

In Figure 2 we have plotted the maintenance rate and the steady-state density coefficient as functions of for the above choices of the parameters and . The allowed interval for appears shaded in that figure. It is encouraging that the observed value Sheldon1972 (); KerrDickie2001 (); JenningsMackinson2003 () is contained within the interval.

## Vi Stability of the steady state

It has been observed via numerical simulations in Law2009 (); Datta2008 () that the power-law steady state is not always stable against small perturbations but rather that the system can undergo a bifurcation in which the steady state becomes unstable and a stable travelling wave solution emerges. This phenomenon was investigated analytically in Datta2010 () through a linear stability analysis. We now perform a similar analysis in our generalised model.

The only other paper that we are aware of that investigates the stability of the power-law steady state is Arino2004 () but it deals, for reasons of simplicity, with a model where the growth due to feeding is independent of the prey density, thus avoiding having to deal with the associated non-linear terms.

In order to discuss stability analytically, we consider the particular case . According to (20), this corresponds to a maintenance rate proportional to the weight and a mortality rate independent of the weight. This is not quite realistic, but simplifies the analysis considerably because it leads to translational invariance in both time and logweight. This will allow us to use the standard Fourier transform.

### vi.1 Perturbation in u0

Before we consider the general, weight-dependent perturbation we take a look at a particular perturbation that affects only the overall population density . So instead of (29) we consider the solution

 u(x,t)=u0(t)e−ρx (42)

where we now allow to depend on time. Substituting this into the evolution equation (27) with gives

 du0dt=−cPu20+(cB−rρ−d)u0, (43)

where and are given in eq. (31). The solutions to this differential equation depend crucially on the strength of reproduction relative to maintenance and mortality. If reproduction is weak, i.e., if , then the non-zero fixed point (steady state) at is unstable. If, however, reproduction is strong enough so that

 cB>rρ+d, (44)

then the fixed point is stable. In between, the system undergoes a bifurcation at which there is a whole line of fixed points exactly when .

In particular, this observation shows that the model without reproduction can never be stable against a perturbation in the overall population density . This instability was already noticed in Datta2010 (), where it was argued that it could be avoided by suitable modifications of the model at the ends of the size-spectrum (plankton dynamics, senescent death). The observation that reproduction can provide a stabilising effect is new to this paper.

As we have seen, the stability against variations in requires that and are positive. In Figure 3 we plot these coefficients together with and as functions of , for two different values of . As discussed in Section V, the parameter is only allowed to lie in a certain interval where the maintenance rate is positive, which appears shaded. Within this interval, both and are seen to be positive, so steady state is stable in this case.

The region where both and are negative corresponds to an unstable steady state, as shown in Figure 4 with the same parameters but . For the following plots we will choose as a suitable value leading to a stable steady state (with ) and for the unstable solution we will choose (for ).

### vi.2 General perturbation

We now consider more general, weight-dependent perturbations. It is convenient to change to the new density related to through

 u(x,t)=u0e−ρxv(x,t), (45)

so that the steady-state solution is just a constant . We add a small perturbation to the steady-state solution

 v(x,t)=1+ϵ(x,t) (46)

and linearise the equation (27). Since the equation is linear and translationally invariant, we can solve it by means of the standard Fourier transform.

We can express any perturbation as a linear combination of plane waves labelled by a wavenumber ,

 ϵk(x,t)=ei(kx+ωt). (47)

In terms of that wavenumber, we get the following non-local dispersion relation

 iω(k)=−u0∫dy s(y)[e(ρ−ik)y+eiky−∫dzeρ(y+z)a(y,z)(eikz+e−ik(y+z)−1)]+∫dy∫dz b(y,z)[e−ρy(eiky−1)+e(ρ−1)z(1−e−y)(e−ikz−1)]+ikr, (48)

where we have used the steady-state condition (30) to eliminate the mortality rate .

The sign of determines stability. If is positive then the amplitude of the plane wave (47) with wavenumber grows exponentially with time, rendering the steady state unstable. We find that where

 λP(k)=u0∫dy s(y)(−(eρy+1)cos(ky)+∫dzeρ(y+z)a(y,z)(cos(kz)+cos(k(y+z))−1)) (49)

and

 λB(k)=∫dy∫dz b(y,z)(e−ρy(cos(ky)−1)+e(ρ−1)z(1−e−y)(cos(kz)−1)). (50)

Note that the maintenance rate parameter and the death rate parameter no longer appear in these expressions. Because a parent always uses weight during spawning, is positive wherever is nonzero, and hence we can see that is negative for any nonzero . This shows that reproduction always has a stabilising effect.

In the remainder of the section we will discuss the consequences on the stability of the steady state for the choices (37), (40) and (41) for the reproduction and predation functions. We show in Figure 5 the combined eigenvalue spectrum for two different values of , corresponding to both a stable (, ) and unstable (, ) fixed point . We find that the spectrum for is everywhere negative, corresponding to a stable steady state. We can see in the inset that the spectrum for is positive for small wavenumbers, leading to an instability of the steady state against very long-wavelength perturbations. At higher the spectrum is more negative, indicating stronger stability against short-wavelength perturbations.

In Figure 6 we show the contribution from predation for two different values of with and for two different values of for . Increasing from the value of used in Datta2010 () has a considerable stabilising effect. Note that we are allowed to increase because in this work we have separated out the losses due to maintenance processes. We can see as well that decreasing the preferred body size ratio between predator and prey has a stabilising effect. Nevertheless, for realistic values of the parameters, the contribution from predation alone is positive at some wavenumbers , showing that reproduction is required to achieve stability.

In order to characterise the effect of variability in we have considered non-zero in Figure 7. We have to keep the standard deviation sufficiently small so that the probability that the absorption efficiency is above is negligible. This leads us to impose that

 σψ≪(1−Q)e−y (51)

for some typical value of the log of the predator:prey mass ratio. In practice, typical values are and , so . Therefore the fact that the predator:prey mass ratio is so large implies that has to be very small, and to check its influence in stability we have chosen values around . The effect is negligible for small wavenumbers, although become slightly appreciable for highly oscillating plane waves. In Figure 7 we plot the difference , being the real part of the eigenvalue for . Although the effect is very small, the variability in the feeding efficiency always enhances the stability of the steady state. In Appendix D we show that the effect of these small fluctuations in the feeding efficiency consist on adding a diffusion term to the model.

We have also studied the contribution that reproduction makes to the eigenvalue spectrum for various values for and . The results are shown in Figure 8. As explained earlier, the spectra are always negative, showing that reproduction has a stabilising effect. Asymptotically converges to a negative constant. Variations in and affect the oscillatory behaviour found for small values of .

## Vii Conclusions

In this paper we made use of the fact that the power-law size spectrum that is observed in the pelagic ecosystem will be predicted by any dynamic model that is invariant under scale transformations. That allowed us to generalise earlier models without spoiling the prediction of a power-law steady state. Where earlier models only included the effects of predation and intrinsic mortality, we added terms modelling maintenance costs and reproduction and also allowed variability in the absorption efficiency. Inclusion of maintenance and reproduction has increased the stability of the steady-state solution.

We did not go into much ecological detail in this paper and made no attempt to determine the parameters of the model directly from ecological data. Nevertheless, by exploiting scale invariance, we made several observations that are of ecological relevance and that were not clearly made in previous work:

1. The power-law exponent for the size spectrum is fixed solely by the scaling properties of the parameters of the model. No detailed investigation of the model and its solutions is required to determine it. The exponent does not depend on details like the preferred predator:prey mass ratio, the feeding efficiency, the variability in feeding behaviour, the absorption efficiency, the maintenance costs, the mortality rate, or the details of reproduction.

This is in contrast to the results of earlier works in which only predation was considered. In that special case there are not enough scaling relations to fix the steady state exponent and it will depend on the details of the model.

It is a crucial aspect of our model that it contains both processes that are density-dependent (predation) and processes that are density-independent (maintenance respiration, intrinsic mortality, reproduction). It is the relative scaling of the rates for these processes that determines the steady-state power-law exponent.

Camacho and Solé CamachoSole2001 () studied the steady-state power-law exponent in a model with intrinsic mortality and reproduction. However they assumed that mortality and reproduction rates were proportional to predation rates. Thus, in effect, all their processes were assumed to be density-dependent and again the steady-state exponent was not determined by scaling arguments alone.

2. The assumption of scale invariance leads to predictions about the scaling behaviour (15) of the various parameters of the model, and these can be tested through observation, as discussed in Section III.

We can also make a prediction that has not yet been tested. The prey selection function can be related to data on the stomach contents of fish as follows. Let be the observed average number of prey with logweight between weight and found in the stomach of a fish of logweight . This stomach content reflects what prey a fish has been eating recently, which is determined by the same predation rates that we used in constructing our model. Thus we get (see equations (22) and (23))

 l(x,x′)=e(ρ−ξ)xs(x−x′)u(x′)T, (52)

where is a time related to the speed of digestion. We can assume that the fish density is close to the steady-state density (29). Substituting the steady-state density into the above equation gives

 l(x,x′)=u0e−ξxs(x−x′)T. (53)

So scale invariance makes a prediction for the form of the stomach content data. There is a lot of data available Barnes2008 () and it should be possible, via a careful analysis of this data, to determine the scaling exponent . We predict that this will confirm that .

3. From the condition that in the steady state the number of individuals must be conserved we deduced a relation between the parameters of the model which in particular restricts the exponent to an interval. We saw that for biologically reasonable values of the parameters this interval is around , which is in agreement with observations.

4. We have seen that the steady state in the model without reproduction, and with , is unstable against perturbations in the overall population density . In fact this was our motivation for including a reproduction term in the model in this paper. We have derived the condition (44) for the magnitude of the reproduction rate with respect to the magnitude of the maintenance and mortality rates which ensures stability against this perturbation.

5. We have studied the stability of the steady state against all small perturbations in the case . The stability is determined by the sign of , which has the two contributions from predation and from reproduction. Maintenance and mortality do not enter these expressions. The contribution coincides with that calculated in Datta2010 () in the case of fixed absorption efficiency, however with the conversion efficiency replaced by the much larger absorption efficiency . We have seen that this has an important stabilising effect. The contribution from reproduction is always negative, thus enhancing stability.

6. We have generalised previous models for the predation process to allow for variability in the absorption efficiency and have found that this does not have a big impact on the stability of the steady state. In Appendix D we show that this small effect can also be modelled by an additional diffusion term in the equations, which leads to simpler expressions.

In future we intend to extract further interesting information from the model by taking more detailed ecological information into account. In particular we intend to use theoretical results from metabolic theories about the energy budget Kooijman2000 (); Nisbet2000 (); Brown2004 () to determine the relative strength of the various processes in our model and observational results to choose appropriate values for parameters.

The analytic derivation of solutions and study of their stability that we performed in this paper should be pushed further. Numerical studies of the jump-growth model Law2009 () has hinted at the existence of travelling wave solutions. It would be nice to learn more about them within the framework of scale-invariance, possibly by extending our analysis in Appendix C. We would also like to investigate the speed at which perturbations move through the size spectrum.

Our stability analysis was restricted to the case . In this case the symmetry algebra , generated by time-translations and scale transformations, simplifies to the abelian symmetry group of translations in and . This allowed us to perform the stability analysis in terms of plane waves. It would be nice if the technique could be extended to the non-abelian case, possibly using the techniques from non-commutative harmonic analysis Kirillov1994 ().

The power-law size spectrum has been observed over many orders of magnitude and covers not only fish but also all types of plankton and even inanimate particles. Our model is appropriate only for organisms that feed by swallowing smaller organisms and that reproduces by spawning a large number of smaller organisms. It is not appropriate for phytoplankton or inanimate particles. Other models are needed for these and it will be interesting to see how these models can be coupled together. We propose that again the guiding principle should be scale invariance.

There is a limit to the amount of detail that can be incorporated into a community size spectrum model like the one described in this paper. In particular, because only the size of individuals is taken into account, their different behaviour in different life stages can not be modelled. A first step in the direction of refining the model was taken in AndersenBeyer2006 () where an individual is described not only by its size but also by one species-specific trait, namely its size at maturity. That allowed individuals to follow more realistic ontogenetic growth curves. An interesting observation of AndersenBeyer2006 () was how the size spectra for the different species, each of which singles out a particular scale in the form of the maturity size, combine into a power-law community spectrum described by a scale-invariant McKendrick–von Foerster equation.

In AndersenBeyer2006 () the size of offspring was encoded through a boundary condition on the species’ size spectrum but when combining these spectra into a community size spectrum the boundary condition had to be ignored in order to achieve a power-law steady-state size spectrum. It would be interesting to see if that work could be extended to include a realistic model of the reproduction process and still produce a power-law community spectrum.

In a pure size spectrum model like ours, that does not take life stages into account, a realistic modelling of reproduction is not possible. For example, the offspring of many marine species start their life in an egg stage during which they are not part of the consumer spectrum but are already preyed upon. In our model we had to neglect this fact. The scaling relation (19) implies that the typical size of offspring is proportional to the size of the parent. However, because we can not model the egg life stage, we can not claim that this gives a reliable prediction for the relation between egg size and parent size. Indeed, while such a relation may hold for copepods HuntleyLopez1992 (), it certainly does not hold for fish fishbase (), where most species lay eggs of a similar size.

One might be tempted to simply replace the reproduction function (19) by one that represents the fact that most fish lay eggs of a similar size. However this would immediately lead to a steady state solution that is not a power law. The steady state spectrum would exhibit a peak at the preferred egg size which is not observed in the size-spectrum data. The correct approach is not to break the scale invariance of the community spectrum model, but instead to study how more detailed models at the species level can combine into a scale-invariant model at the community level, in the spirit of AndersenBeyer2006 (); Hartvig2010 (). The species-level models would have to take into account effects like the duration of the egg stage, the separation between spawning grounds and feeding grounds, the various spawning strategies, like for example the laying of eggs in clusters, and many more. How and why the combination of many such detail-rich species-level models can lead to a scale-invariant community-level model is an intriguing problem and it is to be hoped that the experience that mathematical physics has with the emergence of scale invariance in complex systems can be exploited.

###### Acknowledgements.
We thank Richard Law and José Cuesta for useful comments and discussions. J.A.C. acknowledges funding by projects MOSAICO (Ministerio de Educación y Ciencia), MODELICO (Comunidad de Madrid), COST Action MP0801 (European Science Foundation) and by a contract from Comunidad de Madrid and Fondo Social Europeo.

## Appendix A Systematic expansion of the master equation

Our model is a Markov model and can be described by a master equation that gives the time evolution of the probability that the system is in the state at time . There will be a contribution from each of the processes involved,

 ∂P∂t=(∂P∂t)P+(∂P∂t)B+(∂P∂t)R+(∂P∂t)D. (54)

This appendix will be devoted to the systematic expansion of the master equation in powers of the inverse system volume vanKampen (). This expansion separates the macroscopic behaviour from the fluctuations, and gives a linear Fokker–Planck equation describing the fluctuations. Since the procedure is quite similar for the four processes we are considering, we will get into details for the predation process and will give the contributions of the other processes thereafter. An alternative derivation of the deterministic, macroscopic equation can be found in Datta2008 ().

A concise way of writing the contribution of predation events to the master equation uses the step operator notation vanKampen (). Since the probability to undergo a predation event in is , we have

 (55)

where the step operators act on any function as and .

In Section II we already introduced the split of each variable into a deterministic, macroscopic component describing the density of individuals in weight bracket , and a fluctuation component as

 ni(t)=Ωϕi(t)+Ω1/2ηi(t). (56)

The new stochastic variables have a probability distribution and

 (∂P(n,t)∂t)P=Ω−1/2∂Π(η,t)∂t−∑i∂Π(η,t)∂ηidϕidt, (57)

where we have used that (this follows from (1), taking the time derivative for fixed ). The step operator now transforms to , and can be expanded as

 Ei=I+Ω−1/2∂∂ηi+12Ω−1∂2∂η2i+⋯ (58)

Substituting this expansion into the master equation (55) we arrive at an equation containing different powers of the system volume . The highest order () terms in the expansion contain only macroscopic variables and vanish if these satisfy the macroscopic equation (3).

Terms at next order () give a linear Fokker-Planck equation for the probability distribution of the fluctuations,

 (∂Π∂t)P=−∑ijLPij∂∂ηi(ηjΠ)+12∑ijNPij∂2Π∂ηi∂ηj. (59)

By introducing the symmetric combination

 fijk=12(Pijk+Pjik), (60)

we can give concise expressions for the coefficients in the Fokker–Planck equation,

 LPii=−2∑jlfijlϕj,LPij=−2∑l(fijlϕi−fljiϕl),NPii=∑jl(2fijlϕiϕj+fjliϕjϕl),NPij=2∑l(fijlϕiϕj−filjϕiϕl−fljiϕlϕj). (61)

Terms at higher order in specify how the fluctuations deviate from being Gaussian. Fluctuations are seen to be damped by a factor , and the volume of our system is very large. This justifies that we can concentrate on the study of the macroscopic equation in this paper.

The remaining processes involved in the model can be treated in a similar fashion. The master equation for reproduction can be written as

 (∂P(n,t)∂t)B=∑i,j,kBijk(EiE−1jE−mijkk−I)[niP(n,t)]. (62)

At leading order the expansion in powers of gives the deterministic contribution (4). At next order the fluctuations are governed by a Fokker-Planck equation like (59) but with coefficients given by

 LBii=−∑jlBijl,LBij=∑l(Bijl+Bjlimjli),NBii=∑jl(Bijlϕi+Bjilϕj+Bjlim2jliϕl),NBij=∑l(−Bijlϕi−Biljϕi+Bljimljiϕl). (63)

The effect of maintenance on the time evolution of the probability is given by

 (∂P(n,t)∂t)R=∑iRi(E−1i−1Ei−I)[niP(n,t)], (64)

which leads to (5) at leading order and to the following contributions to the coefficients of the Fokker–Planck equatio