Prey Patchiness, Predator Survival and Fish Recruitment

Bulletin of Mathematical Biology (2001) 63, 527–546 doi:10.1006/bulm.2001.0230 Available online at http://www.idealibrary.com on

Prey Patchiness, Predator Survival and Fish Recruitment JONATHAN WILLIAM PITCHFORD AND JOHN BRINDLEY∗ Department of Applied Mathematics, University of Leeds, Leeds LS2 9JT, England A simple mathematical model is presented for the population dynamics of fish larvae when the main food supply, in this case copepods, is spatially patchy in its distribution. Encounters by an individual predator (larva) with prey patches, and with individual prey within patches, are represented by Poisson processes. It is demonstrated analytically, and confirmed by numerical experiments, that prey patchiness fails to alter mean predator–prey encounter rates from their values for homogeneous prey distributions. Individual variance in encounter rate is, however, much affected. This has significant consequences for the (small) numbers of larvae surviving to metamorphosis and recruitment to the adult fish population. c 2001 Society for Mathematical Biology

1.

I NTRODUCTION

Oceans are not homogeneous, and the process of formulating mathematical models for marine populations inevitably involves, either explicitly or implicitly, averaging spatial inhomogeneities in some way. In the case studied here, we are concerned with modelling the growth of fish larvae; many fish species, notably the gadoids (cod, haddock) progress from the hatched egg through a larval stage before metamorphosis into juvenile and then adult fish. The larval stage is characterised by a huge depletion in numbers due to higher predation or food scarcity, the survival rate being typically O(1%) (Chambers and Trippel, 1997). The commercial importance of these species, and increasing global concern over the decline in exploited fish stocks, mean that it is crucial to develop an understanding of the key processes governing larval growth. To this end, much scientific effort has been focused on modelling the detailed interactions between larvae, their prey, and the environment [see, for example, Fiksen et al. (1998), Heath and Gallego (1998) and Scott et al. (1999)]. The fish larvae depend on zooplankton, principally copepods, for their food and zooplankton populations are well known to exhibit spatial patchiness at length scales ranging from tens of km to less than 10 m (Tsuda et al., 1993; Currie et al., 1998; Tokarev et al., 1998). Patches observed at a given spatial resolution may, on examination at a higher resolution, be seen to be patchy themselves so that there is ∗ Author to whom correspondence should be addressed. E-mail: [email protected]

0092-8240/01/030527 + 20

$35.00/0

c 2001 Society for Mathematical Biology

528

J. W. Pitchford and J. Brindley

a hierarchy of patches within patches (Pinel Alloul, 1995). Whilst the main driving force for zooplankton patchiness is probably physical, and can be effectively modelled using theories of turbulence and diffusion (Powell and Okubo, 1994; Abraham, 1998), there are also biological reasons why patchiness may be favourable for the organisms involved, including predator avoidance, prey exploitation and sexual reproduction (Davis et al., 1991; Folt and Burns, 1999). Since zooplankton patchiness exists at scales smaller than those used in most spatial simulation models of marine ecosystems, it is important to understand how patchiness at the sub-gridscale level might influence the biological processes at work in the ecosystem. The aim of the mathematical model presented here is not to explain zooplankton patchiness, but rather to investigate whether the observed small-scale patchiness can influence the growth of fish larvae which rely on zooplankton for food. It is worth remarking that, though this paper is concerned with fish larvae and zooplankton, its results will be relevant for any predator–prey system with similar encounter probabilities. Section 2 develops a basic model for individual predators encountering patches of prey, and encountering individual prey within the patches. The model assumes that encounters occur as Poisson processes, and estimates are obtained concerning the mean and variance of the number of predator–prey encounters per unit time. In Section 3 the model is applied to the case of fish larvae searching for zooplankton prey in a turbulent patchy environment, using existing theory to provide the necessary Poisson process parameters. Direct numerical simulations are presented to support the analysis, and it is shown that patchiness can cause substantial variation in the number of prey a larva encounters per day. That this variation can be crucial in determining how many fish larvae survive to maturity is demonstrated by individual-based numerical simulations in Section 4. Section 5 considers what features of larval development cause the system to be so sensitive to stochastic fluctuations, and shows that if a larva is aware of, and able to respond to, its environment then it can very effectively exploit patchiness. Finally, some definite conclusions are drawn, concerning the limitations of using deterministic rather than stochastic, and mean-field rather than individual-based, modelling approaches.

2.

A M ODEL FOR P REDATOR – PREY E NCOUNTERS IN A PATCHY E NVIRONMENT

The model for patch and prey encounters presented here was first proposed and developed in Beyer and Nielsen (1996). The central modelling assumption is that predators encounter prey patches, and encounter individual prey, as Poisson processes. Suppose that the prey are found only in disjoint identical patches, and that each individual predator not within a patch encounters patches according to a Poisson process (Pp) with rate parameter α. Thus a predator not inside a patch at time t has a probability α1t + O(1t 2 ) of encountering a patch in the small time interval

Prey Patchiness

529

[t, t + 1t), and its probability of encountering more than one patch in [t, t + 1t) is of O(1t 2 ). As in Beyer and Nielsen (1996), define V (<1) as the fraction of the whole volume taken up by prey patches, so that (1−V ) is the proportion of volume devoid of prey. First consider the situation where there is no interaction between predator and prey within patches. In this case, an initially homogeneous spatial predator distribution would remain homogeneous, since the predators are unaffected by the prey patches. Since the average rate of predators entering any given patch therefore equals the average rate at which predators leave it, defining β as the rate parameter of the Pp governing the event of a predator within a patch leaving that patch gives (1 − V )α = Vβ,

(1)

so that

1−V α. (2) V Now consider interactions between predator and prey within the patch. Assume that a predator, once inside a patch, encounters prey as a Pp with rate γ . If the predator encounters a prey, i.e., if a prey falls within the prey’s perceptive range, then it attempts to capture it. This process (in the case of larval fish this typically involves re-orientation, approach to prey, and a final attack) takes a finite amount of time, which will be assumed to be a constant τ in this model. After the interval τ the predator resumes its hunt for prey. While within a patch, however, a predator which does not encounter prey has a probability of leaving the patch, characterized by a Pp of rate β as above. In this way, events for a predator within a patch while an encounter is not taking place occur as a Pp of rate (γ +β), with any individual event γ having a probability of γ +β of being a predator–prey encounter, and probability of β of being an exit of the prey from the patch. Figure 1 summarizes this scenario. γ +β Events within the patch can be described using a binomial model: ‘success’ γ β probability γ +β corresponds to prey encounter, and ‘failure’ probability γ +β corresponds to leaving the patch. The mean number of prey encountered during one encounter with a patch is then M P , where β=

MP =

∞ X

nP (no. of prey encountered = n)

n=0

 ∞ X = n n=0

=

γ . β

γ γ +β

n

β γ +β (3)

The mean time for a predator to encounter a patch, encounter prey within it, and then exit the patch is TP , made up of the mean time to encounter a patch, the mean

530

J. W. Pitchford and J. Brindley

Figure 1. Schematic diagram of the model for predator prey encounters. The processes of the larva encountering a patch, encountering prey while within a patch, and leaving a patch, are all modelled as Poisson processes, with rate parameters α, γ and β respectively.

time taken to encounter and attempt to capture TP = =

γ β

prey before leaving the patch:

1 γ 1 + τ+ α β β β + α(γ τ + 1) . αβ

(4)

The mean encounter rate for an individual predator is therefore Rprey prey per unit time, where Rprey =

MP TP

=

γα β + α(γ τ + 1)

=

γV 1 + γVτ

(5)

[using (2)]. Note that (5) is independent of α, so the overall predator–prey encounter rate is independent of the rate at which patches are encountered in this model: the patchiness manifests itself only in its effect on the parameter V . Expression (5) shows that Rprey −→ γ V as τ −→ 0, as one would expect in the case where predator–prey encounters are instantaneous. Also Rprey −→ γ as V −→ 1, reproducing the classical Michaelis–Menten grazing function 1+γ τ with handling time τ for a homogeneous prey distribution.

Prey Patchiness

531

According to this simple model, (5) shows that the prey encounter rate depends only on the product γ V and the handling time τ . Indeed, if one defines c as the spatially averaged prey density and c as the prey density within patches, then one has the identity c = cV. Making a further reasonable assumption that the prey encounter rate within patches γ scales linearly with local prey concentration c, so that γ = ηc,

(6)

then implies that γ V = ηc. (The parameter η represents the rate at which volume is searched by a predator, independent of prey density or patchiness.) Then Rprey =

ηc . 1 + ηcτ

(7)

Equation (7) is precisely the encounter rate one would expect for a homogeneous prey distribution with prey density c and encounter duration τ . In other words, the mean encounter rate in a patchy environment is identical to that in a homogeneous environment. This result is applicable more generally. Consider, for example, a situation where prey within patches (of volume fraction V ) are at a concentration c1 while the prey concentration outside the patches is c2 , so that the spatially averaged prey concentration is given by c = V c1 + (1 − V )c2 . Applying exactly the same method as above, but this time allowing encounters both within and outside the patches, as Pp’s with rates ηc1 and ηc2 respectively, leads similarly to equation (7). The model can also be applied to ‘patches within patches’, with implications for the sampling of data in the field. Suppose a survey were to reveal the existence of prey in patches of a certain concentration and length scale. A second survey of the same prey population at the same time, but able to sample at a higher spatial resolution, might reveal that the prey within these patches were not homogeneously distributed, and arrive at a smaller patch length scale and higher prey concentration within the larger patches. If one assumes the above Pp model to hold, and follows exactly the same logic as above, then the estimated mean encounter rate would be that given by (7) in each case, dependent only upon the spatially averaged prey density.

532

J. W. Pitchford and J. Brindley

In exactly the same way, one can allow non-identical prey patches in the model. If some parameter p characterizes the distribution of patches, with α( p), β( p) and γ ( p) obeying the obvious generalizations of (2) and (6) for each patch of type p, then the above reasoning again shows that mean prey encounter rate is unaffected by patchiness. These results, as they stand, appear to lend theoretical credibility to simulation models where sub-gridscale patchiness in ignored, and a spatially averaged prey distribution is employed, thus casting doubt on the scientific worth of carrying out expensive in situ surveys of small-scale patchiness. The work contained in the rest of this paper seeks to show that any such conclusion about ignoring patchiness is misplaced, that patchiness introduces significant variability into the number of prey actually encountered by an individual prey, and that this extra variability can have a large effect on the population dynamics. We can demonstrate this as follows; let TD be the total time an individual predator spends foraging in one day (12 h in the simulations), and TF be the random variable representing the amount of time a predator spends within prey patches during the day. Clearly TF < TD with probability 1. Defining Z as the number of organisms consumed by a predator in a day, one has that the expectation of Z , E(Z ), is given by E(Z ) = TD Rprey . (8) Calculating the variance of Z , Var(Z ), is less straightforward since patch visits are not independent events: if a predator consumes a large amount of prey in the first patch it visits, it is likely to have less time to encounter prey in subsequent patches. However, Var(Z ) can be estimated in the limit as the handling time τ −→ 0 by conditioning on TF . Setting τ = 0, Z | TF is simply a Poisson random variable, so that E(Z | TF ) = Var(Z | TF ) = γ TF . (9) Similarly, setting τ = 0 also allows the process of entering and leaving patches to be regarded as an alternating renewal process with rates α and β. In this case, (Cox, 1962) concisely derives Var(TF ) for such processes, with the results E(TF ) = and Var(TF ) ≈

αTF α+β

(10)

2αβTD (α + β)3

(11)

as TD −→ ∞. Combining these results, Var(Z ) = Var(E(Z | TF )) + E(Var(Z | TF )) = Var(γ TF ) + E(γ TF )

Prey Patchiness

= γ 2 Var(TF ) + γ

533

αTF α+β

  αTD 2γβ ≈γ 1+ α+β (α + β)2   2γ (1 − V ) = γ TD V 1 + . (α + β)

(12)

An alternative approach using hyperexponential distributions, for a system in which τ is identically zero, is given in Beyer and Nielsen (1996). Equation (12) demonstrates that patchiness can have a major influence on the variance. It should be noted that the term γ TD V corresponds to the variance of a Poisson random variable with mean TD E with τ
= 0, i.e., to the spatially homogeneous situation, so (1−V ) that the factor 1 + 2γ(α+β) represents the increase in variance due to patchiness. Increasing the rate at which patches are encountered, α, all other parameters remaining equal, reduces the variance toward the homogeneous value as expected. Similarly, an increase in predator–prey encounter rate within patches, γ , increases the variance since more encounter events occur per patch visit. Finally, note that the product V (1 − V ) governs the size of the effect of patchiness on variance: the variance is highest where V = 12 so that prey patches fill half of the available space. Of course, the exact effect of patchiness on the variance in encounter rates depends on the calculation of α and γ , and the interdependence of these parameters, as seen in Section 3.

3.

N UMERICAL S IMULATION OF E NCOUNTER R ATES IN T URBULENT PATCHY E NVIRONMENTS

The numerical results presented in this section confirm the preceding analytical predictions that patchiness does not alter mean predator–prey encounter rates, but can alter the variance significantly, a phenomenon which will be of importance in the later simulations. The model used to estimate the encounter rates of predators with prey or prey patches (i.e., to calculate γ and α) is that developed by Rothschild and Osborn (1988) from the theory of Gerritsen and Strickler (1977), and discussed more thoroughly in Kjorboe (1997) and Sundby (1997). The model assumes that the predator swims at a constant speed v, encountering prey organisms which are also swimming at a constant speed u in random directions. An encounter occurs when the distance between predator and prey falls below the predator’s perceptive distance R. The model also attempts to take into account the effects of local turbulence in the fluid flow, by including w R , the root-mean-square turbulent velocity on the R length scale. The Rothschild and Osborn model gives, for a predator

534

J. W. Pitchford and J. Brindley

foraging in a homogeneous environment, γ = dπ R 2

u 2 + 3v 2 + 4w2R , 3(v 2 + w2R )1/2

(13)

where d is the number of prey organisms per unit volume. For the purposes of this paper the zooplankton prey are regarded as being essentially non-motile (u = 0), their motion being entirely governed by the local turbulent fluid speed w R . It is also assumed that the size of an individual prey item is small in comparison to R. If one regards the prey patches as coherent spheres of radius L, whose motion is determined by the local turbulence, then by the same logic expression (13) can be applied for the encounter of predators with prey patches, resulting in α = dπ L 2

3v 2 + 4w2L . 3(v 2 + w 2L )1/2

(14)

The necessary assumptions for (14) to hold, namely that R  L and that the only significant non turbulence-driven velocity is that of the predator, are reasonable in the case of fish larvae encountering zooplankton patches. The turbulent length scales a involved [a = R for predator–prey encounters, a = L for predator–patch encounters (Visser and Mackenzie, 1998)] are assumed to be greater than Kolmogorov scale. Thus the effect of fluid viscosity can be neglected, resulting in the simple expression for the relevant turbulent velocities wl = 1.9(a)1/3 ,

(15)

where  is the turbulent kinetic energy dissipation rate [see, Rothschild and Osborn (1988), Hill et al. (1992) and Sundby (1997) for details and justification]. The above model is best regarded as one of many possible methods of arriving at the parameters γ and α governing the Pp’s in question. Other models might produce different parameters, but the results would be qualitatively similar. The results in Fig. 2 serve to illustrate the validity of the the preceding analysis and to quantify the effect of patchiness, rather than being a thorough investigation of the particular encounter model chosen. The predators are taken to be haddock larvae, feeding upon zooplanktonic prey which may be thought of as the younger copepodite stages of Calanus (Cushing and Horwood, 1994). Results are displayed for small (S: length l = 0.44 cm), medium (M: l = 0.80 cm) and large (L: l = 1.20 cm) larvae encountering prey which are at an average density of 15 000 organisms m−3 , but at a density of 30 000 organisms m−3 within patches. The larval swimming velocity and perceptive distance are evaluated according to body length l: v = 1.5 l, and R = 0.75 l (Cushing and Horwood, 1994). The time taken for an encounter, τ , is held constant at 2 s. The turbulent dissipation rate  is set at 7.4 × 10−8 m2 s−3 , appropriate for turbulence generated at 20 m depth

Prey Patchiness

535

mean

std dev. 600

3000 L

2500

500

2000

400

1500

300

1000

M

500

L

200

M

100 S

S

0

0 0

2

4

6

8 10 12 14 16 18 20

0

2

4

6

8 10 12 14 16 18 20

Figure 2. Mean and standard deviation of encounter rate (encounters d−1 ) in a patchy environment, versus patch length scale (m). Results shown for small (S) medium (M) and large (L) haddock larvae. The curves represent estimates using equations (5) and (12), the points are the results of direct numerical simulaton (S ♦, M ×, L +).

by moderate wind velocities of 6 − 7 m s−1 (MacKenzie and Kjorboe, 1995). The horizontal axis represents patch radius, with a zero radius corresponding to a homogeneous prey distribution. Each solid curve shows estimates of the mean and standard deviation of the encounter rate calculated using (5) and (12). The analysis was validated by direct numerical simulation (S ♦, M ×, L +), each point being generated by simulating the encounters of 1000 individuals with prey during a 12 h day. Figure 2 shows that the analytical estimates for encounter rates work well, accurately predicting mean rates. The standard deviation prediction is also accurate for the small and medium larvae. The loss in accuracy for the large larvae is due to the fact that (12) assumes that encounters take up a negligible part of the 12 h day (i.e., Rprey TD τ ≈ 0). This assumption begins to lose validity for the largest larvae, which encounter up to 3000 prey per day. The most important point is that small-scale patchiness, with patch length scales less than 20 m, causes the standard deviation of the number of encounters per day to increase by an order of magnitude, compared to that in the homogeneous situation.

4.

S IMULATING L ARVAL G ROWTH IN A PATCHY E NVIRONMENT

The above results suggest that small-scale prey patchiness is unimportant for mean values of encounter rates. Its major effect on the variability of encounter rates, however, can affect the growth of individual predators by significantly broadening the probability distribution of the weight of an individual at a given time. In the case of larval fish, where only a very small proportion of hatched larvae contribute to recruitment, one is concerned with the tails of such a probability distribution and the effect of variability is amplified.

536

J. W. Pitchford and J. Brindley

Recruitment

0.035

det.

0.03

homog.

0.025

10 m 20 m

0.02 0.015

30 m

0.01 0.005 0.0 16

17

18

19

20 21 22 *103 Food supply

23

24

Figure 3. Graphs of recruitment, defined as the probability of a larva reaching metamorphosis, against food supply. Results are shown for the deterministic case, and for stochastic cases with homogeneous prey, and prey in 10, 20 and 30 m patches. Food supply is held constant for each simulation, with values ranging from 16 000 prey m−3 to 24 000 prey m−3 .

The model investigated numerically is for haddock larvae feeding on small copepods (Cushing and Horwood, 1994). An individual larva grows at a rate proportional to the amount of prey consumed, provided this rate of growth never exceeds some threshold (1.12 times its weight, in the following simulations). A larva grows until it reaches some pre-defined weight of metamorphosis, at which point it is regarded as contributing to recruitment. However, at every time step there is a significant probability that the larva will fall prey to external mortality, via predation by larger organisms for example, and a larva which suffers such a mortality event clearly does not contribute to recruitment. The amount of prey consumed in a day, and the probability of external mortality, are modelled stochastically, so that recruitment can be regarded as the probability of a newly hatched larva reaching metamorphosis. The details of the model are presented in the Appendix. Figure 3 shows the results of such a simulation, showing recruitment as a function of prey density, for the deterministic case, and for stochastic cases with homogeneous prey, and prey in 10, 20 and 30 m patches. The prey density in patches is double the spatially averaged value. Average prey density ranges from 16 000 to 24 000 prey m−3 . Since a larva can never increase its weight by a factor exceeding 1.12 d−1 , there is a minumum duration of growth of around 40 days before metamorphosis, resulting in a maximum recruitment probability of 3.3% due to high external predation during the growth period. The deterministic results show a very abrupt dependence of recruitment on prey density: recruitment is zero for 17 500 prey m−3 , but is near its maximum at 19 000 prey m−3 . That such an abrupt dependence is a feature of deterministic modelling methods, and probably not representative of the real situation, is shown by the results of the stochastic simula-

Prey Patchiness

537

tions. In the homogeneous stochastic case recruitment is broadly similar to that of the deterministic model, but with a slightly lower threshold below which there is no recruitment: 17 500 prey m−3 produces approximately 0.4% recruitment, but there is negligible recruitment at 17 000 prey m−3 . The patchy simulations again show recruitment at food levels where a deterministic model would not, but there is also a decrease in recruitment (relative to the deterministic model) for higher prey concentrations. The effect increases with increasing patch length scale, with recruitment in a 30 m patchy environment being less than half that predicted by the deterministic model at intermediate (20–22 000 prey m−3 ) prey densities. In summary, small-scale patchiness appears to smooth out significantly the sensitive dependence of recruitment on food supply. Reasons behind the positive and negative effects of stochasticity on recruitment are discussed in Section 5.

5.

D ISCUSSION

5.1. Explaining the role of stochasticity. The numerical simulations in Section 4 show that stochasticity, even if it does not affect mean rates, can profoundly effect the output of a model. Variability benefits the larvae at low prey densities, but has the opposite effect at high prey densities. A theoretical basis to explain this phenomenon is borrowed from the literature concerning adaptive behaviour in animals, as summarized in Houston and McNamara (1999). At the heart of the matter lies Jensen’s Inequality which states that, if f (x) is a positive increasing function with positive second derivative (i.e., it is convex), and X is a random variable with mean X , then E( f (X )) > f (X ).

(16)

In other words, the average of the function is greater than the function of the average. Conversely, if f (x) is concave, the inequality in (16) is reversed. The ) ∂2 f difference between E( f (X )) and f (X ) increases approximately as Var(X (X ) 2 ∂x2 (Real, 1980). In the context of animal behaviour (16) gives insight into why an animal may choose a random foraging strategy rather than a more deterministic one, even if the two strategies have an equal mean rate of food acquisition: if the animal’s reproductive value increases as a convex function of food, then the stochastic strategy is preferable. An analogous result has recently been proved in the context of stochastic differential equations for population growth (Alvarez, 2000): if a population X (t) has a deterministic growth rate which is a convex function of X but whose growth is also subject to white noise perturbations of some magnitude σ , then increasing σ will act to increase E(X (t)). No estimate of the magnitude of the effect is given in Alvarez (2000), and the result relies on the deterministic growth rate being everywhere convex. If the growth rate is linear then there is no mean effect on E(X (t)), while if the growth rate is concave then E(X (t)) decreases.

538

J. W. Pitchford and J. Brindley Growth rate

W0

Weight of larva

WT

Figure 4. Schematic diagram showing the dependence of growth rate on weight. The solid curve represents the deterministic growth rate, the dotted line represents a stochastic realization. Note that the stochastic value can never exceed the deterministic value in the upper, linear, part of the graph.

The fish larvae modelled in this paper have no choice but to develop in an inherently stochastic environment, the magnitude of this stochasticity having been estimated in Section 3. The factors governing the growth of an individual are highly nonlinear, and are detailed in the Appendix. The key points are that for each item of prey eaten, in addition to being one step closer to metamorphosis, a larva experiences • • • •

an increase in swimming speed and therefore in mean prey encounter rate, a decrease in external mortality, since it can more easily avoid predators, an increase in digestive efficiency, an increase in metabolic costs.

Only the last of these is detrimental to the larva’s development, the others all acting directly to increase its chance of surviving to recruitment. Thus the positive effects of consuming a prey can outweigh the negative ones by an amount which increases as the larva develops, resulting in a convex growth function. Indeed, for the model presented in the Appendix, one can show that a larva’s average growth rate is a convex function of its weight until some threshold is reached. Above this threshold the growth rate is a linear function of weight. This threshold weight is higher the lower the food supply; larvae in low food environments can expect to have a convex growth rate for a greater proportion of their lives. Figure 4 summarizes the situation schematically. The arguments of Houston and McNamara (1999) and Alvarez (2000) can therefore explain why, at low food levels, Fig. 3 shows recruitment in the stochastic models when there is none in a deterministic model: at low food levels the larval growth rate is convex during most of the larva’s development, and the larvae thereby benefit from stochasticity. The linear section of the graph in Fig. 4 is due to the fact that the larval growth rate saturates: a larva can never increase its weight by a factor exceeding 1.12 per

Prey Patchiness (a)

(b) *10 - 6

*

16

6

10 - 4

14

5

12

Recruitment

Recruitment

539

10 8 6 4

4 3 2

2 1

0 2

4

(c)

6

8 10 12 14 16 18 20 Patch length scale

0

2

4 6 8 10 12 14 16 18 20 Patch length scale

*10 - 3

Recruitment

12 10 8 6 4 2 0

2

4 6 8 10 12 14 16 18 20 Patch length scale

Figure 5. Recruitment against patch length scale (m). (a) 16 000 prey m−3 . (b) 17 000 prey m−3 . (c) 18 000 prey m−3 . Note the difference in orders of magnitude in the vertical axes.

day. This maximum (linear) growth rate is the key to understanding why stochasticity is detrimental to recruitment at higher food levels. When the encounter rate variance is large a larva has a high probability of encountering more prey than it can use for its growth: physiological limitations mean it cannot take full advantage of its good fortune. There is, however, a correspondingly high probability that it will encounter a number of prey significantly below the mean value. This will result in a very small increase, or even a decrease, in the individual’s weight. In the high mortality environment in which a larva finds itself, any temporary setback in its development significantly reduces its chances of reaching metamorphosis. In effect, the fact that the larval growth rate saturates at a certain level serves to decrease the mean prey encounter rate. The schematic representation of stochastic growth rates in Fig. 4 (dotted line) illustrates this: while a larva may benefit from stochasticity during its early (convex) growth phase, it will always fare less well than in a deterministic model during the linear growth phase. Figure 5, showing recruitment against patch length scale for constant food levels

540

J. W. Pitchford and J. Brindley

of (a) 16 000; (b) 17 000 and (c) 18 000 prey m−3 , consolidates the above arguments. At 16 000 prey m−3 there is an increasing trend in recruitment as patch size, and therefore stochasticity, increases. [The erratic nature of the graph is due to that the fact that recruitment is very low (10−5 ) at this food level.] In contrast at 18 000 prey m−3 , and at higher food levels, increased stochasticity causes decreased recruitment due to the saturation effect described above. Figure 5(b), at 17 000 prey m−3 , shows an intermediate situation. Recruitment is negligible in the homogeneous case, and rises to a maximum for patch length scales between 3 and 10 m. For larger patches the increased variability serves to reduce recruitment as in (c). In other words, at this particular concentration of prey, the combined effect of environmental stochasticity and physiological constraints on larval growth is to favour certain scales of prey patchiness. While bearing in mind the limitations discussed below, it is plausible to turn this argument around: larvae developing in an environment whose prey concentration and spatial distribution are predictable and fairly constant from year to year would evolve their swimming behaviour so as to maximize their recruitment probability. 5.2. Modelling individual behaviour. The results concerning means and variances of encounter rates presented in Section 2 are independent of the method used to estimate α, β and γ . While the relationship (2) is a natural way to define β when the predators are moving at random with no perception of their surroundings, other aspects of larval behaviour may lead to a different choice of the encounter parameters, and these alternative choices are easily incorporated into the growth model. That larval behaviour can have a profound impact on recruitment in a patchy environment is demonstrated in Fig. 6. Here the simulation is exactly that used to produce Fig. 3, except that the value of β has been reduced to 80% of the value given by (2). This models a situation where larvae, once within a prey patch, are aware of being in an advantageous environment and attempt to remain there. Of course, a larva would ideally wish to remain within a food patch indefinitely, but its difficulty in detecting whether or not it is within a patch, together with local fluid turbulence, would make this ideal impossible to realize. In many ways, reducing β by a small amount is similar to the approach to model larval growth in a patchy environment used in Davis et al. (1991). In that paper a random walk model is developed resulting in fish larvae swimming, on average, up a gradient of prey. The model in Davis et al. (1991) is deterministic in that each individual larva only encounters an average number of prey for its particular spatial location. However, the main conclusion of Davis et al. (1991), that recruitment is very sensitive to spatial patchiness, is strongly echoed by the research presented here. Figure 6 shows that, by reducing β by only a small amount, a larva can significantly improve its chance of recruitment. Patchiness is of particular benefit to such a larva at low food levels, where significant recruitment (of the order of 1%) is observed where a deterministic or homogeneous model would predict zero recruitment. Again, at higher food levels the larvae in a patchy environment fare less

Prey Patchiness

541

0.035 homog.

0.03 10 m 20 m

0.025 Recruitment

30 m

0.02 0.015 0.01 0.005 0.0 16

17

18

19

20

21

22

23

24

*103

Food supply

Figure 6. Results of simulations as per Fig. 3, but with a reduced value of β, to model the effect of larvae behaving so as to remain within food patches.

well than the deterministic or homogeneous models would predict, but they fare considerably better than identical larvae with a value of β set by (2) (compare with Fig. 3). 5.3. Limitations and improvements. The above model is, of course, rather idealized. Prey patches are assumed to exist as coherent entities on timescales longer than those associated with predator foraging, while the mechanisms acting to form the patches are assumed not to affect the predators. It is also assumed that the prey concentration within a patch does not alter due to predation by the larvae, i.e., that prey concentration is not significantly depleted on the predator foraging timescale. This involves an implicit assumption that the predator concentration be sufficiently small. In addition it is assumed that the prey organisms are all identical and constant in size. The above modelling simplifications are justifiable in that they focus attention upon the effect of small-scale patchiness, and stochasticity in general, on larval growth. In Pitchford and Brindley we address some of these limitations by incorporating more realism in the prey population dynamics, using the approach of Biktashev et al. (submitted). The conclusions, based on numerical simulations, are in agreement with those presented here, but two additional phenomena emerge. Firstly it can be beneficial for the prey to form high density patches, since predators in a patch cannot fully exploit the patch. Secondly the predators can also benefit from prey patchiness since reduced grazing pressure on the prey population allows it to grow, thereby providing food during the later stages of larval development. The patchy encounter model proposed in Section 2 has the strength that it is independent of how one determines the rates α, β and γ . The estimates in Section 3 rely on characterizing turbulence at the various appropriate length scales, but offer

542

J. W. Pitchford and J. Brindley

no mechanism for patch formation, nor explain why patchiness is not destroyed by turbulence. While there are certainly biological driving forces for the formation of patchiness (Folt and Burns, 1999), physical processes almost certainly dominate for the scales of organism and patches considered here (Tokarev et al., 1998). Further modelling work might seek to combine the turbulent processes driving patch formation with those governing predator–prey encounters. However, so long as the patch and encounter processes can be regarded as Pp’s, the results presented here are qualitatively robust. It is also implicitly assumed that every prey encountered is successfully captured. This is an over-simplification, since whilst increased turbulence increases encounter rates, it can have a detrimental effect on the probability of an encountered prey being captured (MacKenzie and Kjorboe, 2000), leading to an effective decrease in γ . Such effects could easily be incorporated into models, but are probably not important at the moderate levels of turbulence used in the simulations in Sections 3 and 4. One apparent weakness of the model is the constraint that predators cannot leave a patch during a prey encounter. This can be addressed as follows: suppose a predator can leave a patch during an encounter as a Pp with rate δ, where δ ≤ β. The probability of the predator leaving the patch during a single encounter is then 1 − e−δτ . Following the same method as that used to find (5), but with the added complication of having to allow for predators leaving patches either during encounters or while seeking prey, results in the expression Rprey =

αγ β + γ (1 −

e−δτ )

+ α(1 + γ τ )

.

The difference between the above rate and (5) is negligible, provided the time spent encountering prey is much less than the time spent searching for prey. This will be the case for most larval fish, whose ingestion of prey would seldom exceed 1000 organisms per day (43 200 s for 12 h foraging) even at times of abundant food. It is likely, as seen in some of the simulations in Section 4, that the predator would become satiated before the basic model loses its validity.

6.

C ONCLUSIONS

The principal conclusions from the model are as follows: • In the absence of ‘intelligent’ predator behaviour, small-scale prey patchiness does not significantly influence an individual predator’s average prey encounter rate. However, the variance in the encounter rate is very much increased by patchiness. • The increased variability caused by patchiness can have a positive or negative effect on the growth of an individual predator organism, compared to the

Prey Patchiness

543

growth of an identical organism in a deterministic environment. Whether patchiness is favourable to growth, or not, depends upon the detailed form of the organism’s growth function. • In the case of the growth of larval fish, the combination of high external mortality rates and physiological constraints on growth act to make the probability of an individual surviving to metamorphosis very sensitive to variability in food supply. In low-food environments increased stochasticity is beneficial to recruitment, but the effect is reversed in higher food environments. • If the predator is aware of being in a prey patch, and is able to behave so as to decrease its probability of leaving the patch, it is able to exploit the stochasticity and survive in much lower food environments than in the deterministic case. At a more general level, the work presented here shows that a simple mean-field deterministic approach to modelling might have misleading results. This is particularly relevant when one is modelling very rare events (such as the probability of a fish larva reaching maturity) since the tails of the probabilistic distributions involved can be greatly affected by environmental stochasticity. In such circumstances, a stochastic individual-based modelling approach seems advisable.

ACKNOWLEDGEMENTS We would like to thank Dr C C Taylor of the Department of Statistics, University of Leeds, Dr J W Horwood of CEFAS, Lowestoft, and Dr J E Beyer of the Danish Institute for Fisheries Research, Charlottenlund, for their valuable inputs to this research. The work was made possible by support for one of us (JWP) by NERC on grant GR3/11671.

A PPENDIX The mathematical model for larval growth used is essentially a stochastic version of that presented in Cushing and Horwood (1994) for the growth of haddock larvae. The key differences are that a constant mean food supply and a weight-dependent mortality rate are employed. Each individual larva hatches with weight W0 , and its subsequent weight W (t) evolves according to the equation dW = b(W )r (W ) − sW n . dt

(A1)

In (A1) the function r (W ) represents the rate at which food is consumed, b(W ) represents the efficiency of the conversion of consumed food into body weight,

544

J. W. Pitchford and J. Brindley

Table 1. Parameters used in the simulations of larval growth. Parameter (larval growth)

Value

Units

Parameter (to calculate Z )

Value

Units

s n rL wprey bmin bmax j W0 WT m1 m2

2.60 0.67 0.12 1.00 0.135 0.480 0.002 33.0 3165.0 0.0892 6.39 × 10−5

µg1−n — d−1 µg — — — µg µg µg (µg)−1

l A B v R  τ — — — —

see (A2) 2.01 × 10−3 0.2234 1.5 l 0.75 l 7.4 × 10−8 2.0 — — — —

m m(µg)−B — m m m2 s−3 s — — — —

and the term sW n models the individual’s metabolic rate. Explicitly,   r L W + sW n r (W ) = min , Z wprey b(W ) and b(W ) = bmax − (bmax − bmin )e− j W . The form of r (W ) deserves further explanation. The random variable Z represents the number of prey encountered in a day. Each prey has the same weight wprey , with wprey held constant at 1 µg to simulate predation on small copepods (Cushing and Horwood, 1994; Alcaraz, 1997), and the number of prey per unit volume is held constant for each simulation. [This approach differs, for pedagogical reasons, from that used in Cushing and Horwood (1994), where prey organisms increased in weight at a constant rate, but where prey concentration was depleted by predation]. The random quantity Z is calculated using the methodology described in Section 3, using the parameters therein, where the larval body length l and weight W are related according to l = AW B (A2) for positive constants A and B (Cushing and Horwood, 1994). Thus the rate of increase in larval weight is basically a linear function of the amount of prey consumed, but this growth rate is never allowed to exceed a maximum value of r L W . A larva continues to grow until it reaches a fixed weight of metamorphosis, WT (or until its weight falls below zero: a decrease in individual weight is theoretically possible, but is not a significant factor in the numerical simulations). However, at each one day time step, the larva of weight W (t) has a finite mortality probability Pmort , where −m 1

Pmort = e 1+m 2 W (t)

(A3)

Prey Patchiness

545

is a decreasing function of weight to represent the fact that larger larvae are slightly less prone to predation. The parameters are chosen to as to ensure that (A3) exactly coincides with the age-dependent formulation in Cushing and Horwood (1994) for a larva growing at its maximum rate r L W . The parameters used in the model are summarized in Table 1.

R EFERENCES Abraham, E. R. (1998). The generation of plankton patchiness by turbulent stirring. Nature 391, 577–580. Alcaraz, M. (1997). Copepods under turbulence: grazing, behaviour and metabolic rates. Sci. Marina 61 (Suppl. 1), 177–195. Alvarez, L. H. R. (2000). On the comparative static properties of the expected population density in the presence of stochastic fluctuations. J. Math. Biol. 40, 432–442. Beyer, J. E. and B. F. Nielsen (1996). Predator foraging in patchy environments: the interrupted poisson process (IPP) model unit. Dana 11, 65–130. Biktashev, V. N., J. Brindley and J. W. Horwood. Phytoplankton blooms and fish recruitment rate. J. Plankton Res., submitted. Chambers, R. C. and E. A. Trippel (1997). Early Life History and Recruitment in Fish Populations, London: Chapman and Hall. Cox, D. R. (1962). Renewal Theory, Monographs on Applied Probability and Statistics, London: Methuen. Currie, W. J. S., M. R. Claereboudt and J. C. Roff (1998). Gaps and patches in the ocean: a one-dimensional analysis of planktonic distributions. Mar. Ecol. Prog. Ser. 171, 15–21. Cushing, D. H. and J. W. Horwood (1994). The growth and death of fish larvae. J. Plankton Res. 16, 291–300. Davis, C. S., G. R. Flierl, P. H. Weibe and P. J. S. Franks (1991). Micropatchiness, turbulence and recruitment in plankton. J. Mar. Res. 49, 109–151. Fiksen, O., A. C. W. Utne, D. L. Aksnes, K. Eiane, J. V. Helvik and S. Sundby (1998). Modelling the influence of light, turbulence and ontogeny on ingestion rates of larval cod and herring. Fisheries Oceanography 7, 355–363. Folt, C. L. and C. W. Burns (1999). Biological drivers of zooplankton patchiness. Trends Ecol. Evol. 14, 300–305. Gerritsen, J. and J. R. Strickler (1977). Encounter probabilities and community structure in zooplankton: a mathematical model. J. Fish. Res. Bd. Can 34, 73–82. Heath, M. R. and A. Gallego (1998). Bio-physical modelling of the early life stages of haddock, melanogrammus aeglefinus, in the North Sea. Fisheries Oceanography 7, 110– 125. Hill, P. S., A. R. M. Nowell and P. A. Jumars (1992). Encounter rate by turbulent shear of particles similar in diameter to the Kolmogorov scale. J. Mar. Res. 50, 643–668. Houston, A. I. and J. M. McNamara (1999). Models of Adaptive Behaviour, Cambridge: Cambridge University Press.

546

J. W. Pitchford and J. Brindley

Kjorboe, T. (1997). Small-scale turbulence, marine snow formation, and planktivorous feeding. Sci. Marina 61 (Suppl. 1), 141–158. MacKenzie, B. R. and T. Kjorboe (1995). Encounter rates and swimming behaviour of pause-travel and cruise larval fish predators in calm and turbulent laboratory environments. Limnology and Oceanography 40, 1278–1289. MacKenzie, B. R. and T. Kjorboe (2000). Larval fish feeding and turbulence: A case for the downside. Limnology and Oceanography 45, 1–10. Pinel Alloul, B. (1995). Spatial heterogeneity as a multiscale characteristic of zooplankton community. Hydrobiologia 300/301, 17–42. Pitchford, J. W. and J. Brindley. Turbulence, the refuge effect, and fish recruitment. Math. Biosci., submitted. Powell, T. M. and A. Okubo (1994). Turbulence, diffusion and patchiness in the sea. Phil. Trans. R. Soc. London Ser. B 343, 11–18. Real, L. A. (1980). Fitness, uncertainty, and the role of diversification in evolution and behaviour. The American Naturalist 115, 623–638. Rothschild, B. J. and T. R. Osborn (1988). Small scale turbulence and plankton contact rates. J. Plankton Res. 10, 465–474. Scott, B., G. Marteinsdottir and P. Wright (1999). Potential effects of maternal factors on spawning stock-recruitment relationships under varying fishing pressure. Canadian J. Fisheries and Aquatic Sci. 56, 1882–1890. Sundby, S. (1997). Turbulence and icthyoplankton: influence on vertical distributions and encounter rates. Sci. Marina 61 (Suppl. 1), 159–176. Tokarev, Y. N., R. Williams and S. A. Piontkovski (1998). Small-scale plankton patchiness in the black sea euphotic layer. Hydrobiologia 375/376, 363–367. Tsuda, A., H. Sugisaki, T. Ishimaru, T. Saino and T. Sato (1993). White-noise-like distribution of the oceanic copepod neocalanus cristatus in the subarctic North Pacific. Mar. Ecol. Prog. Ser. 97, 39–46. Visser, A. W. and B. R. Mackenzie (1998). Turbulence-induced contact rates of plankton: the question of scale. Mar. Ecol. Prog. Ser. 166, 307–310. Received 11 November 2000 and accepted 22 January 2001