## Abstract

Many mathematical models use functions the value of which cannot exceed some physically or biologically imposed maximum value. A model can be described as ‘capped-rate’ when the rate of change of a variable cannot exceed a maximum value. This presents no problem when the models are deterministic but, in many applications, results from deterministic models are at best misleading. The need to account for stochasticity, both demographic and environmental, in models is therefore important but, as this paper shows, incorporating stochasticity into capped-rate models is not trivial. A method using queueing theory is presented, which allows randomness and spatial heterogeneity to be incorporated rigorously into capped rate models. The method is applied to the feeding and growth of fish larvae.

## 1. Introduction

Predators eat prey. When the predator is carnivorous the prey tends to arrive non-deterministically in discrete units. Many mathematical models which purport to describe this process use continuous approximations; typically an ordinary differential equation (ODE) is derived. When uncertainty plays an important role in the process, for example when prey are distributed patchily or the predator has a high risk of mortality, there is a concomitant need to extend these ODE models to incorporate stochasticity. This paper shows that some intuitively plausible approaches can mislead; the derivations of ODEs and their stochastic generalizations rely on infinitesimally small time steps, and are not necessarily applicable to biological systems.

This work springs from the authors' interest in models of the growth and recruitment of fish larvae into the adult population, a process which, for an individual larva, is fraught with peril. Larvae exist in a highly stochastic and patchy environment and possess only limited locomotory and sensory ability. Moreover, there are physiological limits on how fast an individual can grow, imposed by factors such as gut size and metabolism. Even under favourable conditions only a tiny minority of larvae survive to adulthood—the ‘average’ larva is most definitely dead.

Consider a general model describing an individual's change in mass, Δ*M*, during the time interval Δ*t*. Mass will increase according to the number of prey encountered during Δ*t*, but cannot exceed some physiological limit. Thus one may write(1.1)where *M* represents the mass of an individual, *Z*(*M*) is the rate at which prey is encountered and digested (dependent on local prey concentration), *f*_{1}(*M*) describes the efficiency of converting prey into mass gain, *f*_{2}(*M*) represents metabolic costs, and *g*(*M*) is the maximal rate at which the individual can grow. When food is scarce, growth occurs at a rate *f*_{1}(*M*)*Z*(*M*)−*f*_{2}(*M*) which depends (through *Z*(*M*)) on local prey concentration; when food is plentiful growth is ‘capped’ by a rate *g*(*M*) dependent only on *M*. Note that the mass of an individual is not subject to any constraint, only the rate of change of mass has an imposed maximum value.

It makes biological sense that the functions *f*_{1}(*M*), *f*_{2}(*M*) and *g*(*M*) are deterministic. If *Z* is also treated as a deterministic quantity, then allowing the time interval Δ*t* to tend towards zero in equation (1.1) leads to a limiting process described by an ODE. Particular cases, e.g. the Cushing–Horwood model (Cushing & Horwood 1994), have been explored previously (James *et al*. 2003). The main focus of these explorations was the time, *t*_{mat}, taken by an individual to reach a prescribed adult (or maturity) mass *M*_{mat}. This event of reaching a target mass and hence moving into a different life stage is usually described as recruitment.

Previous work (Pitchford & Brindley 2001; Pitchford *et al*. 2005) shows that it is imperative to treat *Z* as a random variable, dependent on the local environment in which the larva finds itself at any instant. In this case allowing Δ*t* to tend towards zero in equation (1.1) does not lead to a limiting process; the problem is ill-posed (see §2). Determining the distribution of *t*_{mat} requires the more detailed treatment based on queueing models described in §3. This approach is first illustrated with simple examples, where it is shown how spatial and temporal heterogeneity can be included, before being applied to the Cushing–Horwood herring recruitment model (Cushing & Horwood 1994). Section 4 contains the discussion and conclusions.

## 2. Two naive approaches to stochasticity

The following notation will be used to describe two important stochastic processes governing changes in an individual's mass, *M*. First, define *M*(*t*) as a diffusion process with constant drift(2.1)where *W*(*t*) is a standard Wiener process (Grimmett & Stirzaker 2001)—illustrated schematically in figure 1*a*. The process *W*(*t*) has instantaneous mean zero, so *M*(*t*) has expectation *λt*.

The second approach(2.2)uses a Poisson noise process (Feller 1950)—illustrated in figure 1*b*. Here *M*(*t*) is a compound Poisson process, increasing by one at each arrival of a Poisson process with mean rate *λ*. As in equation (2.1), *M*(*t*) has expectation *λt*. Either of *W*(*t*) and *N*_{λ}(*t*) could, in principle, be used to generalize equation (1.1). However, where capped rate models are concerned this can generate misleading results, as is described below.

### 2.1 Using *W*(*t*) to derive a stochastic differential equation

When the prey contact rate, *Z*, is defined as a random variable representing a heterogeneous food supply equation (1.1) is no longer an ODE; it becomes a stochastic differential equation (SDE). The field of SDEs is well established and there are many techniques, both analytical and computational, for dealing with such equations (see, for example, Higham 2001; Øksendal 2003). However, when the defining function for the rate of change is not smooth, as in the min(., .) function of equation (1.1), the authors are not aware of any suitable techniques. Indeed, for the case where *Z* follows a Gaussian distribution, it is argued that the problem is ill-posed.

Firstly consider the uncapped deterministic case(2.3)If *Z*(*M*) is no longer deterministic, the infinitesimal rate *Z*(*M*)d*t* should be replaced by *Z*(*M*)d*t*+*σ*_{z}(*M*)d*W*(*t*), where *Z*(*M*) is an average value and *σ*_{z}(*M*) measures the intensity of the stochastic perturbations (assumed to be driven by a Wiener process). The SDE(2.4)is obtained. This is a well-defined SDE that can be solved analytically or numerically depending on the forms of the various functions. The time to maturity in this case has been studied extensively, and some analytical solutions have been found (Karlin & Taylor 1981; Grimmett & Stirzaker 2001).

Now consider the deterministic capped-rate case(2.5)If *Z* is now stochastic, as above, it is tempting to write(2.6)and hope to arrive at a SDE in the limit as Δ*t*→0. However, no such limiting process exists, fundamentally because Wiener process sample paths do not have bounded variation (Grimmett & Stirzaker 2001). To illustrate this, consider a random trajectory as shown schematically in figure 2 where it has been assumed, for simplicity, that *g*(*M*)=*G* is constant. In order to implement equation (2.6) one must compare the gradient of the uncapped trajectory defined by the function with the gradient of the function *G*Δ*t* over a time interval Δ*t*. It is immediately obvious that the minimum of the two gradients depends on the chosen time interval Δ*t*; in (*a*) the chosen time interval results in an increase in mass only slightly below the maximum defined by the capped-rate. In (*b*) the same trajectory is split into three time intervals. Here the net increase in mass is much smaller. Numerical exploration using an Euler–Maryama scheme (Higham 2001) shows that, on average, as the time interval is decreased the change in mass decreases.

### 2.2 Using *N*_{λ}(*t*) to define a stochastic differential equation

Previous approaches (Beyer & Nielsen 1996; Pitchford & Brindley 2001) have remodelled the system as a Poisson process where *Z*(*M*) is the mean arrival rate. At each encounter the mass increases by an amount determined by the efficiency function *f*_{1}(*M*). The mass is also constantly decreasing according to the cost function *f*_{2}(*M*). Once again the uncapped system is well-posed,(2.7)When the metabolic costs are zero, *f*_{2}(*M*)=0, the time to maturity is simply the *N*th arrival of the Poisson process, where *N* is calculated from *M*_{mat} and *f*_{1}(*M*). In the case where *f*_{2}≠0, the problem can be couched as that of a random walk hitting a moving barrier. For a smooth growth function, i.e. an uncapped rate, it can be shown using a generalization of the central limit theorem that this formulation is equivalent to the SDE with , provided the number of food items consumed is large (Feller 1950; Whitt 2002).

Now consider the capped-rate system(2.8)A schematic is shown in figure 3. Again it is argued that no limiting process is reached as Δ*t*→0; the change in mass is very sensitive to the time interval chosen. If Δ*t* is of the same order as or less than the interval between arrivals of the Poisson process, then the increase in mass quickly tends to zero because at each arrival of the Poisson process the gradient is infinite.

## 3. Queueing models

A robust and generic method that can be used to derive an appropriate stochastic process for capped-rate growth is obtained by turning to queueing theory. Consider first the deterministic skeleton given by equation (1.1),In the regime where *M* increases at its maximum rate food is being encountered and digested at a rate ofA new random variable *X* is defined as the number of food items in the individual's gut. This variable is analogous to the number of items in a queue, in this case an M/D/1/N queue where M, the arrivals process, is a Poisson process (§2). The service process, D, is deterministic. The finite limit *N* imposed on *X*(*t*) represents the limited capacity of an individual's gut. As items in the queue are processed the individual's mass increases. The individual's mass will always increase at the capped physiological rate *g*(*M*) provided there is food available in the gut. This implies that, if food is available, it is processed by the gut at the maximum digestion rate *Z*(*M*)_{max}. In addition, the individual's mass will decrease according to the cost function *f*_{2}(*M*) when there is no food available.(3.1)Figure 4 shows a schematic of this process; see Whitt (2002) for a more detailed introduction and further examples.

### 3.1 The queueing model with a Poisson arrivals process

The basic behaviour of the queue model is explored here by considering a scenario in which the organism needs to consume a fixed number of prey, arriving at a constant rate (on average) and in the absence of metabolic cost. Suppose *f*_{1}(*M*)=1, i.e. for every item of food consumed the mass increases by one unit. Also set *f*_{2}(*M*)=0 so there are no metabolic costs. *g*(*M*)=1/*T*, so the time taken to consume a unit of food is *T* and food arrives with a constant mean rate *Z*(*M*)=*λ*. Also, the individual's gut has a finite capacity *N*. This constitutes the well-known M/D/1/N queue, for which results are available in the literature (Karlin & Taylor 1981; Brun & Garcia 2000). The behaviour of the queue falls into three well-defined regimes characterized by the traffic density *ρ*=*λ**T*,

light traffic

*ρ*<1 the service rate is higher than the arrivals rate, figure 5*a*;medium traffic

*ρ*∼1 the service and arrivals rates are of the same order, figure 5*b*;light traffic

*ρ*>1 the service rate is lower than the arrivals rate, figure 5*c*.

Figure 5*a–c* shows the probability distributions of the time taken to consume 1000 food items for each case. In the light traffic scenario, figure 5*a*, the system is dominated by the Poisson arrival process and as gut capacity increases the distribution tends to the equivalent Poisson process model (solid line). By *N*=10, the Poisson process solution and the queueing solution are indistinguishable. In the heavy traffic scenario, figure 5*c*, the system is dominated by the deterministic service process. As the gut capacity increases the distribution tends towards this limit. In between these is the medium traffic case, figure 5*b*, where the queueing model effectively comes into its own and has no equivalent model, either stochastic or deterministic, which it resembles. The probability density functions in all three figures were calculated numerically from 10^{6} simulations of the stochastic process.

It should be noted that, although the categorizations of traffic density here are defined using the notation from queueing theory, i.e. *ρ*=*λ**T*, it is clear from the original formulation, equation (3.1) thatThus one can rewrite the traffic density criteria asFor all three scenarios having a large gut capacity reduces the expected time to maturity, which is intuitively reasonable. Gut capacity can have a large influence on maturity times, and would therefore strongly influence recruitment. Moreover, the queueing model reveals that deterministic models, or interpretations based on capped-rate SDEs, would not necessarily capture the essential features of the process. In a low traffic regime, for sufficiently large gut size, the time to maturity is constrained by the arrivals process, with pdf symmetric about the deterministic average 1000/*λ*. Since the gut is seldom full (Brun & Garcia 2000), an uncapped SDE would capture the same result (Pitchford *et al*. 2005). However, for smaller gut sizes a full simulation of the queueing model is necessary. In a high traffic regime for sufficiently large gut size, the time to maturity is essentially deterministic (at 1000/*T*), and is constrained by the time taken to process food items. However, the queueing model again shows the importance of gut capacity; as capacity decreases both the expectation and variance of the time to maturity increase. Even at large gut capacity the medium traffic scenario has no realistic ODE or SDE analogue, but the same general trends (increased expectation and variance as capacity decreases) are clear.

### 3.2 A comparison with other capped-rate models

The concept of capped-rate models is not new in biology, a simple capped-rate model is the Holling type II (Murray 1989). Consider the following Poisson process based models

*A Holling II rate*. The (stochastic) rate of arrivals is governed by a smooth Holling type II function.(3.2)

*An unsmooth rate*. The arrivals rate is a linear function of mass when the mass is small and a constant for larger values. This gives an unsmooth approximation to the Holling II.(3.3)Notice that this unsmooth Poisson process formulation does not rely on a comparison of gradients as discussed in §2.

*A fully capped-rate*. This is the queueing theory model discussed earlier where the rate is never able to go above the cap.(3.4)For simplicity, metabolic costs have been excluded. For a more realistic approach, these costs would need to be incorporated into the models. Note also that the queueing theory model, in this case, has no constraint on the individual's gut size. Figure 6 shows schematics of these functions and their stochastic variations. Notice that for both of the standard Poisson process models it is possible via a ‘lucky’ random fluctuation for the growth to exceed the models capped-rate. However in the queueing theory model the in built cap is a definite one and the highest growth rate that can be achieved is the deterministic maximum. Observe that, in the case where there is assumed to be a deterministic physiological upper bound to the growth rate, neither the Holling model nor the unsmooth Poisson model correctly captures the underlying process which is successfully modelled by the queueing theory approach. If random fluctuations were to occur in the upper bound, perhaps due to internal physiological variability, the capped-rate model would require modification to include a stochastic rather than deterministic processing process but the underlying concept of a queue would remain valid.

Figure 7 shows the expected probability density functions of time to maturity for a realization of each of these models. In (*a*) the parameters *r* and *k* are chosen so that the growth is dominated by the linear, mass dependent, part of the function, i.e. the light traffic regime. In (*b*) *r* is increased to allow the growth to be dominated by the constant, mass independent rate, i.e. the heavy traffic regime. In both regimes the expected time to maturity is markedly different for each model and in the heavy traffic regime the pdf also has a decidedly different shape illustrating the effect of never allowing the growth rate to fluctuate above its deterministic maximum.

It is clear that the unsmooth rate and the capped rate are capturing different underlying dynamics. In the unsmooth case the organism's growth is determined by one function (a simple linear one here) at the start of its life. Then, when a certain size is reached (*k* in this case) a different rate underlies the dynamics. This is clearly a different scenario than that described by the queueing model.

### 3.3 A Pareto model

An advantage of the queueing theory model is that it is possible to explore the effects on the time to maturity of different types of arrival process. When a Poisson arrivals process is used the arrival times are exponentially distributed with a mean rate of arrival *λ*. A simple variation on this is to use a heavy-tailed distribution to allow the arrivals to have a clustered or patchy distribution. Distributions of this sort have been used frequently when modelling internet traffic which is well known to be ‘bursty’ (Fischer & Harris 1999). An example of this type of distribution is a Pareto distribution.

#### 3.3.1 The Pareto distribution

If prey items arrive as a Poisson process with mean rate *λ* per unit time, *X*, the total number of prey items that arrive in *U* units of time, follows a Poisson distribution with parameter *λU*, defined by the probability function(3.5)It follows that the time between arrivals, *T*, follows an exponential distribution with parameter *λ*, defined by the cumulative distribution function(3.6)The expectation and variance of the time between arrivals are given by(3.7)The Poisson process model can be extended to include the effect of patchy prey arrivals by allowing the mean prey arrival rate, *λ*, to be a random variable, while assuming that the conditional distribution of *X*|*λ* is Poisson (as in equation (3.5)). A plausible distribution for *λ* is the gamma distribution with parameters *r*, *α*>0, defined by the probability density function(3.8)where *Γ*(.) is the ordinary gamma function (see, for example, Abramovich & Stegun 1968). It can be shown that (see, for example, Johnson *et al*. 1992) *X* then follows a negative binomial distribution with parameters *r*, *α*>0, defined by the probability function(3.9)Thus the time between prey arrivals, *T*, has cumulative distribution function(3.10)and hence probability density function(3.11)which defines a Pareto distribution (see, for example, Johnson & Kotz 1970). It follows that(3.12)To quantify the effects of prey patchiness on recruitment, it is convenient to reparameterize the distribution by letting(3.13)to give(3.14)Comparing equations (3.7) and (3.14), observe that the expected time between prey arrivals is now the same as for the ordinary Poisson process model. The variance in the time between arrivals is controlled by the parameter *r*. As *r*↓2, Var(*T*)→∞, giving rise to a patchy arrivals process. As *r*→∞, Var(*T*)↓1/*λ*^{2} (i.e. the same as for the ordinary Poisson process model).

Figure 8*a–c* show the effect of a Pareto arrivals distribution on the different traffic regimes. For each figure the Pareto parameter is kept constant, *r*=2.1, and the mean arrivals rate is varied to change the traffic density.

In a light traffic regime, the introduction of a patchy arrivals process increases the variability in the time to maturity in comparison to the Poisson process, figure 5*a*. Moreover, the effects of reduced gut capacity are accentuated (i.e. expected time to maturity shifts dramatically to the right compared to a Poisson arrivals process with identical cut size). In a medium traffic regime, again, both the expected time to maturity and the variance in the time to maturity increase. The effects of reduced gut size are again amplified in comparison to a Poisson arrivals process. In a heavy traffic regime, for sufficiently large gut size, time to maturity tends towards the deterministic limit (as for the Poisson arrivals process), and again gut capacity is shown to be critical.

### 3.4 A model for fish recruitment

While simple models demonstrate the types of behaviour a system can display it is important to demonstrate the method's capabilities by using a more realistic model. The Cushing–Horwood model of larval fish growth (Cushing & Horwood 1994) is an ideal candidate for this approach. The model gives explicit functions for the capped growth rate, the larva's efficiency for converting food received into biomass, and for the maximum physiological growth rate. Initially, the external mortality rate in the model is ignored and it is assumed that all larvae whose body mass increases will reach maturity. This allows for easy examination of the probability density functions of the time to maturity.

#### 3.4.1 The Cushing–Horwood model

Cushing & Horwood (1994) presented a mathematical model for the growth of herring larvae. Each individual larva hatches with mass *M*_{0}, and then grows according to the equation(3.15)where *b*(*M*) is the larva's efficiency at converting food into biomass(3.16)and the metabolic costs and the capped physiological growth rate are(3.17)respectively. The mean rate of catching prey (in individuals s^{−1}) is found using the Rothschild–Osborne turbulent contact rates (Rothschild & Osborn 1988)(3.18)where *d* is the spatially averaged concentration of prey measured in individuals m^{−3}. Each individual prey item is assumed to have a mass of 1 μg, to simulate predation on small copepods. Larval body length and mass are related according to(3.19)and the larva's perceptive distance is(3.20)The appropriate turbulent fluid velocity for predator–prey encounters is(3.21)Each predator–prey encounter takes a fixed time *τ* (during this time no further encounters are possible) and larvae spend 12 h foraging each day. At this point the model can be used to explore the probability density functions of the time taken to reach maturity (cf. figure 10).

Finally, the larva is subject to a weight dependent external mortality probability(3.22)With this addition the model can be used to calculate the probability of reaching the target recruitment weight (cf. figure 11).

The parameters used are summarized in table 1.

Figure 9 shows schematically how the larval growth rate depends on mass in the Cushing–Horwood model. The model specifies that growth is governed by the minimum of the capped rate and the food dependent rate. The capped rate is necessarily independent of food concentration levels and is a simple linear function of mass. The food-limited rate is shown for three different prey concentration levels. In the low food concentration case the initial growth rate (i.e. at *M*_{0}) is governed by the food limited function and, furthermore, is negative. In this case the deterministic model will always predict that all larvae will die, but the stochastic model allows the ‘lucky few’ to experience an unusually high initial growth rate and hence there is a low but significant survival probability. In the medium prey concentration scenario a small part of the larva's initial time will be spent in a largely food limited environment but as the larva's mass increases it quickly moves into a regime that is dominated by the capped rate. Finally, in a high food environment the larva's growth will always be dominated by its capped rate as this is the minimum of the two functions for all values of mass.

Figure 10 shows the probability density functions for the time to maturity when the Cushing–Horwood model is implemented using the queueing theory model with Poisson arrivals and no external mortality. The pdf is calculated for a range of food concentration levels between *d*=17 000 and *d*=22 000 prey m^{−3}. At each concentration level three gut sizes are considered; 5, 10 and 100 prey items, respectively. It is clear that as the concentration level increases the system is dominated by the physiological cap and the variance of the time to maturity decreases rapidly as the gut size increases. Also, at the lower concentration levels, a larva with a small gut will never reach maturity.

Now, the external mortality function given in equation (3.22) is included and the focus is on the probability of reaching the target recruitment weight. Finally, figure 11 shows the probability of reaching maturity for varying food concentration levels and different gut sizes. Here an external death rate has been included in the model, so that rather than a pdf of time to reach maturity there is now a probability of reaching the given maturity weight. In figure 11*a* the arrivals process is Poisson, in figure 11*b* it is a Pareto process with *r*=2.1. In this more realistic model the modification to a Pareto arrivals process, to model a patchy food environment, has the effect of shifting the mean of the pdf of time to maturity to a higher level thus decreasing the probability of recruitment. This is the same effect as for the Pareto arrivals process in the simple model. For the Pareto arrivals process a gut size of 5 gives zero recruitment for the whole range of food concentration levels, whereas in the Poisson arrivals scenario this gut size gives significant recruitment levels at higher food concentrations. For the Poisson arrivals process the difference between a gut size of 10 and 100, while still significant is not large. However, in the Pareto case there a very large difference between *N*=10 and *N*=100. Even at *N*=20, the difference is still significant. It is interesting to note that for high concentration levels, when an individual resides in a world governed almost wholly by the capped physiological rate, the expected stochastic recruitment equals the deterministic recruitment. This is in direct contrast to previous results for uncapped growth which have predicted that increasing levels of stochasticity will always have a beneficial result on growth rates and recruitment (Pitchford *et al*. 2005).

## 4. Discussion

Although it both increases computational complexity and impedes significant analytical progress, the queueing approach offers several advantages. The generalization of the arrivals process quantitatively to account for spatial heterogeneity has been described above, but other factors such as variable prey size spectra and nutritional value, and variable predation, are readily incorporated. The fact that the key model parameters are individual-based rather than population-based makes the development of testable data-driven models more achievable.

The parameters used in §3.4 were not tuned so as to fit any particular observational data; the large temporal and spatial variations in observed recruitment would render this a meaningless exercise. The model's power is in explaining precisely which ecological processes drive recruitment, highlighting particularly the role of environmental stochasticity and mortality. However, all parameter values are based on published data (Cushing & Horwood 1994). It is interesting to note that larvae with the smallest gut capacity (5 prey items) fare considerably worse in medium traffic simulations where food is patchily distributed at realistic abundances. There would appear to be a strong evolutionary pressure to increase gut capacity from 5 toward 10, but a further (unrealistic) increase to 100 items appears to confer only a slight advantage. Refinements of the model which include a metabolic cost of growing and maintaining large gut capacity could further elucidate this evolutionary aspect, providing insight into how environment and risk should influence physiology; these are beyond the scope of this paper.

In the light of the results presented here it can be argued that ODE and SDE formulations of growth represent two extreme cases, with different traffic regimes in finite capacity queueing models providing a link between the two. For the M/D/1/N queue, as traffic increases and *ρ*→∞ the probability of the queue being empty tends to zero (Brun & Garcia 2000), so that the time to process a fixed number of items is essentially deterministic; ‘growth’ could then be described by an ODEAlternatively, if *N*→∞ and *ρ*≪1, i.e. light traffic and an infinitely long maximum queue length, then the probability density function for the number of items processed after time *t* is the same as that for the equivalent Poisson process. It can be shown, using the central limit theorem, that provided the number of items needed to reach maturity is large this can be approximated by an SDE (Karlin & Taylor 1981; Whitt 2002). Outside these extremes, where growth is stochastic and arrival and processing operate on similar time scales, it is impossible to provide a sound mathematical description without recourse to the individual-based stochastic process (i.e. the queue) governing the dynamics.

SDE formulations remain an important and powerful tool in understanding growth and population dynamics; the message contained herein is that they ought to be used with caution wherever the underlying deterministic process is not smooth. In such cases an iterative halving of time step and master equation transition probabilities may not lead to a meaningful convergent solution. If the existence of a ‘natural’ time scale can be argued on physiological, ecological or physical grounds, then SDE approaches can yield compelling results even when the underlying system is not smooth (see, for example, Truscott & Gilligan 2003, where the authors model the stochastic development of an epidemic under the influence of temperature as an SDE with a fixed time step corresponding to a characteristic time scale for soil heating). However, where an arbitrary time step is artificially imposed, perhaps through simple historical precedent (Cushing & Horwood 1994; Pitchford & Brindley 2001), a model's outputs may have an unrealistic dependence on this choice of time step; simply decreasing the time step does not solve this problem.

While being based on the larval growth models outlined previously, the queueing approach presented here could apply to any scenario where an organism (or cell, or molecule) with a finite capacity has to experience a certain number of similar stochastically driven events, each of which takes a finite time, before exhibiting a response. For example, ongoing work seeks to apply the method to stochastic ‘integrate and fire’ models for nerve impulses (Iyengar 2000).

## Acknowledgments

We would like to thank Martin Hairer for useful advice regarding the existence of limiting stochastic processes and the anonymous referees for their thorough reviews and their clear and constructive comments.

## Footnotes

- Received February 28, 2005.
- Accepted July 8, 2005.

- © 2005 The Royal Society